ROYAL INSTITUTE OF TECHNOLOGY
SENSITIVITY AND UNCERTAINTY ANALYSIS OF BWR STABILITY
IVAN GAJEV
LICENTIATE THESIS
STOCKHOLM, SWEDEN, 2010
i Abstract
Best Estimate codes are used for licensing, but with conservative assumptions. It is claimed that the uncertainties are covered by the conservatism of the calculation.
As Nuclear Power Plants are applying for power up-rates and life extension, evaluation of the uncertainties could help improve the performance, while staying below the limit of the safety margins.
Given the problem of unstable behavior of Boiling Water Reactors (BWRs), which is known to occur during operation at certain power and flow conditions, it could cause SCRAM and decrease the economic performance of the plant. Performing an uncertainty analysis for BWR stability would give better understating of the phenomenon and it would help to verify and validate (V&V) the codes used to predict the NPP behavior.
This thesis reports an uncertainty study of the impact of Thermal-Hydraulic, Neutronic, and Numerical parameters on the prediction of the stability of the BWR within the framework of OECD Ringhals-1 stability benchmark. The time domain code TRACE/PARCS was used in the analysis. This thesis is divided in two parts: Sensitivity study on Numerical Discretization Parameters (Nodalization, Time Step, etc.) and Uncertainty part.
A Sensitivity study was done for the Numerical Parameters (Nodalization and Time step). This was done by refining all possible components until obtaining Space-Time Converged Solution, i.e. further refinement doesn’t change the solution. When the space-time converged solution was compared to the initial discretization, a much better solution has been obtained for both the stability measures (Decay Ratio and Frequency) with the space-time converged model.
Further on, important Neutronic and Thermal-Hydraulic Parameters were identified and the uncertainty calculation was performed using the Propagation of Input Errors (PIE) methodology. This methodology, also known as the GRS method, has been used because it has been tested and extensively verified by the industry, and because it allows identifying the most influential parameters using the Spearman Rank Correlation.
Key Words: BWR Stability, Sensitivity, Uncertainty.
ii
iii Table of Contents
Abstract ... i
List of Figures ... v
List of Abbreviations ... vii
Acknowledgments ... ix
1 INTRODUCTION ... 1
1.1 Boiling Water Reactor Stability ... 1
1.2 System Codes ... 3
1.3 State-of-the-Art in BWR Stability ... 5
1.4 Overview of Uncertainty Methodologies ... 7
1.4.1 Propagation of Input Errors Methodology ... 10
1.4.2 Sources of Uncertainties ... 13
1.5 Aim of the Thesis ... 18
2 BWR STABILITY TEST CASE ... 19
2.1 TRACE/PARCS Model Description ... 19
2.2 Model Validation... 23
3 UNCERTAINTY RESULTS ... 33
3.1 Selection of Input Parameters ... 33
3.2 PIE Calculation Evaluation ... 36
3.3 Spearman Correlation Results ... 38
4 CONCLUSIONS ... 41
References ... 43
Appendix A. Spearman Correlations for Cycle 14. ... 47
iv
v
List of Figures
Figure 1 Oskarshamn-2 1999 event ... 2
Figure 2 Calculation of stability parameters – Decay Ratio and Frequency ... 3
Figure 3 Classification of models adopted in the stability analysis of BWRs. ... 4
Figure 4 Scheme of the Safety Margin ... 8
Figure 5 Propagation of input errors, Propagation of output errors .... 10
Figure 6 Illustration of the PIE method and results interpretation ... 11
Figure 7 Diagram of the sources of uncertainties ... 17
Figure 8 Initial Nodalization Scheme ... 20
Figure 9 Results from the base cases. Cycles 14, 15, 16, 17... 22
Figure 10 Comparison of bundles grouping. ... 24
Figure 11 Numerical Parameters Not Affecting the Solution, Parameters Affecting the Solution ... 25
Figure 12 Vessel Nodalization Comparison. ... 26
Figure 13 Sensitivity Calculations on Decay Ratio for: Vessel, Time Step ... 27
Figure 14 Sensitivity Calculations on Frequency for: Vessel, Time Step 27 Figure 15 Final Model Results ... 29
Figure 16 Numerical Parameters Not Affecting the Solution, Parameters Used to Affect the Solution ... 30
Figure 17 All Cycles, Comparison of Initial and Final Model ... 31
Figure 18 Example of Decay Ratio spread after a PIE calculation ... 36
Figure 19 Uncertainty Calculations for Cycles 14, 15, 16, 17 ... 37
Figure 20 Averaged Spearman Correlation Coefficient for the Decay Ratio from all the points from Cycle 14. ... 39
Figure 21 Spearman Correlation for Cycle 14 Point 01 ... 47
Figure 22 Spearman Correlation for Cycle 14 Point 03 ... 47
Figure 23 Spearman Correlation for Cycle 14 Point 04 ... 48
Figure 24 Spearman Correlation for Cycle 14 Point 05 ... 48
vi
Figure 25 Spearman Correlation for Cycle 14 Point 06 ... 49
Figure 26 Spearman Correlation for Cycle 14 Point 08 ... 49
Figure 27 Spearman Correlation for Cycle 14 Point 09 ... 50
Figure 28 Spearman Correlation for Cycle 14 Point 10 ... 50
vii
List of Abbreviations
ASSUME Automated Software for Sensitivity and Uncertainty MEthod
BWR Boiling Water Reactor
DR Decay Ratio
FR Frequency
ITF Integral Test Facility
LOCA Loss Of Coolant Accident
NPP Nuclear Power Plant
PIE Propagation of Input Errors PDF Probability Density Function
PIRT Phenomena Identification and Ranking Tables SCRAM Safety Cut Rope Axe Man
SET Separate Effect Test
U.S. NRC United States Nuclear Regulatory Commission
viii
ix
Acknowledgments
I would like to gratefully thank to my supervisor Tomasz Kozlowski for the support and the advice during all my work, even during the weekends.
Also I would like to thank to Francesco Cadinu and to Kaspar Kööp for the irreplaceable help with Linux, which I really needed. Thanks to Sean Roshan for the smokes we had and to Roberta Hansson and Joanna Peltonen for the help with administrative issues. To all the people in the Nuclear Power Safety, Nuclear Reactor Technology, and the Reactor Physics divisions at KTH.
Last but not least I would like to thank to my girlfriend, Ilina Atanasova, for the support during all the time, and also to my friends and my family in Bulgaria which I love.
Thank you all!
This work has been performed thanks to the support of the Swedish Radiation Safety Authority (SSM).
x
1
1 INTRODUCTION
1.1 Boiling Water Reactor Stability
BWR plants are prone to instability due to the nature of the two-phase thermal-hydraulics in the reactor core facilitated further by a power feedback through core neutron kinetics, plant system thermal-hydraulics, and the plant control system. In general, BWRs may be subject to coupled thermal-hydraulic - neutronic instabilities while operating at relatively low flow rate relative to power (e.g. 30% of nominal flow and 50% of nominal power). The instability originates in the core and leads to oscillations in flow rate, pressure and void fraction.
A properly designed core should ensure that core oscillations are efficiently suppressed during nominal reactor operation. However, oscillations of local, regional or global nature may occur under a complex, and often unexpected combination of accident transient conditions. In a number of situations, SCRAM and consequent reactor power shutdown are the only tools available to stop the uncontrolled power excursions.
An example illustrating stability transient analysis is a real plant event that occurred in Oskarshamn-2 in 1999 (O-2 1999) feedwater transient.
A loss of feedwater preheaters and control system failure resulted in a situation with a high feedwater flow and low feedwater temperature without reactor SCRAM. In addition to the initiating event, interaction of automatic power and flow control system caused the plant to move into low flow - high power regime. Combination of the above events culminated in diverging power oscillations which triggered automatic scram at high power. The power evolution for the event is shown in Figure 1.
Figure
In order to get more data about this phenomenon, a stability benchmark has been specified by Vattenfall AB, Electricity Generation with leading author Tomas Lefvert [
code developers and analysts in the OECD Member countries to test their codes and also to validate the predictive capability of their respective codes and models for BWR stability analysis.
The measure for stability is u
Natural Frequency (FR) of the reactor power. The Decay Ratio is a measure of the growing or damping of the reactor power signal. If the value of the DR is between 0 and 1 then the oscillations are damping (the reactor is stable), see
oscillation are growing (the reactor is unstable)
Natural Frequency is the frequency of oscillation of the given reactor power signal.
On Figure 2, it can be se
the equation fit of the power signal (red line). The Decay Ratio is defined by equation (1):
2
Figure 1 Oskarshamn-2 1999 event
more data about this phenomenon, a stability benchmark has been specified by Vattenfall AB, Electricity Generation with leading author Tomas Lefvert [1]. The purpose of this benchmark was to enable code developers and analysts in the OECD Member countries to test their codes and also to validate the predictive capability of their respective codes and models for BWR stability analysis.
The measure for stability is usually the Decay Ratio (DR) and the Frequency (FR) of the reactor power. The Decay Ratio is a measure of the growing or damping of the reactor power signal. If the value of the DR is between 0 and 1 then the oscillations are damping , see Figure 2, while if the DR is greater than 1, the oscillation are growing (the reactor is unstable), see Figure
Natural Frequency is the frequency of oscillation of the given reactor
, it can be seen the power signal versus time (blue line) and the equation fit of the power signal (red line). The Decay Ratio is defined by equation (1):
+
=
2 3 1 2
2 1
B B B
DR B
more data about this phenomenon, a stability benchmark has been specified by Vattenfall AB, Electricity Generation with leading se of this benchmark was to enable code developers and analysts in the OECD Member countries to test their codes and also to validate the predictive capability of their
sually the Decay Ratio (DR) and the Frequency (FR) of the reactor power. The Decay Ratio is a measure of the growing or damping of the reactor power signal. If the value of the DR is between 0 and 1 then the oscillations are damping , while if the DR is greater than 1, the Figure 1. The Natural Frequency is the frequency of oscillation of the given reactor
en the power signal versus time (blue line) and the equation fit of the power signal (red line). The Decay Ratio is
(1)
3
Figure 2 Calculation of stability parameters – Decay Ratio and Frequency
1.2 System Codes
For the licensing of a NPP, a safety analysis is conducted which has to show that the proposed design can handle various transient scenarios, and to propose safe design solutions which would guarantee the safe operation of the reactor. All of these analyses need to be estimated using complex computer simulation tools. Those tools have been developed by collaboration between the research community, regulatory bodies, reactor constructors and the NPP operators. A number of complex physical equations and large amount of experimental data are combined in the so called codes to compute the reactor behavior in a transient and normal operation [2].
There is a large diversity of codes that simulate different physics, neutronics, thermal-hydraulics, structural mechanics, etc. One of the most important of these are the thermal-hydraulics system codes.
0 1 2 3 4 5 6 7 8 9 1 0
-1 -0 . 5 0 0 . 5 1 1 . 5 2
B3 B2
B1
Figure 3 Classification of models adopted in
Figure 3 shows the categorization of models adopted in the BWR stability analysis. Simplified Phenomenological Models are mainly addressing only a few of the basic phenomena involved in BWR stability or even considering most of them but using some simplifications at the mathematical or physical level. On the other hand, System
aiming at the detailed simulation of the BWR plant,
up-to date descriptions for each relevant phenomenon [
based on two-phase flow equations which are usually resolved in Eulerian-coordinates. The two
momentum, and energy conservation equations for the liquid and the vapor phases separately, and mass conservation equations for noncondesable gas present in the mixture [
Non
Simplified Phenomenological Models
Models for Basic Phenomena (heated
channel, loop)
Simplified BWR Models (neutron kintetics + TH
4
Classification of models adopted in the stability analysis of BWRs.
shows the categorization of models adopted in the BWR stability analysis. Simplified Phenomenological Models are mainly ssing only a few of the basic phenomena involved in BWR stability or even considering most of them but using some simplifications at the mathematical or physical level. On the other hand, System
aiming at the detailed simulation of the BWR plant, in some cases using to date descriptions for each relevant phenomenon [
phase flow equations which are usually resolved in rdinates. The two-phase flow is described by the mass, momentum, and energy conservation equations for the liquid and the vapor phases separately, and mass conservation equations for noncondesable gas present in the mixture [4].
Non-Linear Models for BWR Dynamics
Simplified Phenomenological Models
Models for Basic Phenomena (heated
channel, loop)
Simplified BWR Models (neutron kintetics + TH
feedback)
System Codes
Codes for Core Dynamics
Codes for Overall Plant
Dynamics
Engineering Models
Specific Codes for BWR Stability
General Purpose codes
he stability analysis of BWRs.
shows the categorization of models adopted in the BWR stability analysis. Simplified Phenomenological Models are mainly ssing only a few of the basic phenomena involved in BWR stability or even considering most of them but using some simplifications at the mathematical or physical level. On the other hand, System Codes are in some cases using to date descriptions for each relevant phenomenon [3]. They are phase flow equations which are usually resolved in phase flow is described by the mass, momentum, and energy conservation equations for the liquid and the vapor phases separately, and mass conservation equations for
Codes for Overall Plant
Dynamics
Engineering Models
Specific Codes for BWR Stability
General Purpose codes
5
System Codes have two classes that are considered in the safety analysis [3]:
1. Frequency Domain Codes, whose purpose is the linear stability analysis of BWRs or other boiling systems; they are based on linearization and Laplace transform of the governing equations.
2. Time Domain Codes, which include analysis tools specifically developed to simulate the transient behavior of plant systems;
these codes have the capability to deal with the non-linear features of BWRs and are based on the simulation techniques.
This work is mainly focused on general purpose time domain system codes (RELAP, TRACE, etc.) coupled with neutronic feedback coming from neutronic codes (PARCS, DYN3D, etc.).
1.3 State-of-the-Art in BWR Stability
Due to the numerical approximations and the empirical nature of the included models in the thermal-hydraulic system codes, a huge effort has been done aiming at validation of the codes. The validation is done from NPPs data and special scaled down test facilities, Integral Test Facilities (ITF) and Separate Effects Test Facilities (SET). SET are experiments where a particular phenomenon is being studied on its own, while ITF are more complex experiments which study the evolution of a system under transient conditions. The literature review in this thesis mainly shows the validation of the codes for stability predictions against real plants, and in some works sensitivity calculations have been performed for parameters that affect it.
In work [5], a BWR stability validation has been performed for the transient system codes RAMONA-3, RAMONA-5, SIMULATE-3K, and POLCA-T. The performance of these codes, against the stability
6
measurements performed in cycle 19 at the Swiss NPP Leibstadt (KKL), a BWR/6 from General Electric, showed that all results have good agreement with the measured values of Decay Ratio and Frequency, within the measurement uncertainties. However, the analysis showed also that the codes consistently underestimate the Frequencies. The same Frequency under-prediction has also been observed in [6], where the same plant, cycles 19 and 23, has been validated against the system code SIMULATE-3K.
Work [7], summarizes part of the extensive validation database (Ringhals-1, cycles 14-17; Oskarshamn-3, cycles 7-17; Olkiluoto-1, cycles 18-26; Leibstadt, cycle 19) for the code SIMULATE-3K, and discusses the influence of fuel pin model parameters on the stability results. It is claimed that the Decay Ratios and Frequencies are well predicted, however, in this work the Frequency prediction is overestimated for Ringhals-1, Oskarshamn-3, and underestimated for Leibstadt. Sensitivity calculations are also presented for the fuel pin parameters, but no quantification of these is presented.
Work [8], presents the results of the time domain coupled 3D neutron kinetics/thermal-hydraulic code POLCA-T analysis of BWR Peach Bottom 2 End of Cycle 2 Low-Flow Stability Tests. The Decay Ratios calculated by the coupled code are in good agreement with the measured data. However, there is a systematic underestimation of the calculated Frequencies. Further, the paper reports a sensitivity study for 4 parameters, which shows that core inlet temperature, carry under, and gas gap heat transfer coefficient have stronger effect on the results, while the steam dome pressure has negligible influence on the stability.
In work [3], parameters affecting the stability have been identified and characterized. Four broad categories have been identified:
7 1. Thermal-hydraulic feedback 2. Neutronic feedback
3. Simultaneously neutronic and thermal-hydraulic feedback 4. Specific pattern of oscillations – (regional or out-of-phase).
No clear quantification of influence of the above has been shown in this work.
Paper [9], describes the validation of TRACE/PARCS for BWR Ringhals-1 stability analysis. Sensitivities are also performed to investigate the impact on the Decay Ratio and Natural Frequency of the gas gap heat transfer coefficient, and the pressure loss coefficients of the spacers. Results show good agreement with measured stability parameters (DR and FR), but again slight under-prediction for the Frequency. Sensitivities showed that the spacers pressure losses impact the Decay Ratio only, while the gas gap heat transfer coefficient impacts both stability parameters.
1.4 Overview of Uncertainty Methodologies
One of the most important terms in reactor safety is the safety margin.
The illustration of the safety margin is shown on Figure 4. Safety margin is the reserve that certain value has before crossing the regulatory acceptance criteria. Licensing of NPPs has been done by assuming conservative values for the computations. Nowadays it is possible to estimate certain parameters using non-conservative data with the complement of uncertainty evaluation, and these calculations can also be used for licensing. It can be seen on the figure that the best-estimate plus uncertainty value can be slightly lower than the conservative approach. Some countries, like Sweden, use this difference for power up- rates and life extension of nuclear power plants. Thus, uncertainty methods within the nuclear field are important to help increase production, while staying below the safety limits [10]. This section
8
discusses sources of uncertainties, various uncertainty analyses methods, and describes the choice of uncertainty method appropriate for this work.
Figure 4 Scheme of the Safety Margin [10]
Several methods have been developed for the uncertainty evaluation of system codes:
Adjoint Sensitivity Analysis Procedure [11]: The theory behind the method is to get sensitivity response for all uncertain parameters and then use it to assess the uncertainty by propagation of error formulas. The method works for both linear and non-linear systems and this makes it appropriate for the case of a Nuclear Power Plant. This method requires significant modifications of the code, and this is out of the scope of this work. That is why this method is not considered further.
9
Propagation of input errors [12]: This approach represents statistical variation of the input parameters, together with their uncertainties, in order to reveal the propagation of errors through the code, Figure 5 (a).
The method relies on a real code calculation, but instead of calculating certain parameters with discrete values, it assumes the given parameter as a probability distribution by performing the calculation a number of times. This approach allows evaluation of correlation coefficients between the uncertain parameters and the target output parameter, which can be used to determine the most influential parameters. The drawback of this method is that all uncertain parameters have to be identified and probability density functions must be assigned to them, which sometimes have to be an engineering judgment. Another issue is that the error propagates through the code, which is not a perfect simulation tool.
Propagation of output errors [13]: This methodology focuses on extracting the errors from a suitable experimental database coming from integral test facilities and real plant data. The calculated results from the system codes are extrapolated with the uncertainty band, to obtain the desired uncertainty in the code prediction, Figure 5 (b). The benefit is that the method does not require the identification of all uncertain input parameters. The drawback is that it relies on non-physical extrapolation of errors coming from all available relevant data and giving the final uncertainty band which is not necessarily directly related to the given case.
The second method (Propagation of Input Errors) is the most used in applications within the industry and it has been extensively verified. It provides the opportunity to identify the influence of each parameter on the results, while this is not possible in the third (Propagation of Output Errors) method. Furthermore, for the propagation of output errors, one needs experimental data to implement the method, which in the case of BWR stability is not available. Moreover, the propagation of input errors
10
(PIE) method relies on actual code calculations without using approximations such as fitted response surfaces. For those reasons, the PIE method has been chosen in the work.
Figure 5 (a). Propagation of input errors [12] (b). Propagation of output errors [12].
1.4.1 Propagation of Input Errors Methodology
The Propagation of Input Errors method (PIE) is based on statistical variation of input parameters and it does not require modifications of the code (with the possible exception of the closure models and correlations). In principle there is no limitation on the number of uncertain parameters and range of variation of the perturbed parameters. In this method, every parameter is considered as a probability distribution, not as a discrete value, as shown in Figure 6. A random set of discrete values based on their Probability Distribution Functions (PDFs), for all “suspicious” parameters needs to be generated. The sources of PDFs are literature (e.g. scientific papers, PIRTs) or an expert judgment. After all important parameters are identified and the random sample is set, N model runs should be performed with these parameter uncertainties as input. The required number of runs (N) is given by the Wilks’ formula, Eq. (2) and (3) [12]:
11
100 / ) 100 / (
1− α N ≥ β (2)
One can be β% confident that at least α% of the combined influence of all the characterized uncertainties is below the tolerance limit.
The formula for the two sided statistical tolerance intervals is given as:
100 / )
100 / ))(
100 / ( 1 ( ) 100 / (
1− α N −N − α α N−1≥β (3)
Here α and β have the same meaning as in the previous equation.
Figure 6 Illustration of the PIE method and results interpretation [12].
Depending on the requirements of the confidence interval, one can estimate how many calculations are needed in order to obtain the desired results. The minimum number of runs can be found in Table I. The U.S.
NRC has recognized the 95%/95% coverage/confidence estimation, and this is usually taken as the best-estimate prediction requirement [12].
In this work a set of 93 runs will be used.
12
Table I Minimum number of calculations N.
One-sided statistical tolerance limits
Two-sided statistical tolerance limits
α
β 0.9 0.95 0.99 0.9 0.95 0.99
0.9 22 45 230 38 77 388
0.95 29 59 299 46 93 473
0.99 44 90 459 64 130 662
In the PIE calculation, all parameters are varied simultaneously. In order to achieve 95/95 coverage/confidence, a total of 93 PIE runs have to be performed as specified in Table I. The PIE Uncertainty calculation provides the uncertainty of the output parameters which are to be examined. Comparing this result with the measurements provides a measure of the confidence of the code prediction.
The Spearman rank correlation coefficient between the output parameter and the input parameter value can be used to identify the most influential parameters. This coefficient provides a quantitative indication if the considered parameter (and its uncertainty, represented by its PDF) has a significant impact on the stability of the reactor. The criterion for “importance” consists of the so called critical value of the Spearman correlation coefficient. If the absolute value of the estimated rank is higher than the critical value, then this parameter is determined to have a statistically significant impact on the result. Otherwise, the parameter will be considered to have low or no importance. The critical values rs for the Spearman rank correlation coefficient are evaluated using the formula (4) [14].
−1
±
= n
rs z
(4)
13
For the confidence interval of 95%, z=1.96 and n is the number of runs (93 in this case) [14]. Therefore, one gets that the Critical values for the Spearman rank correlation coefficient are rs=±0.2.
1.4.2 Sources of Uncertainties
The propagation of input errors method relies on the availability of the input uncertainties Probability Density Functions. Without these functions no PIE uncertainty calculation could be done. Several classifications of uncertainties can be found. In general, uncertainties can be divided in two broad categories.
The first category is called aleatory or stochastic uncertainty. This uncertainty comes from the random nature of how a system would behave. An example of this uncertainty can be the measurement of the temperature in a room. Each time one measures, it will be different within the limits of the measuring instrument. This uncertainty is irreducible (cannot be reduced using the same technology) and is associated with the performance of an experiment or a measurement [15].
The second category is the subjective or epistemic uncertainty. This category stands for the lack of knowledge of a certain parameter and one is not able to specify the exact value of a quantity that is assumed to have a fixed value within a particular analysis and thus is a property of the analysts carrying out the study. An example of this uncertainty can be a certain correlation in the particular analytical model which might not be able to represent the reality in some ranges of other parameters during the simulation. This uncertainty is reducible (it could be reduced if new knowledge is provided) [15].
Three major sources of uncertainties in accident analysis have been identified [16].
14
Code or model uncertainty: Uncertainty associated with the models and correlations, the solution scheme, model options, not modeled processes, data libraries and deficiencies of the computer program.
Representation or simulation uncertainty: Uncertainty in representing or idealizing the real plant, such as that due to the inability to model a complex geometry accurately, three dimensional effects, scaling, control and system simplifications.
Plant uncertainty: Uncertainty in the measuring or monitoring of a real plant, such as reference plant parameters, instrument error, set points, instrument response.
More detailed identification of uncertainties has been done in [17] and categorized in [10] and [13]:
A. Balance (or conservation) equations are approximate [10] and [13].
B. Presence of different fields of the same phase, for example liquid droplets and film thus only one velocity per phase is considered by the codes [10] and [13].
C. Geometry averaging at a cross-section scale - the lack of consideration of the velocity profile (i.e. cross-section averaging) constitutes an uncertainty source of ‘geometric origin’ [10] and [13].
D. Geometry averaging at a volume scale - only one velocity vector (per phase) is associated with a hydraulic mesh along its axis while different velocity vectors may occur in reality, [10] and [13].
E. Presence of large and small vortices or eddies - energy and momentum dissipation associated with vortices are not directly accounted for in the equations at the basis of the codes [10] and [13].
15
F. The second principle of thermodynamics is not necessarily fulfilled by the codes - irreversible processes occur as a consequence of accidents in nuclear reactor systems [10] and [13].
G. Models of current interest for thermal-hydraulic system codes comprise a set of partial differential equations which are solved by approximate numerical methods [10] and [13]. Truncation errors in the codes, together with the numerical calculation of a non-linear system, might lead to huge magnification of the output error.
H. Extensive and unavoidable use is made of empirical correlations [10] and [13].
I. A paradox should be noted: a steady state and fully developed flow condition is a necessary prerequisite or condition adopted when deriving correlations. In other words, all qualified correlations must be obtained under steady state and fully developed flow conditions. However, almost in no region of the nuclear power plant do these conditions apply during the course of an accident [10] and [13].
J. The state and the material properties are approximate [10] and [13].
K.Code user effect exists - different groups of users using the same code and with the same information available for modeling a nuclear power plant do not achieve the same results [10] and [13]:
L. Computer and compiler effect exists - the code is used on different computational platforms and current experience has shown that the same code with the same input deck applied within two computational platforms produces different results. Differences are typically small in ‘smoothly running transients’, but may become noticeable in the case of threshold or bifurcation driven transients [10] and [13].
M. Imperfect knowledge of boundary and initial conditions. Some boundary and initial condition values are unknown or only approximately known [10] and [13]. Local measurements from NPPs are usually not provided by the plants – code users rarely get
16
local measurements to compare their calculation with. Not having this data could give a good overall output which is the result of compensation of errors.
N. Code model deficiencies cannot be excluded. Examples of the phenomena for which some thermal-hydraulic system code models prove deficient are:
a. The heat transfer between the free liquid surface and the upper gas–steam space;
b. The heat transfer between a hotter wall and the cold liquid flowing down inside a steam–gas filled region.
c. These deficiencies are expected to have an importance only in special transient situations [10] and [13].
O.The Nodalization effect exists. The nodalization is the result of a wide range brainstorming process, where user expertise, computer power, and manual play a role [13].
All these categorizations are shown in Figure 7.
Figure 7 Epistemic, Subjective
(Reducible)
Code or model uncertainties
A. Balance equations are approximate
B. Presence of different fields of the same phase
C. Geometry averaging at a cross-section scale
D. Geometry averaging at a volume scale
E. Presence of large and small vortices or eddies
F. The second principle of the thermodynamics is not
necessarily fulfilled by codes
G. The numerical solutions are approximate
H. Use of empirical correlations
I. Correlations to be obtained during steady- state and fully developed
flow conditions.
N. Code model deficiencies cannot be excluded.
17
Diagram of the sources of uncertainties
Uncertainties
Epistemic, Subjective
the thermodynamics is not
G. The numerical solutions
state and fully developed
N. Code model deficiencies
Representation or simulation uncertainties
J. The state and the material properties are
approximate
K. Code user effect exists.
L. Computer and compiler effect exists.
O. Nodalization effect exists
Aleatory, Stochastic (Irreducible)
Plant Uncertainties
M. Imperfect knowledge of boundary
and initial conditions.
Aleatory, Stochastic (Irreducible)
Plant Uncertainties
M. Imperfect knowledge of boundary
and initial conditions.
18
1.5 Aim of the Thesis
In order to use the system codes, one must represent the simulation of a given case (reactor or a plant) in a way that the code will “understand”
what it is supposed to simulate. Part of this representation (explained in detail in Section 1.4.2) is the nodalization, for which there are guidelines, but no clear rules one could follow in order to represent the model to the simulation tool. In this thesis a numerical sensitivity study was performed to verify that the chosen nodalization and time step were appropriate.
Nowadays, Best-Estimate codes plus uncertainty evaluation have been used for Loss of Coolant Accidents (LOCA) predictions and other transient scenarios. However, no uncertainty calculations have been done for Boiling Water Reactor (BWR) Stability Analysis. This thesis is focused on the implementation of an uncertainty method and evaluation of output uncertainties for BWR Stability Analysis. Performing an uncertainty analysis for BWR stability would give better understating of the phenomenon and it would help verify and validate (V&V) the codes used to predict NPPs.
Another issue which is addressed in this thesis is the ranking of parameters which affect the stability. In the literature review of BWR Stability Analysis (Section 1.3), various works showed sensitivity results derived from various cases and parameters were identified and characterized. However, in none of these works a quantitative evaluation of the influential parameters has been done. This current work will use the uncertainty method to rank parameters according to their influence on the stability.
19
2 BWR STABILITY TEST CASE
The Ringhals-1 OECD Stability Benchmark was used in this work, primarily because of the large database of available data, and the lack of scaling effect. The measurements for the DR and the FR are provided for four cycles, with 8 to 11 points per cycle, for a total of 37 points at different power and flow conditions. In addition to the measurements, the uncertainty of the evaluation of the Decay Ratio is also provided in the benchmark. However, this is the uncertainty due to the evaluation method only and it does not include the measurement uncertainty [1].
2.1 TRACE/PARCS Model Description
In this thesis, the time domain coupled code TRACE/PARCS was used to simulate the desired stability scenario. The simulation consists of achieving a Steady-State solution for the desired power and flow condition, and performing a control rod perturbation transient. After the perturbation, the natural response of the reactor has an oscillatory behavior and the Decay Ratio and the Frequency of Oscillation can be evaluated. The stability parameters (DR and FR) are evaluated using the DRARMAX program, which is an ARMA-based method, developed at Purdue University [18], [19], [20].
The TRACE/PARCS model consists of Reactor Pressure Vessel, Fuel Bundles, Steam Separators, Steam Lines, Recirculation Loop, and a Feed Water Pipe. Nodalization diagram is shown on Figure 8. In addition to the thermal-hydraulic model, a neutronics feedback is given by the PARCS code.
20
Figure 8 Initial Nodalization Scheme
The initial results were obtained using the coupled system code TRACE/PARCS without any modification of the model and are presented in Figure 9 and Table II for all cycles [9], [21], [22]. The measured uncertainties are provided in the benchmark, which are attributable to the method used to evaluate the Decay Ratios from the power signal of the Ringhals-1 reactor [1]. For Cycle 14, one can notice that all points, except one, fall in the region of the measurement, within the given uncertainty. The Frequency prediction is off by about 10%, but the trend is obviously correct. For the rest cycles the prediction is acceptable.
21
Table II Measured and Calculated Decay Ratio and Frequency - all Cycles.
Cycle Point
Decay Ratio
Measured Error Calculated Abs. Difference
Frequency Measured Calculated
C14P01 0.30 0.12 0.34 0.04 0.43 0.39
C14P03 0.69 0.06 0.60 0.09 0.43 0.40
C14P04 0.79 0.05 0.74 0.05 0.55 0.49
C14P05 0.67 0.06 0.71 0.04 0.51 0.47
C14P06 0.64 0.07 0.61 0.03 0.52 0.47
C14P08 0.78 0.05 0.80 0.02 0.52 0.47
C14P09 0.80 0.05 0.80 0.00 0.56 0.51
C14P10 0.71 0.06 0.69 0.02 0.50 0.46
C15P01 0.23 0.15 0.43 0.20 0.44 0.39
C15P02 0.24 0.15 0.46 0.22 0.42 0.37
C15P03 0.21 0.15 0.43 0.22 0.43 0.38
C15P04 0.33 0.10 0.38 0.05 0.44 0.40
C15P05 0.43 0.09 0.46 0.03 0.44 0.40
C15P06 0.59 0.07 0.61 0.02 0.47 0.42
C15P08 0.77 0.05 0.72 0.05 0.55 0.49
C15P09 0.67 0.06 0.81 0.14 0.53 0.50
C15P10 0.60 0.07 0.68 0.08 0.54 0.50
C16P01 0.54 0.07 0.60 0.06 0.48 0.43
C16P02 0.54 0.07 0.65 0.11 0.48 0.43
C16P03 0.69 0.06 0.74 0.05 0.47 0.43
C16P04 0.71 0.06 0.60 0.11 0.52 0.46
C16P05 0.67 0.06 0.70 0.03 0.49 0.45
C16P06 0.79 0.05 0.83 0.04 0.49 0.43
C16P07 0.72 0.06 0.71 0.01 0.50 0.45
C16P08 0.82 0.05 0.78 0.04 0.49 0.44
C16P09 0.87 0.04 0.90 0.03 0.48 0.43
C16P10 0.65 0.06 0.69 0.04 0.50 0.44
C16P11 0.66 0.06 0.81 0.15 0.48 0.45
C17P02 0.24 0.15 0.32 0.08 0.46 0.36
C17P03 0.22 0.15 0.38 0.16 0.44 0.36
C17P04 0.32 0.10 0.32 0.00 0.46 0.37
C17P05 0.28 0.10 0.35 0.07 0.42 0.37
C17P06 0.34 0.09 0.41 0.07 0.46 0.37
C17P07 0.33 0.09 0.38 0.05 0.46 0.40
C17P08 0.41 0.09 0.42 0.01 0.48 0.40
C17P09 0.57 0.07 0.52 0.05 0.47 0.40
C17P10 0.49 0.08 0.44 0.05 0.49 0.42
22
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
(a) Ringhals Cycle C14
Calculated
Measured
Decay Ratio Frequency
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Decay Ratio Frequency
Calculated
Measured
(b) Ringhals Cycle C15
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Decay Ratio Frequency
Calculated
Measured (c) Ringhals Cycle 16
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Decay Ratio Frequency
Calculated
Measured
(d) Ringhals Cycle 17
Figure 9 Results from the base cases. Cycles (a) 14, (b) 15, (c) 16, (d) 17.
23
2.2 Model Validation
Prior to the application of the uncertainty method, sensitivity study was performed for the Numerical parameters for the given TRACE/PARCS model in order to identify if Nodalization scheme, Time Step, and Bundles Grouping are appropriate.
The validation was done by refining all possible components and obtaining a point where Nodalization refinement doesn’t change the solution i.e. space-time converged solution was obtained. Given the initial model, Figure 8, there are certain numerical parameters that can be refined, that are listed in Table III. The third column shows what refinement levels were evaluated. Only Cycle 14 has been tested for the numerical discretization.
Table III Numerical Uncertain Parameters
Parameter Initial Model Refined Model
Bundles Grouping Parameters Bundles Grouping Symmetric ~ 324
Bundles
648 Bundles
(One Channel per Fuel Assembly) Nodalization Parameters
Recirculation Loop
Nodalization 8 Cells 16 Cells
Feed Water Pipe
Nodalization 3 Cells 6 Cells
Steam Line Pipe
Nodalization 1 Cell 2; 4 Cells
Separators Nodalization 5 Cells 9 Cells
Vessel Nodalization 11 Levels 20; 29; 38 Levels Bundles Nodalization 28 Cells 14 Cells; 56 Cells
Time Step Parameters
Time Step 0.04 s 0.01; 0.004; 0.001; 0.0008
24
In order to be sure that the model is space-time converged, all the parameters in Table III were refined. If there is a change in the solution with the refinement, then the model is not space-time converged. This means that all parameters have to be refined, until the solution is very close (or the same) as the previous refined level.
The initial model consists of 324 fuel bundles arranged in rotational symmetry in order to get the desired 648 fuel assemblies in the Ringhals- 1 reactor. This means that one channel is used to represent two physical bundles. When one applies the “One Channel per Fuel Assembly”
model, which consists of 648 different bundles, the calculation time doubles and the result is practically the same, see Figure 10. Therefore, this is a good grouping scheme, which leads to reducing the time of calculation by a factor of 2 without noticeable impact on the result.
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
DR Base Case DR 648 FR Base Case FR 648
Calculated
Measured
Figure 10 Comparison of bundles grouping.
Results show that the next four parameters (Recirculation Loop, Feed Water Pipe, Steam Lines, and Separators) are already space-time converged, see Figure 11 (a). Any of them refined have no, or very small, impact on the final solution. The only exception is Point 01 of Cycle 14, which has the lowest Decay Ratio. In this point, there is a very
25
small change in the solution, but because Decay Ratio is low, this can be neglected.
The parameters that showed significant changes in the solution when they are refined are Bundles and Vessel Nodalizations, and Time Step, see Figure 11 (b). One can notice that the Frequency has a good agreement with the initial solution, but the Decay Ratio has changed significantly. Refinement of Bundles gives a positive correlation for the DR, while Refinement of Vessel Nodalization and Time Step give negative correlation.
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
DR Base Case
DR Recirculation Loop Refined x2 DR Feed Water Pipe Refined x2 DR Steam Lines Refined x2 DR Separators Refined x2 FR Base Case
FR Recirculation Loop Refined FR Feed Water Pipe Refined FR Steam Lines Refined x2 FR Separators
Calculated
Measured
(a)
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
DR Base Case DR Bundles Refined x2 DR Vessel x2
DR Time Step Refined x5 FR Base Case FR Bundles Refined x2 FR Vessel x2
FR Time Step Refined x5
Calculated
Measured
(b)
Figure 11 (a) Numerical Parameters Not Affecting the Solution (b) Parameters Affecting the Solution
A proper evaluation of those phenomena was done by performing space-time convergence studies. Refinement of the Bundles could be done only by a factor of two due to coupled code limitations. In this case, no further refinement could be obtained from the code and because of the difference between the initial and the refined solution, this concludes that the refined solution is space-time converged.
26
Figure 13 and Figure 14 show convergence studies for the Vessel and Time Step refinement. In these figures it is obvious that the space-time converged solution has not been achieved with the initial model and there are numerical errors in both parameters (Decay Ratio and Frequency).
For the Vessel, a Nodalization finer by a factor of 3 was chosen since the differences in refinements by 3 and by 4 are negligible.
On Figure 12 one can see how the Vessel has been renodalized. Vessel Level 7 and Level 8 have not been renodalized because this is the part where the separators are connected and the code requires a vessel nodalization which is consistent with the separator.
Figure 12 Vessel Nodalization Comparison.
Furthermore, 10 times smaller time step (0.004s) was chosen. This is already small enough to achieve time-converged solution, as shown on Figure 13 (b) and Figure 14 (b).
27
x1 x2 x3 x4
0.2 0.3 0.4 0.5 0.6 0.7 0.8
(a) Vessel Nodalization
DR C14P01 DR C14P06 DR C14P09
Decay Ratio
Vessel Nodalization Factor
1E-3 0.01 0.1
0.2 0.3 0.4 0.5 0.6 0.7 0.8
(b) Time Step
DR C14P01 DR C14P06 DR C14P09
Decay Ratio
Time Step (s)
Figure 13 Sensitivity Calculations on Decay Ratio for: (a) Vessel (b) Time Step
x1 x2 x3 x4
0.35 0.40 0.45 0.50 0.55
(a) Vessel Nodalization
FR C14P01 FR C14P06 FR C14P09
Frequency
Vessel Nodalization Factor
1E-3 0.01 0.1
0.35 0.40 0.45 0.50 0.55
(b) Time Step
FR C14P01 FR C14P06 FR C14P09
Frequency
Time Step (s)
Figure 14 Sensitivity Calculations on Frequency for: (a) Vessel (b) Time Step
28
The final space-time converged model was achieved with the following parameters:
• Vessel refined x3, 29 nodes
• Bundles refined x2, 56 nodes
• Time Step refined x10, 0.004 seconds
In addition, the gas gap heat transfer coefficient had to be changed in order to fit the measurements from the Ringhals-1 Benchmark. The initial value was 4800 W/(m2K), and the new value was set to 9500 W/(m2K). In work [23, Figure 14], one can see the Gap Conductance as a function of the burn up for the bundle of type BWR 8x8 which is the fuel type used in Ringhals-1. The trend is that the Gap Conductivity increases with the burn up and it starts from 4000 W/(m2K) and increases up to 10220 W/(m2K). Given that the uncertainty of this parameter is very large (see Table VI), about 35%, the Gap Conductance used in this calculation, 9500 W/(m2K), is still realistic. This value was set for Cycle 14 and used for all other cycles.
On Figure 15, the new space-time converged base case (red circle and green triangle) for Cycle 14 is shown in comparison with the initial one (black rectangle and blue triangle). One can see that the Decay Ratio (DR) matches better in comparison with the initial, with the exception of the point with the lowest DR. This could be explained by the fact that this point always has different behavior in relation to the rest points.
In [21] and [22], this point has opposite Spearman Correlation Coefficient for some parameters than the other points, and in [22] the calculated DR uncertainty of this point is 6 times larger than the others.
Moreover, the Frequency of Oscillation (FR) matches the measurements of Ringhals-1 Benchmark very well, which gives more confidence in the new model. For those reasons, it is suspected that Point 01 of Cycle 14 could have inaccuracy in measurements or wrong specification data.
29
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
DR Initial Base Case DR New Base Case FR Initial Base Case FR New Base Case
Calculated
Measured
Comparison of Initial and Final Nodalization Cycle 14
Figure 15 Final Model Results
Further, a study was performed to show that the new base case model is space-time converged by refining all the components again, with the exception of the Bundles due to the code limitations.
The first group of parameters which did not have large effect on the solution was Recirculation Loop, Feed Water Pipe, Steam Lines, and Separators. With the new base case, those have similar negligible behavior see Figure 16 (a).
The second group of parameters was Bundles and Vessel Nodalizations, and Time Step. Bundles Nodalization scheme has reached 56 cells for a 4 meter assembly, the distribution of the cells was non-uniform, and the smallest node length was 0.029 m. This gives confidence that the number of cells and their length are enough for this type of calculation.
The next two parameters Vessel and Time Step are shown on Figure 16(b). Vessel refinement x4 and Time Step refinement x250 means refinement factor based on the initial model. It can be seen that the differences between the final and refined solutions are small. This means that the new model is Space-Time Converged.
30
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
DR Final Nodalization DR Recirculation Loop Refined x2 DR Feed Water Pipe Refined x2 DR Steam Lines Refined x2 DR Separators Refined x2 FR Final Nodalization FR Recirculation Loop Refined x2 FR Feed Water Pipe Refined x2 FR Steam Lines Refined x2 FR Separators Refined x2
Calculated
Measured (a)
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 0.0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
DR Final Nodalization DR VesselRefined x4 DR Time Step Refined x250 FR Final Nodalization FR VesselRefined x4 FR Time Step Refined x250
Calculated
Measured (b)
Figure 16 (a) Numerical Parameters Not Affecting the Solution (b) Parameters Used to Affect the Solution
When Cycle 14 of the Ringhals-1 Benchmark was set, the other cycles (C15, C16, and C17) were run with the new Nodalization Scheme and the new gas gap heat transfer coefficient. Comparison with the initial model is shown in Figure 17. In all cases the Final Model’s Frequency matches better with the measured trend. Moreover, the Decay Ratio in all cycles is closer to the diagonal line of perfect agreement. Especially for Cycle 17, the Decay Ratio has a much better trend than the initial model. The new space-time refined model will be used as the base case for all subsequent calculations.
In order to quantify the accuracy of the new model, the standard deviations of the Decay Ratio and the Frequency of all the points of all cycles have been calculated, see Table IV. One can see that the Standard Deviations for the DR are lower for the Final (space-time converged) Model, moreover the space-time converged model’s Frequency matches much better than the initial. This proves that the model accuracy has been improved, and the prediction error has been reduced with the space-time converged model.
