Validation of the Critical Flow Modelsin POLCA-T

DEGREE PROJECT, IN NUCLEAR REACTOR TECHNOLOGY , SECOND LEVEL

STOCKHOLM, SWEDEN 2015

Validation of the Critical Flow Models

in POLCA-T

JONATHAN WÄNG

Master of Science Thesis

Validation of the Critical Flow Models in

POLCA-T

Jonathan W¨

ang

Nuclear Reactor Technology, Department of Physics, School of Engineering Sciences

Royal Institute of Technology, SE-106 91 Stockholm, Sweden

Typeset in LA_TEX

Examensarbete inom ¨amnet fysik f¨or avl¨aggande av civilingenj¨orsexamen inom utbildnings-programmet Teknisk Fysik.

Graduation thesis on the subject of Physics for the degree of Master of Science (MSc) in Engineering from the School of Engineering Physics.

TRITA-FYS-2015:38 ISSN 0280-316X

ISRN KTH/FYS/--15:38-SE

c

Abstract

A theoretical study of the critical flow phenomenon is performed and the critical flow models that are analysed in the thesis are introduced and discussed in terms of their strengths and weaknesses. The models are the Homogeneous Equilibrium Model (HEM), Moody, and the Frozen Flow Model. The validation process of nuclear power safety computer codes is explained and the system code POLCA-T is described. We also discuss the Marviken Critical Flow Tests that provides the experimental data that the thesis uses as a validation basis. Two methods to estimate the uncertainty in the models are also discussed.

Next, we discuss the geometrical model and the validation matrix. Each studied case is described in detail. The sensitivity analyses of two parameters are also explained. These parameters are the initial temperature profile and the geometrical nodalisation. The results and the estimated uncertainty from each transient are then presented.

Finally, the results are discussed. It is found that the Moody model gives conservative results in the two-phase flow region in all cases except when the exhaust pipe diameter is small. HEM generally gives accurate results throughout most transients and has an esti-mated uncertainty in the flow rate calculations of between 6% and 10%. The uncertainty increases further if the exhaust pipe diameter or the exhaust pipe length is small. Addi-tionally, larger uncertainties are found in studied cases where the fluid is heavily subcooled or where the pressure is lower. HEM is observed to be sensitive to the initial temperature profile and many of the uncertainties can be decreased if the temperature in the vessel is accurately known. The Frozen Flow model shows high capability to calculate single phase flow, which is the model’s main area of implementation. The nodalisation only has a minor impact on HEM and Moody and simplification of the geometry is therefore possible in order to save computational time.

Key words: Critical Flow, POLCA-T, validation, Homogeneous Equilibrium Model, Moody, Frozen Flow, Marviken CFT

Sammanfattning

En teoretisk studie av fenomenet kritisk str¨omning ges, och de analyserade modellerna som beskriver kritisk str¨omning introduceras och diskuteras utifr˚an deras styrkor och svag-heter. Dessa modeller ¨ar Homogeneous Equilibrium Modellen (HEM), Moody, och Frozen Flow Modellen. Vi diskuterar ¨aven de experiment f¨or kritisk str¨omning fr˚an Marviken som anv¨ands som experimentell data i valideringen. Tv˚a metoder f¨or att uppskatta os¨akerheten hos modellerna diskuteras ocks˚a.

D¨arefter diskuterar vi den geometriska modellen och valideringsmatrisen. Varje fall som analyseras f¨orklaras i detalj. K¨anslighetsanalyserna som utf¨orts f¨or tv˚a parametrar f¨orklaras ocks˚a. Dessa parametrar ¨ar den initiella temperaturprofilen och den geometriska nodalise-ringen. Resultaten och den uppskattade os¨akerheten fr˚an varje transient presenteras sedan. Slutligen diskuteras resultaten. Vi har funnit att Moody ger konservativa resultat f¨or tv˚afasfl¨ode i alla analyserade fall f¨orutom n¨ar diametern p˚a utfl¨odesr¨oret ¨ar liten. HEM ¨

overensst¨ammer generellt v¨aldigt bra med experimentdata och os¨akerheten i fl¨odesber¨akningen uppskattas till mellan 6% och 10%. Os¨akerheten ¨okar ytterligare om utfl¨odesr¨ordiametern eller utfl¨odesr¨orl¨angden ¨ar liten. Dessutom observeras st¨orre os¨akerheter om vattnet ¨ar kraftigt underkylt eller om trycket ¨ar l¨agre. Vi har funnit att HEM ¨ar v¨aldigt k¨anslig

iv Sammanfattning

f¨or den initiella temperaturprofilen och m˚anga av os¨akerheterna kan minskas om tempera-turf¨ordelningen i tanken ¨ar v¨alk¨and. Frozen Flow modellen visar h¨og kapacitet f¨or enfas-fl¨odesber¨akningar, vilket ¨ar modellens huvudsakliga implementeringsomr˚ade. Nodalisering-en har Nodalisering-endast Nodalisering-en mycket litNodalisering-en inverkan p˚a b˚ade HEM och Moody och geometrin kan d¨arf¨or f¨orenklas f¨or att spara ber¨akningstid.

Preface

This thesis is the product of my work at the Department of Physics at the Royal Institute of Technology (KTH), in collaboration with Westinghouse Electric Sweden AB, during the period of January 2015 to June 2015. The topic of this thesis is thermal hydraulics and in particular critical flow. A few critical flow models are introduced and the implementation of these models in the coupled code POLCA-T is validated against experimental data from the Marviken Full Scale Critical Flow Tests.

Overview of the thesis

This thesis is structured as follows: Chapter 1 gives an introduction to nuclear reactor technology and nuclear power safety in general. Nuclear power computer codes are discussed with the main focus on POLCA-T. The validation procedure of such codes is described and the Full Scale Marviken Critical Flow tests are discussed.

In Chapter 2 the theoretical foundation of the thesis is laid. The critical flow phe-nomenon is explained and the three critical flow models that are implemented in POLCA-T are discussed. In addition, two uncertainty estimation models are described.

In Chapter 3 the methodology of the validation procedure is described. The geometrical model is explained and the experiment matrix is studied in detail.

The results of the simulations are shown in Chapter 4. These results are then discussed separately for each case in Chapter 5.

Finally, Chapter 6 concludes the thesis. The thesis is summarized and conclusions are drawn based on the observations made.

Acknowledgements

First of all, I would like to thank my supervisor at KTH, Jan Dufek, for helping with all administrative parts and for setting me up with the contacts that enabled me to do my thesis at Westinghouse. Secondly, to my supervisors at Westinghouse, Milan Tesinsky and Ulf Bredolt, for providing help and a guiding hand throughout my work in addition to proof reading my thesis. I have learned so much thanks to you. The rest of the people at SET have my deepest gratitude for welcoming me to the group and for answering my questions whenever I had any.

I would like to extend a huge thank you to my family and friends for all the love and support and for taking my mind of things in times of stress. My time at KTH would not have been the best time of my life, were it not for you. Finally, I would like to extend the very largest of thanks to Erika for always supporting me in all my endeavours and for putting up with me for all these years. None of this would have been possible without you.

Contents

1.2 Nuclear Power Safety . . . 2

1.2.1 LOCA . . . 3

1.2.2 Computer Codes . . . 5

1.3 Marviken . . . 8

1.3.1 Marviken Nuclear Power Station . . . 8

1.3.2 Marviken Full Scale CFT . . . 9

1.4 Objective . . . 11

1.5 Scope . . . 11

2 Theoretical Backgound 12 2.1 Critical Flow . . . 12

2.1.1 Single Phase Flow . . . 12

2.1.2 Two-Phase Flow . . . 14

2.1.3 Critical Flow Models . . . 15

viii Contents

3.2.8 Analytical Solution . . . 32

3.2.9 Temperature Profile Variation . . . 33

3.2.10 Nodalisation . . . 35 4 Results 37 4.1 MXC-15 . . . 37 4.2 MXC-21 . . . 39 4.3 MXC-20 . . . 40 4.4 MXC-22 . . . 41 4.5 MXC-5 . . . 42 4.6 MXC-14 . . . 44 4.7 MXC-23 . . . 45 4.8 Analytical Solution . . . 46

4.9 Temperature Profile Variation . . . 47

4.10 Nodalisation . . . 50 4.11 Code Uncertainty . . . 53 5 Discussion 56 5.1 MXC-15 . . . 56 5.2 MXC-21 . . . 57 5.3 MXC-20 . . . 57 5.4 MXC-22 . . . 58 5.5 MXC-5 . . . 58 5.6 MXC-14 . . . 59 5.7 MXC-23 . . . 59 5.8 Analytical Solution . . . 60

5.9 Temperature Profile Variation . . . 60

5.10 Nodalisation . . . 61

6 Summary and Conclusions 62

Chapter 1

Introduction

1.1

Nuclear Reactor Technology

The main idea behind any nuclear power plant (NPP) is to transform the heat released from nuclear fission into electricity via a fluid. The fluid acts as a coolant for the nuclear fuel and transports the heat to turbines that generate electricity through rotational motion. Nuclear fission in nuclear power plants is most commonly induced by a thermal neutron being absorbed by a uranium-235 nucleus. The nucleus then fissions into new, lighter atoms and a few neutrons, causing a release of heat. The heat is generated due to the difference in the binding energy of the fissioned nucleus and the combined binding energies of the fission products.

Some of the fission neutrons leak out of the core, some are lost due to absorption in non-fuel material and some are absorbed in the non-fuel. Neutrons that are absorbed in the non-fuel and cause a new fission event are called fission inducing neutrons. The reactor is called critical if for every fission inducing neutron, one new fission inducing neutron is created. This means that the neutron population is kept constant. Furthermore, the neutron population is directly proportional to the reactor power and is the means by which the reactor power is controlled.

In Boiling Water Reactors (BWRs) the chosen coolant is light water which also serves as the moderator, effectively slowing down the neutrons in the core so that they can more readily induce fission. Water is an excellent moderator since the majority of the nuclei in water molecules are hydrogen nuclei, which most often consist of a single proton. Protons are of similar mass to the neutrons and in inelastic collisions, when an incoming particle collides with another particle, the incoming particle loses more kinetic energy if the target particle is of similar mass.

As the coolant is heated, it starts to evaporate and the moderation drops, causing a decrease in the neutron population. Fresh coolant thus has to constantly flow through the core to sustain criticality. This is the purpose of the main recirculation pump. The power of the core can hereby be controlled by changing the speed of the main recirculation pump. Another way to control the neutron population is the insertion and withdrawal of control rods. The control rods consist of neutron absorbing material such as hafnium or boron. When control rods are inserted into the core, the amount of neutron absorbing material in the core increases. This means that the number of free neutrons decreases and the neutron population can be rigorously controlled. In the case of a serious event, a quick reactor shutdown can be initialised by rapid insertion of all control rods - a so called SCRAM.

2 Chapter 1. Introduction

Too much moderation can cause the reactor to become supercritical which means that the neutron population increases with every generation. In the event of supercriticality, the power will increase and in turn more water will evaporate due to the increased heat produced. Since BWRs mostly operate with a negative void coefficient, this means that the neutron population will eventually decrease, returning the reactor to a critical state. The void coefficient describes how much the reactivity increases due to the formation of void in the coolant. A negative void coefficient means that the reactivity decreases with increased void.

Increased heat also increases the stress on the fuel elements. If the temperature is too high for an extended period of time, fuel damage could occur. Sustaining proper core cooling is thus one of the primary concerns in terms of safety of the reactor.

1.2

Nuclear Power Safety

Radioactivity is a natural part of the environment and society benefits from radioactive substances in a great number of areas ranging from medicine to power production. Radiation could, however, be harmful both to the public and the environment. Radiation doses above 0.1-1 Sv can cause nausea, hair loss and damage to the blood, the bone marrow, intestines, and the central nervous system [1]. The risk for various types of cancer can also be increased by increased exposure to radiation [1]. Safety must therefore be ensured in all activities where radioactive material is being handled.

Nuclear power today provides around 11% of the worlds electricity and as much as 40% of the electricity in Sweden [2, 3]. Nuclear fuel is a highly dense source of energy, with a single fuel pellet corresponding to over 800 kg of coal [2]. An NPP can therefore provide a steady supply of electricity for several months between refuellings. However, the safety of the reactor must be ensured at all times. Radioactive release could pose a risk to both the public and the environment. This is highlighted by the consequences in the well known severe accidents at TMI (1979), Chernobyl (1986) and Fukushima (2011).

Even though the regulation of nuclear power is handled by each country separately, there is a potential that the radiological effects can cross borders. For example, the Chernobyl accident was discovered when increased levels of radioactivity was measured by workers at the Swedish power plant Forsmark [4]. International cooperation and common safety standards are therefore vital in nuclear power safety. The International Atomic Energy Agency (IAEA) has, for this purpose, developed a set of safety objectives for nuclear power, the General Safety Objective being [5]:

”To protect individuals, society and the environment by establishing and main-taining in nuclear power plants an effective defence against radiological hazard”.

This means that measures must be taken to ensure minimal probability of an accident having harmful consequences. One of the ways in which the safety objectives are reached is through accident prevention and mitigation. First of all, actions are taken for the prevention of accidents that can potentially cause a loss of control of the reactor core. Secondly, the escalation of accident sequences is to be prevented. Thirdly, the spreading of radioactive substances must be minimised [6].

1.2. Nuclear Power Safety 3

implementation of defence in depth, the goal is to ensure that no single failure can lead to serious accident sequences.

A second means to prevent and mitigate accidents is diversification and redundancy as exemplified by the usage of several backup power generators. Not only should there be several generators but there should be several different types of generators located at separate parts of the facility. If a flooding causes a station blackout and simultaneously floods a generator room, a generator located somewhere else could still start. Forsmark, for example, has a gas turbine located far away from any reactor building that can supply the power plant with electricity in case of a station blackout [7].

Safety assessments are performed to evaluate all safety-significant systems and protection aspects regarding siting, design and operation of nuclear facilities. The safety assessment evaluation process aims to ensure that all requirements set by the regulatory body are met. The safety assessment includes safety analyses of how the reactor behaves during fault conditions. A safety analysis is performed using two complementary methods: Deterministic Safety Analysis (DSA), and Probabilistic Safety Analysis (PSA).

During a DSA, the behaviour of a reactor is studied after an assumed initiating event has occured. The physical processes in the reactor during fault conditions of the normal operating and control systems are considered. Here, the probabilities of different event sequences are disregarded [8]. Any potential releases and consequences are to be shown to not exceed permissible values.

In a PSA the probability of core damage and subsequent damages are evaluated using estimated probabilities for different accident sequences. The PSA is used to identify possible failure scenarios and therefore works as a complement to DSA. Not only are internal hazards recognised, but also external threats such as earthquakes or flooding. The difference between DSA and PSA is that a DSA assumes a postulated scenario and calculates the consequences regardless of the probability of the scenario, whereas a PSA estimates the probabilities and consequences for each step in the accident sequence.

1.2.1

LOCA

4 Chapter 1. Introduction

Figure 1.1. BWR design evolution driven mainly by LOCA. The figure displays, from left to right: an external pump reactor, a jet-pump reactor, and an internal pump reactor. The red cross marks the limiting break location for each reactor design.

In early BWR designs, the main recirculation circuit was located outside the reactor pressure vessel, as exemplified by the leftmost image in Fig. 1.1. Water was taken from the downcomer and pumped by the external recirculation pump into the lower plenum. Such reactors are known as external pump reactors and include the models GE-2, ABB-I and ABB-II. In the case of a LOCA in external pump reactors, the limiting break is a guillotine break located on the main recirculation circuit, after the pump and close to the vessel. In such a case, the core would rapidly become uncovered and the emergency core cooling system (ECCS) therefore had to be fitted with a dedicated spray system to cool the core [10].

The potential of such LOCAs drove the development of BWRs away from external pump reactors. In GE-3, jet-pumps were installed to remove the need for external pumps and allow for the connections of the main recirculation circuit to be placed above the lower plenum [10]. The design is shown in the middle image in Fig. 1.1. This meant that the core could be partially reflooded in the event of a limiting break LOCA. The limiting break location was now moved to the suction line of the jet pump.

Instead, internal pumps were installed in the ABB-III and ABB-IV reactor designs. By integrating the recirculation pumps inside the reactor pressure vessel, the need for large pipes in direct connection to the vessel was eliminated [10]. This design is displayed in the rightmost image in Fig. 1.1. The previously limiting LOCA problem was thereby effectively solved.

Completely rethinking a reactor design is a large and costly venture. However, both General Electric and ABB found it was necessary when faced with the potential problems caused by limiting break LOCA on external pump reactors.

1.2. Nuclear Power Safety 5

will then be completely uncovered and the ECCS has to be fitted with a dedicated spray system above the core in order to adequately cool the exposed fuel elements [10].

0 10 20 30 40 50 60 0 2000 4000 6000 8000 10000 12000 14000 16000

O2: Mass Flow Rate

Time [s]

Break Flow [kg/s]

Total Tank side MCP side

Figure 1.2. Mass flow rate through a break in the main recirculation line for O2.

0 20 40 60 80 100 0 20 40 60 80 100

O2: Reactor mass Inventory

Time [s]

Mass [%]

Figure 1.3. Reactor vessel mass inventory for a break in the main recirculation line for O2.

The O2 LOCA transient shows that a break in the main recirculation line for an external pump reactor drains the coolant very rapidly. The project exemplifies the potential problems with external pump reactor and illustrates the main reason why new types of BWRs were developed.

The governing phenomenon for the flow through a LOCA break is called critical flow, or sometimes choked flow. Critical flow is the main topic of this thesis and will be described in detail in Chapter 2.1.

1.2.2

Computer Codes

In order to ascertain that the safety systems in an NPP function efficiently, several accident analyses are required. For example, the U.S.NRC requires that the licensee shows that the ECCS can handle design basis accidents that have an anticipated frequency of 10−4_per

6 Chapter 1. Introduction

are associated with a nuclear reactor. Performing experiments on such large and complex systems as a nuclear power plant is extremely costly since any danger the experiments could pose must be eliminated. The ability to use computer simulations to ensure plant safety is thus vital for the nuclear power industry.

Two distinct types of code are used in nuclear power safety calculations: deterministic and stochastic codes. In reactor physics, the stochastic codes use the Monte Carlo method to mainly solve particle transport equations [12]. Examples of such codes are MCNP, TRIPOLI and SERPENT. The Monte Carlo method provides a way of solving deterministic problems using a stochastic approach and random numbers. For neutron transport problems in a reactor, a number of independent neutron histories, N , are simulated. The result is then derived from the averaged observation. This method can simplify complex calculations without reduction of the geometry [12]. When calculating the neutron flux in a reactor using the Monte Carlo method, the neutron transport equation does not have to be explicitly formulated. At the same time accurate geometries can be maintained and continuous-energy cross sections can be used. The results are given with statistical errors and for complex problems, the Monte Carlo method can require a high computational cost. The errors generally decay as √1

N so if more accurate results are required, the computational cost

increases.

Deterministic codes instead solve the governing integro-differential equations with bound-ary conditions and empirical correlations using numerical methods [13, 14]. The uncertain-ties in the input data and the accuracy of the models then directly determine the uncertainty in the results. Determinsitic codes, and specifically POLCA-T, are the main focus of this thesis. The computational cost of formulating and solving complex equations can be large. Therefore, simplified models and correlations as well as geometries are often used. Deter-ministic codes thereby give an exact solution to a simplified problem. Instead, stochastic codes can tackle the problem in full complexity but with results that have a statistical uncertainty.

The codes in nuclear power safety calculations can be either conservative or best-estimate. Conservative means that the input parameters are chosen conservatively and that the equations are solved conservatively. This approach is supposed to encompass the reality with an added necessary margin for uncertainties.

By using realistic input parameters and realistic equations instead, realistic results with an associated uncertainty are obtained. The uncertainty in the results is obtained directly from the uncertainty in the input data and the accuracy of the equations and models used. This is what best-estimate codes strive towards. The safety margin is then defined as the margin from the uncertainty limit to the regulatory body acceptance criteria. The margin obtained from conservative analysis is often larger than the uncertainty in the results from best-estimate analysis [13]. Using best-estimate codes for a reactor analysis instead of conservative codes therefore often means that the power of the reactor can be increased without decreasing the safety margin.

1.2. Nuclear Power Safety 7

it arises outside the break location. This means that this model is used whenever the conditions for critical flow is met, if no other model is specified. This model is called the Frozen Flow model and will also be discussed briefly.

Code Assessment

When designing and implementing a computer code, there is a possibility that the developers will make mistakes. If these mistakes result in undetected system faults they could lead to potential safety hazards down the line. It is therefore vital that any potential mistakes are found and corrected before the code is used in real analyses. Currently, the main parts of the code licensing procedure that aim to seek out such mistakes are the verification and validation. Here, the following core questions are to be answered:

• Verification: Are the correct equations solved? • Validation: Are the equations being solved correctly?

During the verification process the computational model is checked to see if it implements the mathematical model correctly. In practice this means that the source code is checked so that it is written correctly and so that there are no possibilities for negative roots, divisions by zero, negative volumes and so forth. During the verification process the source code is also checked so that it matches the code documentation.

The validation process involves comparing simulated results to real life experience and analysing how well the code estimates the measured data. The simulated results can also be compared to the results from previously licensed codes. The latter case is called code comparison. During code comparison, an already validated code is used to simulate the same transient with the same input data. The results are compared to evaluate if the new code performs equivalently or better than the older code. It is also of great importance to evaluate the code uncertainty during the validation. The uncertainties arise from direct code uncertainties, experiment uncertainties and user effect such as nodalisation choice and input errors.

POLCA-T and Other Codes

One deterministic code for nuclear analyses is the Westinghouse developed system code POLCA-T. It is a coupled thermal hydraulic and neutron kinetics code that can be used to simulate 2D thermal hydraulics and 3D kinetics. The hydraulic part solves the mass and energy conservation equations for each phase and control volume. Additionally, POLCA-T uses a mixture momentum conservation equation and a drift flux correlation that describes the transfer of momentum between the two phases. This gives a total of 5 + 1 equations. POLCA-T is therefore called a non-equilibrium code since the phases are not assumed to be in equilibrium. The mass conservation equation for phase i and control volume j can be written

∂ρi,j

∂t = −∇(ρi,jui,j), (1.1) where ρ is the density, t is time, and u is the flow velocity [14]. Similarly, the energy conservation equation for phase i and control volume j is given by

∂

∂t(ei,jρi,j) + ∇(ei,jρi,ju¯i,j) + ∇( ¯qj

00_{) − q}000

j − ∇(τju¯i,j) = 0, (1.2)

where e = ekinetic+ einternal+ epotential is the total specific energy, ¯u is the velocity field

8 Chapter 1. Introduction

stress tensor [14]. Finally, the mixture momentum conservation equation for control volume j can be written ∂ ∂t(ρl,j(1 − αj) ul,j+ ρg,jαjug,j) + ∇ ρl,j(1 − αj) u 2 l,j+ ρg,jαju2g,j
= −∇Pj− ∇τj− ¯ρjg, (1.3)

where index l represents the liquid phase, index g represents the gaseous phase, α is the void fraction, g is the acceleration due to gravity, and ¯ρj = ρl,j(1 − αj) + ρg,jα [14, 18].

POLCA-T also includes correlations for calculations of pressure drop and critical flow rate, as well as equations for boron concentration and non-condensable gases. The critical flow rate correlations will be explained and studied further throughout this thesis. The thermal model solves the material heat conduction equation and calculates the heat transfer from heat structures to the coolant. The code also includes power generation models to calculate the heat generated from fission. The fission power being determined by the two group 3D kinetics model. Finally, various reactor components are modelled in detail in the system model [14].

Several other codes exist. One of the most common is the U.S.NRC developed light water reactor (LWR) non-homogeneous transient analysis code RELAP5 [13]. The code can be used to simulate a large variety of thermal hydraulic transients both within and outside the nuclear power industry. RELAP5 has a limited number of thermal hydraulic channels, however, and is currently being phased out. As a successor the U.S.NRC have developed TRACE based on the codes TRAC-P, TRAC-B, RELAP5 and RAMONA [13].

TRACE does not use a mixture momentum equation, like POLCA-T, but rather handles the momentum conservation of the two phases separately. It is therefore a six equation non-equilibrium thermal hydraulic code. TRACE can be used for transient calculations for both pressurised water reactors (PWRs) and BWRs. Unlike POLCA-T, TRACE has no neutron kinetics besides point kinetics. It is therefore most often coupled with the the 3D neutron kinetic code PARCS. TRACE is a best-estimate transient analysis code and is currently the U.S.NRC’s primary best-estimate tool for predictions of postulated LWR accidents [13]. TRACE also offers a graphical user interface using the tool SNAP. SNAP provides a layout of all of the components in the nuclear power plant model and at the same time gives the user the ability to modify and monitor all of these components.

Another Westinghouse thermal hydraulic code is GOBLIN. GOBLIN is a three equation equilibrium code that can simulate tough transients. The thermal hydraulics model in GOBLIN solves the mass, energy and momentum conservation equations for one phase, meaning GOBLIN is a homogeneous code. In POLCA-T, the liquid and gaseous state are handled by separate equations. Unlike POLCA-T, GOBLIN only uses point kinetics. GOBLIN has previously been validated for critical flow calculations using three cases from Marviken CFT. For more details on the validation of the critical flow models in GOBLIN see the validation documentation [16]. The results from POLCA-T will, for some transients, be compared to GOBLIN results in order to increase the validation basis of POLCA-T.

1.3

Marviken

1.3.1

Marviken Nuclear Power Station

1.3. Marviken 9

also plans to be able to produce weapon grade plutonium in the Marviken reactor [19]. The heavy water was therefore to be imported from Norway because importing it from the U.S. would lead to demands by the U.S. to inspect the reactor so that it was not used for plutonium production [19]. By 1968 the reactor had been fully assembled and non-nuclear testing began.

The reactor design of Marviken was shown to be inherently unstable. In the early cal-culations, the void coefficient as well as the temperature coefficients for fuel and moderator were all shown to be negative. The uncertainty in these calculations was, however, very large. When more accurate calculations were made in the latter part of the 1960s the total reactivity coefficient was shown to become positive after only a small burnup [19]. During the non-nuclear testing, several more issues were found. For example the efficiency of the ECCS was insufficient, the SCRAM system needed replacing and the reactor power was required to be lowered [19]. The decision was made that the changes required to fulfil cur-rent safety standards would be too costly. The project was cancelled in 1970 before fuelling and Sweden instead shifted its eyes to the rapid development of the light water reactor technology.

In January 1970, Sweden ratified the nuclear weapons non-proliferation treaty and thereby abandoned all plans to produce nuclear weapons [20]. The ability to produce weapon grade plutonium was incorporated in the reactor design of Marviken and the ratifi-cation of the nuclear proliferation treaty was therefore another reason to why the Marviken project was abandoned [19]. Left was a fully built reactor that never produced any nuclear power. The plant was later rebuilt into an oil power plant that was still used as a backup power plant in 2007 [19].

Marviken did, however, end up contributing greatly to the nuclear power industry in other ways by than producing power. In addition to the experience gained by the Swedish engineers and nuclear companies during the Marviken NPP project, several multinational full scale experiment projects were conducted between the years 1972-1985. The experi-ments provided data both for simulated LOCAs and for the spreading of aerosols in the primary circuit. One of these experiment projects regarded critical flow.

1.3.2

Marviken Full Scale CFT

Between 1977 and 1979 the Marviken Full Scale Critical Flow Tests (CFTs) were conducted at the Marviken power station. The reactor tank and containment was used to conduct the largest critical flow experiments to date [19]. Other CFTs have been conducted on smaller scales which require extrapolation to be applied to nuclear reactor systems. With nine different exhaust pipe geometries and a wide range of initial conditions a total of 27 CFTs were conducted and thoroughly documented [21]. The Marviken CFT has provided the bulk experimental data used in this thesis.

In the Marviken CFTs, two-phase water mixtures were discharged from the reactor vessel. The vessel had a net volume of 425 m3 _{and a maximum design temperature and}

10 Chapter 1. Introduction

and discharge pipe during blowdown were recorded using a pulse code modulation (PCM) system from 180 s before test initiation until the experiment was terminated. The tests were terminated when steam entered the discharge pipe or when the ball valve in the discharge pipe began to close. Error limits were evaluated for all parameters [21, 23].

Table 1.1. Summary of the nine discharge pipe geometries used in the Marviken Full Scale CFT. The tolerance in the values is ±4 mm [21].

Exhaust pipe diameter, D [mm] Exhaust pipe length, L [mm] L/D 200 590 3.0 300 290 1.0 300 511 1.7 300 895 3.0 300 1116 3.7 500 116 0.3 500 730 1.5 500 1809 3.6 509 1589 3.1

During the tests, the fluid was exhausted into the reactor containment. The containment consisted of a drywell with a volume of 1934 m3_{, a wetwell with a volume of 2144 m}3 _and

the fuel transport hall with a volume of 303 m3 _{[22]. The pressure in the containment was}

relieved by discharging some of the fluid through exhaust tubes located both at ground level and above the steam dome. These venting tubes were connected directly to the atmosphere. The objective of the Marviken CFT project was to obtain the critical mass flow rate ˙

mcritical as a function of stagnation pressure P0, stagnation enthalpy h0, nozzle length L,

nozzle diameter D and air concentration c. The air concentration was determined to have no effect [21]. Efforts were made to keep all other parameters fixed whilst one of them was varied for each test.

In each test, the vessel was filled with de-ionised water and then drained until the required elevation at room temperature was reached. The water was then heated to the specified initial temperature profile by circulating the water from the vessel bottom through an electric heater and back into the steam dome. Some water was also taken directly from the vessel bottom and injected into the steam dome through spray nozzles. The steam dome was periodically vented during the heat-up phase to remove air from the vessel. This produced a temperature profile at the specified pressure with a saturated zone towards the top of the vessel, a transition zone in the center and a subcooled zone in the vessel bottom. The discharge pipe temperatures were somewhat lower than the temperatures in the vessel due to increased heat losses. Hot water was therefore circulated through the discharge pipe from around four hours prior to the blowdown initiation in an effort to minimise this temperature difference. In all cases before test 14, the circulation was ended when the test was initiated. This was observed to produce a non-homogeneous temperature zone in the vessel bottom. After test 13, the circulation was therefore ended 1.5 hours before test initiation.

1.5. Scope 11

gamma densitometer system. These measurement devices were installed in various posi-tions along the vessel and discharge pipe. For more details on the experimental setup see the experiment documentation. [21, 22]. The data was processed and evaluated in order to obtain other parameters such as the mass flow rates and the stagnation conditions.

1.4

Objective

This thesis aims to validate the critical flow models in POLCA-T by comparing the sim-ulation results to the Marviken CFT experimental data and the simsim-ulation results from GOBLIN so that LOCA calculations can be performed using the code in the future. Sec-ondly the thesis provides a stepping stone for moving from conservative calculations to best estimate by providing an estimation of the code uncertainty for the critical flow case.

1.5

Scope

Chapter 2

Theoretical Backgound

2.1

Critical Flow

In thermal hydraulics, when a compressible fluid flows from a high pressure region with a given set of stagnation conditions to a region of lower pressure, the mass flow rate increases with decreasing external pressure. There is, however, an upper limit to how large the mass flow rate can become. The maximum mass flow rate is reached when the flow velocity reaches the local speed of sound at some location in the duct. This location is called the critical cross section. For a straight duct, the critical cross section is usually located at the exit whereas for diverging or converging nozzles the critical cross section can be located both upstream or downstream of the exit [18]. The corresponding flow rate is called critical flow rate and depends only on the stagnation parameters, which means it is independent of the external pressure. For external pressures below a certain critical pressure, the pressure dependence of the mass flow rate thus vanishes.

This phenomenon occurs for both single phase, and two-phase flow systems and has been extensively studied due to its important role in a variety of industrial implementations. One of the largest incentives to the experimental and theoretical study of two-phase critical flow has been its involvement during LOCAs in nuclear power plants [24, 25]. If critical flow is not predicted correctly, any simulated LOCA transients will correspond badly to real scenarios. Without proper implementation of critical flow conditions and models the pressure increase in the containment, the time until dryout, and time when conditions are met so that the low pressure ECCS can start to cool the core can all be predicted wrongly. This could lead to false predictions of the reactor behaviour during transient and hence improper definitions of safety margins. In order to accurately predict the reactor response to LOCA scenarios, sufficient critical flow models are therefore vital.

2.1.1

Single Phase Flow

The easiest way to illustrate single phase critical flow is to consider the analysis of isentropic critical flow in a horizontal tube in one dimension [25]. First of all we write down the mass and momentum conservation equations

2.1. Critical Flow 13

where ˙m is the mass flow rate, ρ is the fluid density, V is the fluid velocity, A is the tube cross sectional area, z is the location along the tube length and P is the pressure.

Critical flow implies that the mass flow rate reaches a maximum value, independently of any pressure changes downstream. Hence,

d ˙m

dP = 0. (2.3)

Differentiating Eq. (2.1) with pressure gives

d ˙m dP = V A dρ dP + ρA dV dP, (2.4)

which together with Eq. (2.3) leads to to

dV dP = − V ρ dρ dP. (2.5)

By applying the critical flow condition from Eq. (2.5) onto Eq. (2.2) the critical mass flow rate ˙mcriticalcan be defined as

˙ mcritical= ρA V dP dρ. (2.6)

The speed of sound of an isentropic fluid, c, is given by

c2= dP dρ



s

. (2.7)

Upon observing Eqs. (2.6) and (2.7) it is clear that the critical flow is exactly the mass flow obtained at sonic speed.

Eq. (2.6) is, however, not always optimal for practical usage since it requires knowledge of the local conditions. By considering a blowdown from a tank through a pipe and relating the fluid velocity at the throat to the upstream enthalpy at stagnation conditions, h0, a

more applicable equation can be obtained [25]. The fluid in such a tank is usually stagnant and at known conditions. In the case of isentropic flow there is no heat addition from friction and hence the stagnation enthalpy is constant along the pipe. This can be expressed as

h0= h +

V2

2 . (2.8)

Inserting Eq. (2.8) into Eq. (2.1) gives

˙

m = ρV A = ρp2(h0− h). (2.9)

For an ideal gas the enthalpy h and temperature T can be related using the specific heat capacity cP by

dh = cPdT. (2.10)

For isentropic expansion of an ideal gas

14 Chapter 2. Theoretical Backgound and T T0 = P P0 (γ−1)/γ , (2.12)

where γ = cP/cV is a ratio between the heat capacities for constant pressure and volume.

By applying Eqs. (2.10), (2.11) and (2.12) to Eq. (2.9) a more applicable formulation of critical flow than Eq. (2.6) can be obtained:

˙ m = ρ0A s 2cPT0  1 − T T0   P P0 2/γ = ρ0A v u u t2cPT0  P P0 2/γ − P P0 (γ+1)/γ! . (2.13)

For values of downstream pressure Pd below some critical pressure, Pcritical, the mass

flow is constant. Above Pcritical, the mass flow varies with the pressure ratio PPd0. To find

the critical pressure Eq. (2.13) can be differentiated with respect to pressure using Eq. (2.3) to find  Pd P0  critical =  ₂ γ + 1 γ−1γ . (2.14)

2.1.2

Two-Phase Flow

In the case of single phase flow, critical flow is reached when the flow velocity reaches the local speed of sound. This is not quite applicable in the more complex two-phase flow. For two-phase fluids, the speed of sound is not the same for the liquid and the gaseous phase. Neither are the velocities of the two phases necessarily the same depending on the flow regime. Instead, consider flow from stagnant conditions to some point in a pipe [16, 24]. Using the indices g for the vapour phase and l for the liquid phase, the energy balance can then be written ˙ mh0= ˙mg hg+ V2 g 2 ! + ˙ml  hl+ V2 l 2  . (2.15)

By introducing steam quality x = m˙g

˙

m and slip factor S = Vg

Vl Eq. (2.15) can be rewritten

as h0= h + 1 2  ˙m A  x ρg +1 − x ρl S 2 x + 1 − x S2  . (2.16)

Eq. (2.16) is further simplified by defining the specific volume v as

2.1. Critical Flow 15

where G = m_A˙ is the mass flux. For frictionless flow the entropy s remains constant and the impulse equation is hence written

s = s0. (2.19)

This means that if all thermodynamic quantities are viewed as functions of pressure and enthalpy, they only depend on pressure. Since critical flow is equivalent to maximum mass flux, we find

dG

dP = 0, (2.20)

which together with Eq. (2.18) gives

0 = 2GdG dP = −2dh dPv 2_{− 2 (h} 0− h) 2v_dPdv
v4 = −2dh dP − G 2 _2vdv dP
v2 . (2.21)

The flow velocity can be written V = Gv so Eq. (2.21) can be rewritten as

Vcritical= s v_dPdh −dv dP = c, (2.22)

2.1.3

Critical Flow Models

There are numerous models of varying complexity developed for critical two-phase flow. The discussion of these models has been limited to the ones that are examined in the thesis.

Homogeneous Equilibrium Model (HEM)

16 Chapter 2. Theoretical Backgound o ..--< * o -""" Rapport 97-12-05

Kritisk strömning enligt HEM

PA 97-165 Figur 3.3-1 0~~--t-~--~~~==~=i==~=f~~~~==~~ 0.0 0.4 0.8 1.2 1.6 Enthalpy (J/kg) 2.0 *10 6 2.4 2.8

Figure 2.1. The critical mass flux as a function of stagnation enthalpy and pressure as calculated by HEM from Ref. [16].

Fig. 2.1 shows how the critical mass flux calculated by HEM behaves for varying stag-nation pressure and enthalpy. The saturation lines are represented by dashed lines perpen-dicular to the mass flux lines. It is evident that the flow of subcooled fluids is much greater than for saturated fluids. An increased stagnation pressure also means an increased flow. For all cases analysed in this thesis, the pressure is between 40 and 50 bar and the lowest enthalpy is around 1 MJ/kg meaning a maximum mass flux of almost 105 kg

m2sis expected in

the HEM calculations. In the saturated region, the mass flux is very small, and differences in pressure and enthalpy are less significant.

HEM usually gives good agreement with data in nuclear applications [16]. This is especially true for pipes of longer length were the time is sufficient to achieve equilibrium. The error can be larger for shorter pipes or if the flow regime allows large flow velocity difference between phases. An example of such a flow regime is annular flow. HEM then generally gives an underestimation of the flow [16, 28].

2.1. Critical Flow 17

data during these years. HEM has therefore been chosen by Westinghouse to serve as the best-estimate model for critical flow in POLCA-T.

Moody

Moody’s model for critical flow was proposed by F. J. Moody in 1965 and will henceforth be referred to simply as Moody [29]. Unlike HEM, Moody does not assume total phase equilibrium. Instead Moody views the steam and water as saturated and allows the slip factor to be a free parameter. The slip factor is then chosen to maximise the flow according to dG dS    _P = 0, (2.23)

which leads to the condition

Smax= Vg Vl    _max =r ρ3 l ρg. (2.24)

Given saturation conditions for the enthalpy and pressure, Moody first determines the stag-nation entropy and chooses the pressure so that the flow velocity is equal to the local speed of sound.

18 Chapter 2. Theoretical Backgound C') o ,..., * o o N o 00 ,..., o \0 ,..., o 'T ,..., o N ,..., o \0 o N Rapport 97-12-05

Kritisk strömning enligt Moody

PA 97-165 Figur 3.2-1

o~~~~~~~~~~~~

0.0 0.4 0.8 1.2 1.6 2.4 2.8 Enthalpy (J/kg)

Figure 2.2. The critical mass flux as a function of enthalpy and pressure as calculated by Moody from Ref. [16].

Fig. 2.2 shows how the critical mass flux calculated by Moody behaves for varying stag-nation pressure and enthalpy. The saturation lines bound the region of validity for Moody. In comparing Fig. 2.1 and Fig. 2.2, it is obvious that Moody calculates a higher flow in the saturated region than HEM does. This is a direct effect of allowing the slip to be a free, maximised parameter. For subcooled fluids the flow is single-phase, meaning the slip is undefined and therefore cannot be maximised. Moody is therefore invalid for single-phase flow and a separate flow model should be used. Fig. 2.2 therefore shows the critical mass flux predicted by HEM for subcooled flow.

Frozen Flow

2.2. Uncertainty Analysis 19

2.2

Uncertainty Analysis

In order to move from conservative to best estimate methodology, the accuracy of the model results must be known. In this thesis the uncertainties have been estimated using two separate methods. The methods are to be seen as complementary since they both have strengths and weaknesses, as will be discussed once the methods have been described.

2.2.1

Normality Method

The first method is derived from the assumption that the errors are distributed normally. This method will henceforth be called the Normality Method. Results from an arbitrary case are shown in Fig. 2.3.

0 10 20 30 40 50 60 0 2000 4000 6000 8000 10000 12000 Time [s] Break flow [kg/s] Experimental Result, r(t) Model Result, m(t) Experimental Uncertainty, r(t) ±σ_r(t)

Figure 2.3. The experimentally measured and simulated mass flow rate from an arbitrary case. The standard deviation interval of the experimental results is shown as dashed red lines.

Fig. 2.3 shows mass flow rate, but the following argument is valid for any parameter. The time dependent experimental results of some property r(t), with the corresponding time dependent standard deviation, σr(t) and the results from the simulation m(t) are

known. For the Marviken CFT, the standard deviations for every parameter in the exper-imental results have been evaluated [23]. The problem is to find the standard deviation, σm, associated with the simulated results. The models perform differently in different flow

regimes and for different phases and each blowdown is characterised by different zones of various degrees of subcooling for different times during the discharge. This means that the standard deviations change with each timestep i. Consequently, σm,i has to be found for

each timestep separately and an average then has to be calculated according to

20 Chapter 2. Theoretical Backgound

We know that, for a certain timestep i, the correct result lies within 1σ of the ex-perimental result with a probability of P1σ = 68.269%. The interval can be written as

(ri− σr,i, ri+ σr,i). The probability for normally distributed errors is defined as

P = Z ri+σr,i ri−σr,i 1 σr,i √ 2πe −(x−ri)2 2σ2_r,i dx. (2.26)

We want to know the standard deviation σm,ifor the model result miso that the correct

result (note that the correct result is not the same as the experimental result since the experimental result is associated with errors) falls within the interval (mi− σm,i, mi+ σm,i)

with a probability of P1σ. This effectively means that we want to solve the following equation

for σm,i: P1σ= Z mi+σm,i mi−σm,i 1 σr,i √ 2πe −(x−ri)2 2σ2_{r,i dx. (2.27)}

I have chosen to solve Eq. (2.27) iteratively. In the iteration process, a starting guess is made for σm,i and the probability is then calculated. The probability is checked to see if

P = P1σ. Otherwise, if the resulting probability is P > P1σ, a smaller interval is chosen.

If, on the other hand, the probability is P < P1σ, a larger interval is chosen. The iterative

solution of Eq. (2.27) is performed for each timestep i separately.

The above described iteration process will yield a set of results σm,i. This is a different

standard deviation for each time step in the simulation. As previously stated, the mean of the standard deviations is taken according to Eq. (2.25) before being presented in this thesis as a final result. The relative error δmis also presented and is calculated as a mean

of the relative error of each timestep δm,i according to

δm= 1 N N X i=1 δm,i, (2.28) where δm,i= σm,i ri . (2.29)

The Normality Method assumes that the errors are distributed normally in each timestep, an assumption that may not necessarily be correct. There are several normality checks available but neither can be implemented here since there is only one measurement point for each timestep and each timestep has to be treated separately due to the time dependence of σm and σr. Therefore, a separate uncertainty estimation method is also used.

2.2.2

Ratio Method

2.2. Uncertainty Analysis 21

whether the difference between the simulated result and the real result is smaller or larger than the difference between the simulated result and the experimental result.

This means that if the results from the Normality Method and the Ratio Method are similar, it can be assumed that the results from the Normality Method are correct. Since the Ratio Method does not include the experimental uncertainty in the estimation, the uncertainty in the estimation of the model error using the Ratio Method is directly deter-mined by the experimental uncertainty. Therefore, if the absolute difference between the estimated relative uncertainty by the Ratio Method and the Normality Method is larger than the relative uncertainty of the experiment (|δRM− δNM| > δr) the Normality Method

estimation cannot be trusted. The error can then not be assumed to be distributed nor-mally. If, instead, the difference is smaller than the relative uncertainty of the experiment (|δRM− δNM| < δr) the Normality Method estimation is within the uncertainty limits of

the Ratio Method, meaning that the more accurate Normality Method should be trusted. The mean relative uncertainty in the experiment measurement of the flow is around 6%, depending on case, whilst the relative experimental uncertainty for the pressure is around 2%, depending on case [21].

The main idea behind the Ratio Method is to calculate the error by summing the absolute distance between the experimental result and the model result for every timestep. The relative uncertainty is then given as a ratio according to

Chapter 3

Methodology

3.1

Geometrical Model

The detailed description available from the Marviken CFT meant that a highly accurate geometrical model could be built for the simulations [22]. The geometrical model is shown in Fig. 3.1. The colour scale in the figure represents void, where red means pure steam and green means pure liquid. A total of 61 cylindrical nodes were used. These nodes were divided into three parts, namely VESSEL, NOZZLE and TEST as seen in Fig. 3.1. The diameter of every node was adapted so that the total volume of the reactor pressure vessel was equal to the real version. The geometrical model includes the internal structures and the plot in Fig. 3.1 displays the free volume. This is the reason why the diameter of the nodes in the middle if the vessel in Fig. 3.1 is not constant.

The height of each node is 0.48 m for the VESSEL, 0.90 m for the NOZZLE and between 0.05 and 0.60 m for the TEST, depending on which discharge pipe was studied. The ther-modynamic quantities were defined nodewise for each analysis case according to the initial data described in each respective experiment report [30, 31, 32, 33, 34, 35, 36].

3.2. Experiment Matrix 23 VESSEL_C36 VESSEL_C23 VESSEL_C49 VESSEL_C03 VESSEL_C19 NOZZLE_C04 TEST_C02 VESSEL_C32 VESSEL_C29 VESSEL_C27 TEST_C03 VESSEL_C08 VESSEL_C44 VESSEL_C34 VESSEL_C10 NOZZLE_C01 VESSEL_C47 VESSEL_C35 VESSEL_C42 VESSEL_C39 VESSEL_C07 NOZZLE_C03 VESSEL_C41 VESSEL_C45 VESSEL_C48 VESSEL_C21 VESSEL_C05 VESSEL_C02 VESSEL_C12 VESSEL_C14 VESSEL_C17 VESSEL_C20 NOZZLE_C02 NOZZLE_C07 VESSEL_C24 VESSEL_C18 VESSEL_C04 VESSEL_C06 NOZZLE_C06 VESSEL_C01 VESSEL_C26 VESSEL_C37 VESSEL_C30 TEST_C01 VESSEL_C50 VESSEL_C13 VESSEL_C16 VESSEL_C15 VESSEL_C46 VESSEL_C22 VESSEL_C25 NOZZLE_C05 VESSEL_C33 VESSEL_C28 VESSEL_C11 VESSEL_C43 VESSEL_C51 VESSEL_C09 VESSEL_C31 VESSEL_C40 VESSEL_C38 Time 0.000e+00 s 1.000e+00 0.000e+00 5.000e−01

Figure 3.1. The geometrical vessel model for MXC-15 built according to experiment docu-mentation [21, 30]. The models for the other cases vary in terms of water level, nodalisation, and the geometry of the TEST-pipe.

3.2

Experiment Matrix

24 Chapter 3. Methodology

exhaust pipe geometries are similar to pipe sizes in real nuclear power plants. Additionally, the different cases should ideally only change one parameter per case. This makes it easier to identify which parameters impact the accuracy of the critical flow models the most. A summary of each case in the Marviken CFT and some of their key features are presented in Tab. 3.1.

Table 3.1. Summary of the Marviken CFT matrix [21].

Case number Exhaust pipe diameter D [mm] Exhaust pipe length L [mm] Nominal sub-cooling, Tsc [◦C] Experimental peak mass flow rate [kg/s] MXC-1 300 895 20 3750 MXC-2 300 895 38 440 MXC-3 500 1589 19 11250 MXC-4 500 1589 37 12500 MXC-5 500 1589 33 10950 MXC-6 300 290 31 4560 MXC-7 300 290 18 4500 MXC-8 500 1589 35 12300 MXC-9 500 1589 2 8930 MXC-10 500 1589 3 8800 MXC-11 500 1589 35 12300 MXC-12 300 895 33 4380 MXC-13 200 590 31 2390 MXC-14 200 590 3 2180 MXC-15 500 1809 31 11900 MXC-16 500 1809 33 12400 MXC-17 300 1116 31 4470 MXC-18 300 1116 32 4650 MXC-19 300 1116 4 5180 MXC-20 500 730 7 10050 MXC-21 500 730 33 12180 MXC-22 500 730 52 13700 MXC-23 500 166 3 11700 MXC-24 500 166 33 13000 MXC-25 300 511 6 5390 MXC-26 300 511 34 4530 MXC-27 500 730 33 12900

3.2. Experiment Matrix 25

Table 3.2. Summary of the cases chosen to analyse in this thesis [21]. *Nominal superheating [◦_C]. Case Number Exhaust pipe diameter D [mm] Exhaust pipe length L [mm] L/D Steam Dome pres-sure [MPa] T = Tsat−Tsc [◦C] Initial vessel water level [m] Experimental peak mass flow rate [kg/s] MXC-15 500 1809 3.6 5.04 233=264-31 19.93 11900 MXC-21 500 730 1.5 4.94 230=263-33 19.95 12180 MXC-20 500 730 1.5 4.99 257=264-7 16.65 10050 MXC-22 500 730 1.5 4.93 211=263-52 19.64 13700 MXC-5 500 1589 3.1 4.06 218=251-33 17.44 10950 MXC-14 200 590 3.1 4.97 261=264-3 18.10 2180 MXC-23 500 166 0.3 4.96 260=263-3 19.85 11700 ANA 500 1809 3.6 5.05 600=265+335* - 1000

The cases here have been chosen so that, to the extent that it is possible, any one case differs in only one parameter to another case. For example MXC-15 and MXC-21 only differ in exhaust pipe length, and MXC-21 and MXC-22 only differ in subcooling temperature.

The initial state of the vessel was recorded in the Marviken CFT using 25 temperature measurements at varying elevations in the vessel. This means that a very accurate initial temperature profile could be created for each case [21]. Input temperature profiles were created that followed the experimental data as closely as possible. However there were no temperature measurements performed in the final part of the vessel geometry called TEST. For simplicity, the temperature in the three TEST nodes were assumed to have a constant temperature equal to the temperature of the bottommost node of the NOZZLE, NOZZLE C07. This may have been a systematic error in the test procedure and is discussed further in Chapter 3.2.9.

The pressure, on the other hand, was only measured at the top and bottom of the vessel so some approximations had to be made. It was assumed that the pressure was constant in the steam dome and then increased linearly to the bottom of VESSEL. See for example Fig. 3.3. In NOZZLE and TEST, the pressure was assumed constant at the same value as the final node in VESSEL.

3.2.1

MXC-15

26 Chapter 3. Methodology 180 190 200 210 220 230 240 250 260 270 −5 0 5 10 15 20 25 Temperature [ °C] Elevation [m] Temperature Profile MXC−15 Measurement Vessel Nozzle Test

Figure 3.2. Initial temperature profile in the setup for case MXC-15.

5.04 5.06 5.08 5.1 5.12 5.14 5.16 5.18 5.2 −5 0 5 10 15 20 25 Pressure [MPa] Elevation [m] Pressure Profile MXC−15 Vessel Nozzle Test Measurement

Figure 3.3. Initial pressure profile in the setup for case MXC-15.

3.2.2

MXC-21

3.2. Experiment Matrix 27 180 190 200 210 220 230 240 250 260 270 −5 0 5 10 15 20 25 Temperature [ °C] Elevation [m] Temperature Profile MXC−21 Measurement Vessel Nozzle Test

Figure 3.4. Initial temperature profile in the setup for case MXC-21.

4.95 5 5.05 5.1 −5 0 5 10 15 20 25 Pressure [MPa] Elevation [m] Pressure Profile MXC−21 Vessel Nozzle Test Measurement

Figure 3.5. Initial pressure profile in the setup for case MXC-21.

3.2.3

MXC-20

28 Chapter 3. Methodology 180 190 200 210 220 230 240 250 260 270 −5 0 5 10 15 20 25 Temperature [ °C] Elevation [m] Temperature Profile MXC−20 Measurement Vessel Nozzle Test

Figure 3.6. Initial temperature profile in the setup for case MXC-20.

5 5.02 5.04 5.06 5.08 5.1 5.12 −5 0 5 10 15 20 25 Pressure [MPa] Elevation [m] Pressure Profile MXC−20 Vessel Nozzle Test Measurement

Figure 3.7. Initial pressure profile in the setup for case MXC-20.

3.2.4

MXC-22

3.2. Experiment Matrix 29 170 180 190 200 210 220 230 240 250 260 270 −5 0 5 10 15 20 25 Temperature [ °C] Elevation [m] Temperature Profile MXC−22 Measurement Vessel Nozzle Test

Figure 3.8. Initial temperature profile in the setup for case MXC-22.

4.9 4.95 5 5.05 5.1 −5 0 5 10 15 20 25 Pressure [MPa] Elevation [m] Pressure Profile MXC−22 Vessel Nozzle Test Measurement

Figure 3.9. Initial pressure profile in the setup for case MXC-22.

3.2.5

MXC-5

30 Chapter 3. Methodology 210 220 230 240 250 −5 0 5 10 15 20 25 Temperature [ °C] Elevation [m] Temperature Profile MXC−5 Measurement Vessel Nozzle Test

Figure 3.10. Initial temperature profile in the setup for case MXC-5.

4.04 4.06 4.08 4.1 4.12 4.14 4.16 4.18 4.2 4.22 −5 0 5 10 15 20 25 Pressure [MPa] Elevation [m] Pressure Profile MXC−5 Vessel Nozzle Test Measurement

Figure 3.11. Initial pressure profile in the setup for case MXC-5.

3.2.6

MXC-14

3.2. Experiment Matrix 31 170 180 190 200 210 220 230 240 250 260 −5 0 5 10 15 20 25 Temperature [ °C] Elevation [m] Temperature Profile MXC−14 Measurement Vessel Nozzle Test

Figure 3.12. Initial temperature profile in the setup for case MXC-14.

4.96 4.98 5 5.02 5.04 5.06 5.08 5.1 5.12 −5 0 5 10 15 20 25 Pressure [MPa] Elevation [m] Pressure Profile MXC−14 Vessel Nozzle Test Measurement

Figure 3.13. Initial pressure profile in the setup for case MXC-14.

3.2.7

MXC-23

32 Chapter 3. Methodology 50 100 150 200 250 −5 0 5 10 15 20 25 Temperature [ °C] Elevation [m] Temperature Profile MXC−23 Measurement Vessel Nozzle Test

Figure 3.14. Initial temperature profile in the setup for case MXC-23.

4.96 4.98 5 5.02 5.04 5.06 5.08 5.1 5.12 −5 0 5 10 15 20 25 Pressure [MPa] Elevation [m] Pressure Profile MXC−23 Vessel Nozzle Test Measurement

Figure 3.15. Initial pressure profile in the setup for case MXC-23.

3.2.8

Analytical Solution

Although not a case from the Marviken CFT, the same geometrical model as in MXC-15 was also used for a single phase case. The tank was filled with pure superheated steam at a pressure of 5.05 MPa and a temperature of 600◦C, meaning 335◦C above the saturation temperature. This was analysed for the critical flow models HEM and Frozen flow, but not for Moody. Moody assumes a two-phase fluid and tries to maximise the slip. It is therefore not applicable to the single phase case.

3.2. Experiment Matrix 33

compared to an experiment. The governing differential equation for the pressure evolution P (t) during blowdown of a perfect gas is [37]:

 P (t) P0 1−3γ2γ _d(P/P 0)) dt = γA ρ0V p γP0ρ0  ₂ γ + 1 _2(γ−1)γ+1 , (3.1)

with initial condition

P (0) P0

= 1. (3.2)

In Eq. (3.1) P0is the stagnation pressure, γ = cP/cV is the specific heat capacity ratio

(γ ≈ 1.3 for steam [37]), A is the exhaust pipe exit area, ρ0 is the stagnation density, and

V is the vessel volume. A solution to Eq. (3.1) for isothermal blowdown is given by [37]:

P (t) = P0exp −  ₂ γ + 1 2(γ−1)γ+1 s γP0 ρ0 At V ! . (3.3)

The pressure solution can then be inserted into the single phase mass flow equation as presented in Eq. (2.13).

3.2.9

Temperature Profile Variation

The purpose of this study was to perform a sensitivity analysis on the initial temperature profile in the vessel. The experiment reports describe in detail the initial temperature in all parts of the setup except for the last 0.1-1.5 m, meaning the majority of the TEST sec-tion [30, 31, 32, 33, 34, 35, 36]. The exhaust pipe had no heating elements and less insulasec-tion than the rest of the setup [21]. Combined with the fact that, at these temperatures, the density of water increases with decreasing temperature, cold water should agglomerate in the exhaust pipe. This phenomenon is observed in some of the temperature profile figures above, see especially Figs. 3.14 and 3.12. In the experiments, water was circulated through the discharge pipe prior to each test initiation in order to equalise the temperature in the bottom of the vessel [21]. The circulation was terminated at the start of the test for all cases before MXC-14 and 1.5 hours before the test initiations for the remaining cases.

In the sensitivity analysis the temperature profile was varied in the area where no ex-perimental data was available. In all other cases, the temperature was assumed constant in the TEST part. One of the cases was chosen where the temperature in the TEST part was varied in order to study the effects of the previous simplification.

The case chosen to perform this parameter study on was MXC-5. This case had a particularly long time between the establishing of initial conditions and the start of the CFT. For MXC-5 this time difference was 53 minutes whereas for the other cases the time difference was around 20-30 minutes [30, 31, 32, 33, 34, 35, 36]. This extended time period means that the effects described above should be even more predominant for this case than any other. This decision was made before knowing that the recirculation in the discharge pipe was stopped at different times for early and late cases. In hindsight, a case after MXC-13 could additionally have been analysed. However, the early cases should have a more non-homogeneous temperature in the bottom part of the vessel and discharge pipe which could lead to other temperature effects.

34 Chapter 3. Methodology

in Fig. 3.10. In the second profile, shown in Fig. 3.16, the temperatures of the TEST nodes come from a linear extrapolation of the three bottommost experimental measurements.

50 100 150 200 250 −5 0 5 10 15 20 25 Temperature [°C] Elevation [m]

Extrapolated Temperature Profile Measurement

Vessel Nozzle Test

Figure 3.16. Initial temperature profile for case MXC-5 using linear extrapolation.

In the three temperature profiles shown in Figs. 3.17, 3.18 and 3.19, the hypothesis that there is a lump of cold water in the end of the nozzle was tested. The slope between the temperature of the three nodes is assumed linear in all three cases, but the coldest node is chosen to have a different temperature for each case. The coldest temperatures were chosen to be 150, 105 and 50◦_C. 50 100 150 200 250 −5 0 5 10 15 20 25 Temperature [°C] Elevation [m]

Temperature Profile with Subcooled Lump at 150 °C Measurement

Vessel Nozzle Test

3.2. Experiment Matrix 35 50 100 150 200 250 −5 0 5 10 15 20 25 Temperature [°C] Elevation [m]

Temperature Profile with Subcooled Lump at 105 °C Measurement

Vessel Nozzle Test

Figure 3.18. Initial temperature profile for case MXC-5 using a cold lump of water in TEST with the coldest temperature at 105◦C.

50 100 150 200 250 −5 0 5 10 15 20 25 Temperature [°C] Elevation [m]

Temperature Profile with Subcooled Lump at 50 °C Measurement

Vessel Nozzle Test

Figure 3.19. Initial temperature profile for case MXC-5 using a cold lump of water in TEST with the coldest temperature at 50◦_C.

3.2.10

Nodalisation

36 Chapter 3. Methodology VESSEL_C23 VESSEL_C20 NOZZLE_C02 VESSEL_C03 VESSEL_C19 VESSEL_C24 VESSEL_C18 NOZZLE_C04 VESSEL_C06 VESSEL_C04 TEST_C02 VESSEL_C01 TEST_C01 VESSEL_C08 VESSEL_C10 NOZZLE_C01 VESSEL_C13 VESSEL_C07 VESSEL_C16 NOZZLE_C03 VESSEL_C15 VESSEL_C25 VESSEL_C22 VESSEL_C11 VESSEL_C21 VESSEL_C05 VESSEL_C02 VESSEL_C12 VESSEL_C09 VESSEL_C14 VESSEL_C17 Time 0.000e+00 s 1.000e+00 0.000e+00 5.000e−01

Figure 3.20. The geometrical vessel model for MXC-15 built according to experiment docu-mentation using a total of 31 cylindrical nodes [21, 30].

VESSEL_C044 VESSEL_C087 TEST_C05 NOZZLE_C10 VESSEL_C067 NOZZLE_C04 VESSEL_C063 TEST_C02 VESSEL_C091 TEST_C03 VESSEL_C098 VESSEL_C065 TEST_C04 VESSEL_C096 VESSEL_C090 VESSEL_C017 NOZZLE_C11 VESSEL_C050 VESSEL_C018 TEST_C06 VESSEL_C041 VESSEL_C073 VESSEL_C088 VESSEL_C051 VESSEL_C040 VESSEL_C043 VESSEL_C009 NOZZLE_C07 VESSEL_C007 VESSEL_C014 VESSEL_C099 VESSEL_C013 VESSEL_C022 NOZZLE_C06 VESSEL_C039 TEST_C01 VESSEL_C059 VESSEL_C066 VESSEL_C057 VESSEL_C083 VESSEL_C037 VESSEL_C010 VESSEL_C097 VESSEL_C068 NOZZLE_C05 VESSEL_C101 VESSEL_C062 VESSEL_C020 VESSEL_C092 VESSEL_C089 VESSEL_C012 VESSEL_C005 VESSEL_C072 VESSEL_C084 VESSEL_C035 VESSEL_C003 VESSEL_C011 VESSEL_C100 VESSEL_C036 VESSEL_C102 VESSEL_C026 VESSEL_C038 NOZZLE_C12 VESSEL_C024 VESSEL_C093 VESSEL_C002 VESSEL_C086 VESSEL_C032 VESSEL_C054 VESSEL_C094 VESSEL_C079 VESSEL_C008 VESSEL_C015 VESSEL_C027 VESSEL_C028 VESSEL_C006 VESSEL_C030 VESSEL_C046 VESSEL_C085 VESSEL_C077 VESSEL_C058 NOZZLE_C01 VESSEL_C034 NOZZLE_C08 VESSEL_C029 NOZZLE_C03 VESSEL_C070 VESSEL_C082 VESSEL_C056 VESSEL_C019 VESSEL_C031 VESSEL_C060 VESSEL_C095 NOZZLE_C13 NOZZLE_C09 VESSEL_C081 VESSEL_C071 VESSEL_C016 VESSEL_C033 VESSEL_C053 VESSEL_C004 NOZZLE_C02 VESSEL_C061 VESSEL_C001 VESSEL_C025 VESSEL_C078 VESSEL_C064 VESSEL_C074 VESSEL_C049 VESSEL_C080 VESSEL_C075 VESSEL_C042 VESSEL_C048 NOZZLE_C14 VESSEL_C076 VESSEL_C047 VESSEL_C045 VESSEL_C055 VESSEL_C023 VESSEL_C021 VESSEL_C069 VESSEL_C052 1.000e+00 0.000e+00 5.000e−01 Time 0.000e+00 s

Chapter 4

Results

The parameters chosen to be analysed in the results are limited to the mass flow rate through the break and the pressure in the bottom of the vessel. The pressure in other parts of the vessel follows the same evolution as in the bottom and other parameters such as the enthalpy, and water level depend directly on the pressure and mass flow rate. Graphs displaying the experimental results and the simulated results from POLCA-T using both HEM and Moody are presented for break flow and vessel bottom pressure in each case. When simulated results from GOBLIN are available, these are also presented. Finally, the estimated model uncertainty is also presented for each case.

4.1

MXC-15

Fig. 4.1 shows the mass flow rate through the break during case MXC-15 for both POLCA-T and GOBLIN. POLCA-The experimental data shows a high flow of the subcooled liquid initially present in the vessel bottom until around 20 s. After this, a flow region with two-phase flow is observed until the experiment termination. GOBLIN is seen to generally produce a higher mass flow than POLCA-T. Moody overestimates the flow in the two-phase region and underestimates the flow in the subcooled region in both codes. HEM seems to predict the flow very well for both codes. GOBLIN using Moody overestimates the flow so much that the mass inventory of the vessel runs out 5 s before the end of the transient.

38 Chapter 4. Results 0 5 10 15 20 25 30 35 40 45 50 55 0 2000 4000 6000 8000 10000 12000 Marviken MXC−15 Time [s] Break flow [kg/s] Experimental data POLCA−T HEM POLCA−T Moody GOBLIN 4.0.0 HEM GOBLIN 4.0.0 Moody

Figure 4.1. Mass flow rate through the break for case MXC-15.

The vessel bottom pressure of case MXC-15 for both codes is shown in Fig. 4.2. HEM describes the pressure drop very accurately in both codes. Moody overestimates the pressure in POLCA-T throughout the transient whereas in GOBLIN it follows the experimental data closely until a few seconds after the two-phase flow region begins.

0 5 10 15 20 25 30 35 40 45 50 55 0 1 2 3 4 5 6 Marviken MXC−15 Time [s]

Pressure at reactor bottom [MPa]

Experimental data POLCA−T HEM POLCA−T Moody GOBLIN 4.0.0 HEM GOBLIN 4.0.0 Moody

4.2. MXC-21 39

4.2

MXC-21

Figs. 4.3 and 4.4 show the break flow and bottom pressure for case MXC-21, respectively. Both POLCA-T and GOBLIN simulations are shown. A single phase flow region is observed in the experimental data until after around 30 s. In the single phase region, all simulated transients underestimate the flow with the exception of GOBLIN HEM for the first 5 s. Moody for the two codes differ rather substantially. GOBLIN using Moody predicts a flow of around 1000 kg/s greater than POLCA-T Moody. In the two-phase flow region, HEM slightly underestimates the flow whereas Moody overestimates the flow.

0 10 20 30 40 50 60 0 2000 4000 6000 8000 10000 12000 Marviken MXC−21 Time [s] Break flow [kg/s] Experimental data POLCA−T HEM POLCA−T Moody GOBLIN 4.0.0 HEM GOBLIN 4.0.0 Moody

Figure 4.3. Mass flow rate through the break for case MXC-21.