# Analyses of Instability Events in the Peach Bottom-2 BWR Using

Authors: Antonella Lombardi Costa; Alessandro Petruzzi; Francesco D'Auria; Walter Ambrosini
Journal: Science and Technology of Nuclear Installations
 (2008)

---

## Abstract

Boiling water reactor (BWR) instabilities may occur when, starting from a stable
operating condition, changes in system parameters bring the reactor towards an
unstable region. In order to design more stable and safer core configurations,
experimental and theoretical studies about BWR stability have been performed to
characterise the phenomenon and to predict the conditions for its occurrence. In
this work, contributions to the study of BWR instability phenomena are
presented. The RELAP5/MOD3.3 thermal-hydraulic (TH) system code and the
PARCS-2.4 3D neutron kinetic (NK) code were coupled to simulate BWR transients.
Different algorithms were used to calculate the decay ratio (DR) and the natural
frequency (NF) from the power oscillation predicted by the transient
calculations as two typical parameters used to provide a quantitative
description of instabilities. The validation of the code model set up for the
Peach Bottom Unit 2 BWR plant is performed against low-flow stability tests
(LFSTs). The four series of LFST have been performed during the first quarter of
1977 at the end of cycle 2 in Pennsylvania. The tests were intended to measure
the reactor core stability margins at the limiting conditions used in design and
safety analyses.


---

## Body

## 1. Introduction

In the last four decades, the nuclear power industry has been upgrading and
developing light water reactor technology, and preparing to meet the future demand
for energy. The presently operating BWRs contribute with about 21% of the total
produced nuclear power worldwide. These plants have reached very ambitious goals of
safety and reliability, together with high availability factors, notwithstanding the
flow instability and thermal-hydraulic oscillations that may affect BWRs under
particular operating conditions.The stability of BWR systems has been of great concern from the safety and the design
point of view at the beginning of the nuclear era; nowadays, the design of reactors
having appropriate stability margins, the adoption of operating procedures avoiding
possible unstable regions, and the development of mitigation strategies to cope with
inadvertent instability occurrences have strongly limited safety concerns in this
regard. This is a direct consequence of the large operating experience gained with
BWRs and of the increased knowledge of instability phenomena obtained from both
experimental and computational activities aimed at simulating reactor behaviour.BWR instabilities occur when an operating condition becomes unstable after some
change in system parameters. As a consequence, state variables identifying the
reactor working conditions are observed to oscillate in different ways depending on
the modalities of the departure from the stable operating point.Figure1 shows an example of the power-flow
map for a BWR. The lower right side of the plot marks the allowed operating region,
the grey region can only be entered if special measures are taken, and finally, the
black regime is forbidden due to stability concerns. In the figure, the 100%
rod line, also called 100% flow-control line, is a line on the BWR power-flow
map which passes through the normal operating point of the BWR (100% power and
100% flow rate). The control rods do not move on this line. Therefore,
operating points on the line are identified by a fixed reactivity value. BWRs
operate on this line during startup and shutdown operations by flow control via the
recirculation pumps [1].Figure 1
Instability region in the power-flow map for a typical BWR.Conditions corresponding to about 50% core nominal power and 30% core inlet
flow rate may be seen as the area in the power-flow map where the highest
probability of oscillations occurs. The operation in this area is avoided by means
of adequately defined control and trip conditions. Nevertheless, certain perturbed
transient conditions can still result in time windows, in which operation in this
area will occur, usually accompanied with the observation of oscillating core
behavior. The core two-phase flow itself provides a potential for oscillatory
behavior and the strong feedback between moderator coolant density and core power
enhances the effect under certain conditions. In-phase and out-of-phase power
oscillations have been actually observed and both modes, for large amplitudes, can
have an unwanted influence on the fuel integrity. Control systems for the RPV
(reactor pressure vessel) pressure and downcomer levels can influence the
oscillatory behavior in unfavourable way. Therefore, the main aims of BWR stability
analyses could be summarized as follows [2]: (i)
to assess the stability margins in the reactor plant, including normal
and off-normal conditions;(ii)
to predict the transient behavior of the reactor, an unstable condition
should occur;(iii)
to help in designing and assessing the effectiveness of countermeasures
adopted to prevent and mitigate the consequences of instabilities.


## 2. Coupled Neutronic/Thermal-Hydraulic Instabilities

BWR transient scenarios, that involve considerable reactivity changes, are described,
for example, in [2, 10]. These documents address overpressurisation events, large
break loss of coolant accidents (LBLOCAs), feedwater temperature decrease, increase
of core flow, main circulation pump flow rate increase, anticipated transient
without scram (ATWS), turbine trip (TT), and control rod removal.From the point of view of the BWR safety, the most important type of power
instability is the reactivity oscillations excited by thermal-hydraulic mechanisms.
Two types of instability by reactivity have been characterized.(i)
In-phase (core-wide) instability. In this case, all the
variables (power, mass flow, pressure, etc.) oscillate in phase
determining a limit cycle; from the point of view of safety, this type
of instability has relatively small relevance, unless it is associated
with an ATWS.(ii)
Out-of-phase instability. In this case, the
instabilities occur when a neutronic azimuthal mode is excited by
thermal-hydraulic mechanisms causing asymmetric power oscillations; at a
given time, while part of the reactor operates at high-mass flow and
low-power level, in the other part the opposite happens; this behaviour
must be studied in detail because of safety implications.A very complex type of power instability in BWRs consists of out-of-phase regional
oscillations, in which, normally, subcritical neutronic modes are excited by
thermal-hydraulic feedback mechanisms. The out-of-phase mode of oscillation is a
very challenging type of instability and its study is relevant because of the safety
implications related to the capability to promptly detect any such inadvertent
occurrence by in-core neutron detectors, thus triggering the necessary
countermeasures in terms of selected rod insertion or even reactor shutdown.Power oscillations can, for large amplitudes, have an unwanted influence on the fuel
integrity. In the fuel temperature limitation, it is essential to prevent the
exceeding of the melting point (3073.15 K forUO2). Fuel elements subjected to temperatures sufficiently high to induce
centreline melting will experience a significantly higher probability of failure
(loss in the functional behavior caused by a change in the physical properties).
Furthermore, the low thermal conductivity of ceramic fuels leads to high temperature
gradients that can cause fuel cracking and swelling [12].


## 3. Inadvertent and Induced BWR Instabilities

Some of the several occurred instability events in BWR plants were inadvertent and
other ones were induced intentionally as experiments. These instabilities were
identified as periodic oscillations of the neutron flux via instrumentation
readings. Essentially, neutronic power signals from local power range monitors
(LPRMs) and average power range monitors (APRM) have been used to detect and study
the power oscillations.One of the inadvertent instability events happened in LaSalle NPP, in 1988. During a
routine surveillance test, an instrument technician inadvertently caused the
shutdown of both recirculation pumps. As a consequence, the core flow rate was
rapidly reduced from 76% to 29% of the rated value, corresponding to
natural circulation conditions; this, in turn, led to the isolation of some of the
steam extraction lines leading to the pre-heaters. The result of this action was a
colder FW supply to the core. Between four and five minutes after the pump trip, the
operators observed power oscillations with amplitude range from 25% to 50%
of the rated value. The reactor scram occurred automatically on high neutron flux at
118% of the rated power at about seven minutes after the pumps tripped. This
accident was analyzed in many works, as in [3,
4], for example.In 1995, in Laguna Verde BWR/5, an instability event occurred during the startup
process. The analysis of neutron noise showed that the transition from stability to
instability is a gradual process that can be stopped by an earlier alarm indication
[5]. No damage to the plant was reported.
In addition, the Oskarshamn-3 BWR experienced power oscillations in February 1998. A
review of the possible causes suggested that the oscillations resulted from the
particular used control rod sequence and the power distribution obtained as a result
[6].In November 2001, an in-phase neutron flux oscillation occurred at Philippsburg-1 NPP
after an FW temperature transient. A similar event occurred at the Swedish BWR
Oskarshamn-2 in February 1999. In both events, a scram terminated the neutron flux
oscillation, but only at the fixed scram set point at 120% and 132% power
levels, respectively. In both cases, control rod insertion was activated too late to
limit the oscillations effectively before reaching the scram set points [7].After the first instability events, authorities in all countries with BWRs required a
review of the stability features of their reactors. The authorities include the
requirements of analyses in the safety analysis reports and changes in the
procedures and plant safety systems. The major safety concern associated with
instability is the cooling of the fuel and cladding integrity.Several occurred instability events in BWR plants were induced intentionally as
experiments. For example, series of turbine trips and stability tests were conducted
in the Peach Bottom-2 BWR in 1977. The low flow stability tests (LFST) were done
along the low-flow end of the rated power-flow line, and along the power-flow line
corresponding to a minimum recirculation pump speed in the region of the power-flow
map where oscillations have a higher probability [8]. The main aim of these tests was to provide a database for the
qualification of transient design methods used for reactor analyses at operating
conditions.A stability test was performed at Forsmark-1 BWR in January, 1989, during a startup
operation after the scram due to a turbine trip. In total, 36 signals including
those from APRM and LPRM, total core flow, and local channel flows were recorded on
a digital computer. These data were useful to the validation of analysis techniques
as verified, for example, in [9]. In 1990, a
stability test was conducted in the BWR Leibstadt to test the ability of the
monitoring system to cope with demanding operation situations; the power
oscillations were transformed from the in-phase mode into the out-of-phase mode by
the removal of some control rods. When the maximal oscillation amplitude was
reached, it was suppressed by decreasing the power [10].The European Commission's NACUSP project, started in December 2000, investigates
natural circulation and stability performance of BWRs. One of its main aims is
understanding the physics of the phenomena involved during the startup phase of
natural-circulation-cooled BWRs, providing a large experimental database and
validating state-of-the-art thermal-hydraulic codes in the low-pressure low-power
operational region of these reactors [11].


## 4. Methods and Tools to Study BWR Instabilities

To study instabilities in BWRs, many numerical models and computer codes have been
developed. The methods have been validated using data provided by signals of several
experimental tests and by inadvertent events. The good agreement found between
computational analyses and the available experimental data has contributed to better
understanding of the BWR instability phenomena.


### 4.1. Mathematical Models

Many literature works are related to numerical methods to study the BWR
instabilities. Several authors use the modal method to investigate the
phenomena, as in [13–15], for example. The modal analysis is a
process conceived to determine the dynamic characteristics of a system, in the
form of natural frequencies, decay factors, and oscillating modes. This
information is used to formulate a mathematical model of the dynamic behavior of
the system, denominated as the “modal model.”The decay ratio (DR) and the natural frequency (NF) of the oscillations are
typical parameters used to evaluate the instabilities. Other parameters can also
provide valuable information, such as the Lyapunov exponents associated to the
time series. In fact, Lyapunov exponents are also used as a measure of the
stability of the neutronic time series [16].Parametric or nonparametric methods can be used to evaluate the decay ratio. For
nonparametric methods, DR is evaluated from the autocorrelation function of the
signal. For parametric methods, it is evaluated from the impulse response of the
system or from its effective transfer function [17]. Different parametric models are actually used, being that the
autoregressive moving average (ARMA), the autoregressive (AR) or the
moving-average (MA) are the most common ones. For the same time series signal,
DR can have significant variation on its result depending on the method selected
for its calculation.Two different algorithms, ADRI (analysis of decay ratio instability), a
parametric method, and DRAT (decay ratio), a nonparametric method, described in
the two next sections, respectively, were used in this work to perform DR and NF
calculations.


#### 4.1.1. ADRI

ADRI code [18] is a package written
with MATLAB scripts and it can evaluate DR and natural frequency for both
signals obtained from transient initiated by short perturbation and noise
during stationary operation. It may be applied to all types of signal, with
or without noise, with high or low DR. ADRI was applied to calculate DR and
NF of Ringhals-1 NPP stability benchmark presenting results in good
agreement with the experimental DR and NF.At the base of ADRI, there are (i) the “AR” MATLAB function
that estimates the parameters of an autoregressive (AR) model, (ii) the
“IDSIM” MATLAB function simulating a dynamic system, and (iii)
the evaluation of DR and NF as the average of decay ratios and natural
frequencies, respectively, of the considered couple of peaks.


#### 4.1.2. DRAT

The method proposed here by Ambrosini is based on the form of the general
differential equation for a second-order system in free evolution, being(1)y¨+by˙+cy+d=0. The basic idea is to extract from the available transient
data estimates of y and its time derivatives to optimise the parameters b, c, and d with the minimum square technique. This amounts to solve,
by the minimum square technique, the overspecified system of equations
(2)fk(x1,x2,x3)=x1y˙k+x2yk+x3+y¨k=0,(k=2,…,N−1), in which N is the number of the available data of yk; the elements x1,x2,x3 are the coefficients of the second-order model, and also
the unknown of the problem that will permit to find the angular frequency
and the damping factor of the system. Estimates of the involved derivatives
can be found as finite difference approximations, for example, of the kind
(3)y¨k≈2tk+1−tk−1(yk+1−yktk+1−tk−yk−yk−1tk−tk−1);y˙k≈12(yk+1−yktk+1−tk+yk−yk−1tk−tk−1).Making use of the minimum square technique, the following definitions apply:(4)Q(x1,x2)=∑k=2N−1fk2(x1,x2,x3);∂Q∂xi=∑k=2N−1fk∂fk∂xi=0,(i=1,2,3);∂fk∂x1=y.k,∂fk∂x2=yk,∂fk∂x3=1, leading to the linear algebraic system (5)∑k=2N−1y˙k[x1y˙k+x2yk+x3+y¨k]=0∑k=2N−1yk[x1y˙k+x2y^k+x3+y¨k]=0∑k=2N−1[x1y˙k+x2y^k+x3+y¨k]=0 or (6)[∑k=2N−1y˙k2]x1+[∑k=2N−1y˙kyk]x2+[∑k=2N−1y˙k]x3=−[∑k=2N−1y˙ky¨k][∑k=2N−1y^ky˙k]x1+[∑k=2N−1yk2]x2+[∑k=2N−1yk]x3=−[∑k=2N−1y^ky¨k],∑k=2N−1[y˙k]x1+∑k=2N−1[yk]x2+∑k=2N−1[1]x3=−[∑k=2N−1y¨k] which represents the minimum square version of the
overspecified system introduced above.Solution of this system allows calculating the damping factor and the
frequency of oscillation of the system, interpreted as a second-order one.
In fact, from the characteristic equation of the original differential
equation(7)z2+x1z+x2=0, it is found that (8)z1,2=−x12±x124−x2, and then the general solution of the equation becomes
(9)y(t)=Ae[−x1/2+(x12/4)−x2]t+Be[−x1/2−(x12/4)−x2]t.As usual, two cases can be considered.(i)x12/4−x2>0, that is, nonoscillatory behaviour, putting (10)zR=−x12,z^R=x124−x2, the general solution is (11)y(t)=ezRt[A¯sinh(z^Rt)+B¯cosh(z^Rt)].(ii)x12/4−x2<0, that is, oscillatory behaviour, putting (12)zR=−x12,zI=|x124−x2|. From the above, it can be clearly understood that the
computed time evolution is consistent with the following theoretical, purely
second-order time evolution: (13)y(t)=ezRt[A¯sin(zIt)+B¯cos(zIt)].It is clear that the algorithm can be applied only while considering the
behavior of a system that, for small oscillations, can be considered
approximately linear after a perturbation, with no explicit forcing.


### 4.2. System Analysis Computer Codes

Computer programs developed for the modelling and the transient simulation of a
complete nuclear-power plant with a high degree of detail are called system
codes. Different choices are adopted for neutron kinetics and two-phase flow
modelling. The application of thermal hydraulics (TH) and neutron kinetics (NK)
codes to LWR analyses was discussed, for example, in the three volumes edited by
the project CRISSUE-SThe acronym CRISSUE-S project stands for Critical Issues in Nuclear
Reactor Technology, a state-of-the-art report. [10, 19, 20].
Specifically, the project CRISSUE-S treated the interactions between neutron
kinetics and thermal hydraulics that affect neutron moderation and influence the
transient performance of the NPPs.Nowadays, the nuclear industry and the scientific community turned their
attention to the development of coupled 3D neutron-kinetics and
thermal-hydraulic system codes to investigate BWR instabilities, in particular,
the regional (out-of-phase) type. The coupled system codes can model accurately
not only reactivity-initiated accidents (RIAs), but also typical reactor
operational transients as turbine trips. These programs are often called
“best-estimate” analysis tools and describe, in a more realistic
way, the local core effects, coupled reactor core, and plant dynamics
interactions.Different coupling code methodologies have been used as, for example,
TRAC-BF1/ENTREÉ [21], RELAP5-3D
[22, 23], TRAC-BF1/RAMONA [24],
MARS/MASTER [25], RETRAN-3D [26], TRAC-BF1/NEM [27], RELAP5/PANBOX/COBRA [28], and RELAP5/PARCS [29,
33].In this work, simulations of in-phase and out-of-phase instabilities in a BWR are
being presented. The thermal-hydraulic system code RELAP5 [30] and the 3D neutron-kinetic code PARCS [31] have been used in a coupled way for
performing the transient simulation. In particular, the PARCS code is used to
evaluate the 3D space-time core power history; it uses a nonlinear nodal method
to solve the two energy group neutron diffusion equations. In the calculation,
PARCS makes use of the moderator temperature and density and of the fuel
temperature calculated by RELAP5 to evaluate the appropriate feedback effects in
the neutron cross sections. Likewise, RELAP5 takes the space-dependent power
calculated in PARCS and solves the heat conduction in the core heat structures.
The coupling process between RELAP5 and PARCS codes is done through a parallel
virtual machine (PVM) environment.The temporal coupling of RELAP5 and PARCS is explicit in nature, and the two
codes are locked at the same time step. For this implementation, the RELAP5
solution lags the PARCS solution by one time step. Specifically, the advancement
of the time step begins with RELAP5 obtaining the solution to the hydrodynamic
field equations using the power from the previous time step. The property data
obtained from this solution is then sent to PARCS and the power at the current
time step is computed.The two processes are loaded in parallel and the PARCS process transfers the
nodal power data to the TH process. The TH process then sends the temperature
(fuel and coolant) and density data back to the PARCS process.The adopted calculation sequence is sketched in Figure2. The user must run two programs simultaneously; the
following sequence can be used during RELAP5/PARCS coupling: (i)
RELAP5 runs in the stand-alone mode for flow initialization (invoking
no PARCS calculations) and generates a restart file at the end of
the calculation (RELAP5 steady state stand alone);(ii)
PVM is launched; using the above restart file, the coupled
steady-state case runs and generates the steady-state restart files
for both PARCS and RELAP5 (RELAP5/PARCS coupled steady state);(iii)
using the restart files, the coupled transient case is launched
(RELAP5/PARCS coupled transient).Figure 2
Scheme of coupling among RELAP5 and PARCS codes.The original mapping between neutronic and thermal-hydraulic codes was explicit
in that the fractions of different TH nodes belonging to a neutronic node had to
be specified in the MAPTAB file for all the neutronic nodes. The postprocessor
WinGraf has been used to read the unformatted binary output data from the
RELAP5. Two different algorithms were then used to calculate the DR from the
power oscillation signals obtained from the transient calculations.


## 4.1. Mathematical Models

Many literature works are related to numerical methods to study the BWR
instabilities. Several authors use the modal method to investigate the
phenomena, as in [13–15], for example. The modal analysis is a
process conceived to determine the dynamic characteristics of a system, in the
form of natural frequencies, decay factors, and oscillating modes. This
information is used to formulate a mathematical model of the dynamic behavior of
the system, denominated as the “modal model.”The decay ratio (DR) and the natural frequency (NF) of the oscillations are
typical parameters used to evaluate the instabilities. Other parameters can also
provide valuable information, such as the Lyapunov exponents associated to the
time series. In fact, Lyapunov exponents are also used as a measure of the
stability of the neutronic time series [16].Parametric or nonparametric methods can be used to evaluate the decay ratio. For
nonparametric methods, DR is evaluated from the autocorrelation function of the
signal. For parametric methods, it is evaluated from the impulse response of the
system or from its effective transfer function [17]. Different parametric models are actually used, being that the
autoregressive moving average (ARMA), the autoregressive (AR) or the
moving-average (MA) are the most common ones. For the same time series signal,
DR can have significant variation on its result depending on the method selected
for its calculation.Two different algorithms, ADRI (analysis of decay ratio instability), a
parametric method, and DRAT (decay ratio), a nonparametric method, described in
the two next sections, respectively, were used in this work to perform DR and NF
calculations.


### 4.1.1. ADRI

ADRI code [18] is a package written
with MATLAB scripts and it can evaluate DR and natural frequency for both
signals obtained from transient initiated by short perturbation and noise
during stationary operation. It may be applied to all types of signal, with
or without noise, with high or low DR. ADRI was applied to calculate DR and
NF of Ringhals-1 NPP stability benchmark presenting results in good
agreement with the experimental DR and NF.At the base of ADRI, there are (i) the “AR” MATLAB function
that estimates the parameters of an autoregressive (AR) model, (ii) the
“IDSIM” MATLAB function simulating a dynamic system, and (iii)
the evaluation of DR and NF as the average of decay ratios and natural
frequencies, respectively, of the considered couple of peaks.


### 4.1.2. DRAT

The method proposed here by Ambrosini is based on the form of the general
differential equation for a second-order system in free evolution, being(1)y¨+by˙+cy+d=0. The basic idea is to extract from the available transient
data estimates of y and its time derivatives to optimise the parameters b, c, and d with the minimum square technique. This amounts to solve,
by the minimum square technique, the overspecified system of equations
(2)fk(x1,x2,x3)=x1y˙k+x2yk+x3+y¨k=0,(k=2,…,N−1), in which N is the number of the available data of yk; the elements x1,x2,x3 are the coefficients of the second-order model, and also
the unknown of the problem that will permit to find the angular frequency
and the damping factor of the system. Estimates of the involved derivatives
can be found as finite difference approximations, for example, of the kind
(3)y¨k≈2tk+1−tk−1(yk+1−yktk+1−tk−yk−yk−1tk−tk−1);y˙k≈12(yk+1−yktk+1−tk+yk−yk−1tk−tk−1).Making use of the minimum square technique, the following definitions apply:(4)Q(x1,x2)=∑k=2N−1fk2(x1,x2,x3);∂Q∂xi=∑k=2N−1fk∂fk∂xi=0,(i=1,2,3);∂fk∂x1=y.k,∂fk∂x2=yk,∂fk∂x3=1, leading to the linear algebraic system (5)∑k=2N−1y˙k[x1y˙k+x2yk+x3+y¨k]=0∑k=2N−1yk[x1y˙k+x2y^k+x3+y¨k]=0∑k=2N−1[x1y˙k+x2y^k+x3+y¨k]=0 or (6)[∑k=2N−1y˙k2]x1+[∑k=2N−1y˙kyk]x2+[∑k=2N−1y˙k]x3=−[∑k=2N−1y˙ky¨k][∑k=2N−1y^ky˙k]x1+[∑k=2N−1yk2]x2+[∑k=2N−1yk]x3=−[∑k=2N−1y^ky¨k],∑k=2N−1[y˙k]x1+∑k=2N−1[yk]x2+∑k=2N−1[1]x3=−[∑k=2N−1y¨k] which represents the minimum square version of the
overspecified system introduced above.Solution of this system allows calculating the damping factor and the
frequency of oscillation of the system, interpreted as a second-order one.
In fact, from the characteristic equation of the original differential
equation(7)z2+x1z+x2=0, it is found that (8)z1,2=−x12±x124−x2, and then the general solution of the equation becomes
(9)y(t)=Ae[−x1/2+(x12/4)−x2]t+Be[−x1/2−(x12/4)−x2]t.As usual, two cases can be considered.(i)x12/4−x2>0, that is, nonoscillatory behaviour, putting (10)zR=−x12,z^R=x124−x2, the general solution is (11)y(t)=ezRt[A¯sinh(z^Rt)+B¯cosh(z^Rt)].(ii)x12/4−x2<0, that is, oscillatory behaviour, putting (12)zR=−x12,zI=|x124−x2|. From the above, it can be clearly understood that the
computed time evolution is consistent with the following theoretical, purely
second-order time evolution: (13)y(t)=ezRt[A¯sin(zIt)+B¯cos(zIt)].It is clear that the algorithm can be applied only while considering the
behavior of a system that, for small oscillations, can be considered
approximately linear after a perturbation, with no explicit forcing.


## 4.1.1. ADRI

ADRI code [18] is a package written
with MATLAB scripts and it can evaluate DR and natural frequency for both
signals obtained from transient initiated by short perturbation and noise
during stationary operation. It may be applied to all types of signal, with
or without noise, with high or low DR. ADRI was applied to calculate DR and
NF of Ringhals-1 NPP stability benchmark presenting results in good
agreement with the experimental DR and NF.At the base of ADRI, there are (i) the “AR” MATLAB function
that estimates the parameters of an autoregressive (AR) model, (ii) the
“IDSIM” MATLAB function simulating a dynamic system, and (iii)
the evaluation of DR and NF as the average of decay ratios and natural
frequencies, respectively, of the considered couple of peaks.


## 4.1.2. DRAT

The method proposed here by Ambrosini is based on the form of the general
differential equation for a second-order system in free evolution, being(1)y¨+by˙+cy+d=0. The basic idea is to extract from the available transient
data estimates of y and its time derivatives to optimise the parameters b, c, and d with the minimum square technique. This amounts to solve,
by the minimum square technique, the overspecified system of equations
(2)fk(x1,x2,x3)=x1y˙k+x2yk+x3+y¨k=0,(k=2,…,N−1), in which N is the number of the available data of yk; the elements x1,x2,x3 are the coefficients of the second-order model, and also
the unknown of the problem that will permit to find the angular frequency
and the damping factor of the system. Estimates of the involved derivatives
can be found as finite difference approximations, for example, of the kind
(3)y¨k≈2tk+1−tk−1(yk+1−yktk+1−tk−yk−yk−1tk−tk−1);y˙k≈12(yk+1−yktk+1−tk+yk−yk−1tk−tk−1).Making use of the minimum square technique, the following definitions apply:(4)Q(x1,x2)=∑k=2N−1fk2(x1,x2,x3);∂Q∂xi=∑k=2N−1fk∂fk∂xi=0,(i=1,2,3);∂fk∂x1=y.k,∂fk∂x2=yk,∂fk∂x3=1, leading to the linear algebraic system (5)∑k=2N−1y˙k[x1y˙k+x2yk+x3+y¨k]=0∑k=2N−1yk[x1y˙k+x2y^k+x3+y¨k]=0∑k=2N−1[x1y˙k+x2y^k+x3+y¨k]=0 or (6)[∑k=2N−1y˙k2]x1+[∑k=2N−1y˙kyk]x2+[∑k=2N−1y˙k]x3=−[∑k=2N−1y˙ky¨k][∑k=2N−1y^ky˙k]x1+[∑k=2N−1yk2]x2+[∑k=2N−1yk]x3=−[∑k=2N−1y^ky¨k],∑k=2N−1[y˙k]x1+∑k=2N−1[yk]x2+∑k=2N−1[1]x3=−[∑k=2N−1y¨k] which represents the minimum square version of the
overspecified system introduced above.Solution of this system allows calculating the damping factor and the
frequency of oscillation of the system, interpreted as a second-order one.
In fact, from the characteristic equation of the original differential
equation(7)z2+x1z+x2=0, it is found that (8)z1,2=−x12±x124−x2, and then the general solution of the equation becomes
(9)y(t)=Ae[−x1/2+(x12/4)−x2]t+Be[−x1/2−(x12/4)−x2]t.As usual, two cases can be considered.(i)x12/4−x2>0, that is, nonoscillatory behaviour, putting (10)zR=−x12,z^R=x124−x2, the general solution is (11)y(t)=ezRt[A¯sinh(z^Rt)+B¯cosh(z^Rt)].(ii)x12/4−x2<0, that is, oscillatory behaviour, putting (12)zR=−x12,zI=|x124−x2|. From the above, it can be clearly understood that the
computed time evolution is consistent with the following theoretical, purely
second-order time evolution: (13)y(t)=ezRt[A¯sin(zIt)+B¯cos(zIt)].It is clear that the algorithm can be applied only while considering the
behavior of a system that, for small oscillations, can be considered
approximately linear after a perturbation, with no explicit forcing.


## 4.2. System Analysis Computer Codes

Computer programs developed for the modelling and the transient simulation of a
complete nuclear-power plant with a high degree of detail are called system
codes. Different choices are adopted for neutron kinetics and two-phase flow
modelling. The application of thermal hydraulics (TH) and neutron kinetics (NK)
codes to LWR analyses was discussed, for example, in the three volumes edited by
the project CRISSUE-SThe acronym CRISSUE-S project stands for Critical Issues in Nuclear
Reactor Technology, a state-of-the-art report. [10, 19, 20].
Specifically, the project CRISSUE-S treated the interactions between neutron
kinetics and thermal hydraulics that affect neutron moderation and influence the
transient performance of the NPPs.Nowadays, the nuclear industry and the scientific community turned their
attention to the development of coupled 3D neutron-kinetics and
thermal-hydraulic system codes to investigate BWR instabilities, in particular,
the regional (out-of-phase) type. The coupled system codes can model accurately
not only reactivity-initiated accidents (RIAs), but also typical reactor
operational transients as turbine trips. These programs are often called
“best-estimate” analysis tools and describe, in a more realistic
way, the local core effects, coupled reactor core, and plant dynamics
interactions.Different coupling code methodologies have been used as, for example,
TRAC-BF1/ENTREÉ [21], RELAP5-3D
[22, 23], TRAC-BF1/RAMONA [24],
MARS/MASTER [25], RETRAN-3D [26], TRAC-BF1/NEM [27], RELAP5/PANBOX/COBRA [28], and RELAP5/PARCS [29,
33].In this work, simulations of in-phase and out-of-phase instabilities in a BWR are
being presented. The thermal-hydraulic system code RELAP5 [30] and the 3D neutron-kinetic code PARCS [31] have been used in a coupled way for
performing the transient simulation. In particular, the PARCS code is used to
evaluate the 3D space-time core power history; it uses a nonlinear nodal method
to solve the two energy group neutron diffusion equations. In the calculation,
PARCS makes use of the moderator temperature and density and of the fuel
temperature calculated by RELAP5 to evaluate the appropriate feedback effects in
the neutron cross sections. Likewise, RELAP5 takes the space-dependent power
calculated in PARCS and solves the heat conduction in the core heat structures.
The coupling process between RELAP5 and PARCS codes is done through a parallel
virtual machine (PVM) environment.The temporal coupling of RELAP5 and PARCS is explicit in nature, and the two
codes are locked at the same time step. For this implementation, the RELAP5
solution lags the PARCS solution by one time step. Specifically, the advancement
of the time step begins with RELAP5 obtaining the solution to the hydrodynamic
field equations using the power from the previous time step. The property data
obtained from this solution is then sent to PARCS and the power at the current
time step is computed.The two processes are loaded in parallel and the PARCS process transfers the
nodal power data to the TH process. The TH process then sends the temperature
(fuel and coolant) and density data back to the PARCS process.The adopted calculation sequence is sketched in Figure2. The user must run two programs simultaneously; the
following sequence can be used during RELAP5/PARCS coupling: (i)
RELAP5 runs in the stand-alone mode for flow initialization (invoking
no PARCS calculations) and generates a restart file at the end of
the calculation (RELAP5 steady state stand alone);(ii)
PVM is launched; using the above restart file, the coupled
steady-state case runs and generates the steady-state restart files
for both PARCS and RELAP5 (RELAP5/PARCS coupled steady state);(iii)
using the restart files, the coupled transient case is launched
(RELAP5/PARCS coupled transient).Figure 2
Scheme of coupling among RELAP5 and PARCS codes.The original mapping between neutronic and thermal-hydraulic codes was explicit
in that the fractions of different TH nodes belonging to a neutronic node had to
be specified in the MAPTAB file for all the neutronic nodes. The postprocessor
WinGraf has been used to read the unformatted binary output data from the
RELAP5. Two different algorithms were then used to calculate the DR from the
power oscillation signals obtained from the transient calculations.


## 5. Analysed Events and Results

In this work, data from experimental low flow stability tests (LFSTs) have been
compared with results obtained with coupled 3D simulations. Other transient cases,
requiring the use of a 3D coupled analysis, have been also simulated (feedwater
temperature decrease, recirculation pump trip, and control rod banks movement) using
the same operating conditions of the LFST. These cases must be considered as
sensitivity analyses with no possibility of comparison with measured data. All the
simulated events are being described next, as well as the obtained results.


### 5.1. Thermal-Hydraulic and Neutronic Model

Peach Bottom Unit 2 is a direct-cycle BWR/4 of General Electric type that has
been subjected to stability testing. Three turbine trip tests and four series of
low flow stability tests (LFST) have been performed during the first quarter of
1977 at the end of cycle 2. The LFST were performed in the region of the
power-flow map where the highest probability of oscillations occurs.The Peach Bottom nodalisation for RELAP5 and PARCS was based on the benchmark
specification document for the turbine trip test (TT) [32], and on data in the related tests report [8]. Details of the adopted nodalization
methodology, developed by the University of Pisa, are described in [36]. The methodology was validated in
relation with the TT test [33, 36] and also for pressure perturbation
stability tests [29]. The Peach Bottom
NPP core was divided into 33 heated regions representing the 764 real core fuel
assemblies, modelled according to the RELAP5 code requirements; channels with
common characteristics were grouped together. In particular, each channel groups
a certain number of fuel assemblies; they were chosen according to their
thermal-hydraulic and kinetic properties, taking into account the lattice type,
the relative power, the inlet flow area, and the relative position within the
core.Figure3 shows the main elements of the
RELAP5 nodalization. Figure 4 represents
part of the nodalization corresponding to the reactor core; in the figure, the
identification number is related to the pipe component in the nodalization. The
core active zone was axially subdivided into 24 meshes with 15.24 cm
each. In a recent study, Ambrosini and Ferreri [39] investigated stability boundaries obtained from a RELAP5 model
for a boiling channel of 3.6 m with 48 and 24 meshes. The results showed
that the stability boundaries predicted with 48 and 24 nodes are very similar.
Therefore, the use of 24 meshes limits the complexity of the model reducing the
calculation time and conserving the accuracy of the results.Figure 3
Peach bottom BWR modelled by RELAP5.Figure 4
Detail of the plant nodalisation with the 33 TH channels in the reactor
core.To represent the reactor core neutronic behavior by the PARCS code, the reactor
core was discretized into parallelepipedal nodes, where the nuclear properties
are assumed to be constant. Radially, 18 fuel types and one reflector node were
defined, whereas, axially, the core was subdivided into 26 axial nodes; the
first and the last nodes represent the reflector zones. In total, 435
compositions or neutronic nodes were considered to represent the kinetic
behavior of the core.


### 5.2. Low Flow Stability Tests, Pressure Perturbation (PP)

The LFST were intended to measure the reactor core stability margins at the
limiting conditions used in design and safety analyses. Table1 gives the test condition for the four
stability points (PT1, PT2, PT3, and PT4), and Figure 5 shows the location of the respective test points in the
Peach Bottom power-flow map.Table 1
LFST-experimental conditions.
Tests
Power
Mass flow rate (Core Inlet)
Enthalpy (Core Inlet)
Pressure (Core Inlet)
DR
NF
%
MW
kg/s
(%)
(kJ/kg)
(MPa)
(Hz)
PT1
60.6
1995
6753.6
52.3
1184.6
7.06
0.121
0.441
PT2
51.7
1702
5657.4
43.8
1187.8
7.01
0.121
0.471
PT3
59.2
1948
5216.4
40.4
1184.6
7.10
0.344
0.437
PT4
43.5
1434
5203.8
40.3
1183.8
7.06
0.296
0.402Figure 5
Peach bottom-2 Low-Flow stability tests in the power/flow map.For all the four stability points, the steady-state simulations were firstly
performed using the RELAP5 code stand alone in order to estimate the
thermal-hydraulic operating conditions under the assumption of fixed and uniform
axial power distribution. These initial conditions are then used to perform the
coupled calculations. In the coupled steady-state calculation, results of the
axial power profile were obtained to the four cases and compared with the
available experimental curves with good agreement in all the four cases as it
can be verified in Figure6.PT1, PT2, PT3 and PT4—experimental and calculated mean axial power
profile.
(a)
(b)
(c)
(d)In the transient experiments, the magnitude of the pressure set point steps was
selected at approximately 8 psi (0.055 MPa), which gave a good
signal-to-noise ratio in the neutron flux response and did not cause operational
difficulties during the testing. Then the series of small pressure perturbation
tests conducted at each of the LFST conditions were composed of pseudorandom
binary switching of small step inputs to the pressure-regulator reference set
point.Typical reactor-core and vessel-pressure responses, and the average and local
neutron flux signals, were taken. The neutron flux to pressure transfer
functions was estimated from the data using the fast Fourier transform (FFT)
algorithm. From this transfer function, the stability margin of the core, in
terms of the decay ratio of the fundamental oscillatory mode of response, was
determined [8].Therefore, DR and NF calculated from the experimental power oscillations data are
available (Table1) and being compared,
in this work, with the data calculated from the simulations.The coupled transient calculations were performed considering boundary conditions
in which the reactor is disturbed with one pressure spike of 0.055 MPa in
the turbine. The turbine corresponds to atmpdvol
(time-dependent volume) component type in the RELAP5 input deck (element number
675 in the nodalization) that permits to impose pressure variation in time.The pressure perturbation propagates through the steam line (SL) and reaches the
core disturbing the mass flow rate and bringing the power to an oscillatory
behavior. Figures7, 8, 9 and 10 show the power behavior after the
pressure perturbation for the cases PT1, PT2, PT3, and PT4, respectively. In all
the cases, the transient begins at the time zero of calculation.Figure 7
Power evolution—PT1.Figure 8
Power evolution—PT2.Figure 9
Power evolution—PT3.Figure 10
Power evolution—PT4.


#### 5.2.1. Case PT1

As it is shown in Figure7, power
oscillates after the pressure wave perturbation disturbs the core flow. The
power oscillation reaches the maximum value of 64% (7.4% higher
than the initial condition) and after about 40 seconds, the steady-state
conditions are re-established.The core mass flow rate oscillates reaching the maximum amplitude of
0.7%. The variation in the mass flow rate is about ten times smaller
than that of power, that is, a small perturbation of mass flow rate can
cause a great variation in the power.DR and NF from the power signal have been calculated using both algorithms
ADRI and DRAT and different time intervals. It has been observed that both
algorithms give results of DR and NF that can vary largely in accordance
with that of the time window considered in the calculation [29].The oscillation in the case PT1 presents a behavior very different in
comparison with the other three cases, PT2, PT3, and PT4, where the
oscillations involve few peaks and are quickly damped. Observing, in detail,
the final oscillation in power for PT1, it is more similar to the linear
oscillatory behavior of the cases PT2, PT3, and PT4. Therefore, DR and NF
were recalculated considering only the last power peaks for the case PT1
because the initial oscillations could be interpreted as a noise signal.
Table2 gives calculation and
experimental data of DR and NF; calculated values by ADRI and DRAT
correspond to different time windows. As it can be seen, DR predicted by
ADRI and DRAT tends to values closer to the experimental one when only the
last peaks are considered (time window from 37.5 to 70.0 s) in spite
of DRAT giving values of DR slightly higher than that from ADRI.Table 2
DR and NF results in comparison with the experimental data for
PT1.
Calculation
Data from experimental
results
ADRI
DRAT
Time window from 37.5 to
70.0 s
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT1
0.121
0.441
0.216
0.355
0.323
0.377
Time window from 0.0 to
55.0 s
0.729
0.416
1.001
0.383
Time window from 36.5 to
55.0 s
0.351
0.424
0.540
0.431


#### 5.2.2. Case PT2

A fast decrease in the amplitude of power oscillations is observed and, after
approximately 30 seconds, oscillations are terminated (Figure8). In this case, the system has a very
stable behavior, the oscillations are completely damped and the power and
the other analysed parameters return to the steady state values.DR and NF, from the power signal, were calculated using both algorithms ADRI
and DRAT. Table3 presents the
numerical results. As it can be observed from the table, both algorithms
give values of DR and NF very close to each other. However, the calculation
results are not in agreement with the experimental data. In fact, the
experimental DR for the point PT2 is very small and it is exactly the same
as for the point PT1 (DBexp=0.121) though PT1 and PT2 are operating points relatively far
from each other in the power-flow map.Table 3
DR and NF results in comparison with the experimental data for the
PT2.
Data from
experimental results
Calculation (from 3.8 to
30 s)
ADRI
DRAT
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT2
0.121
0.471
0.498
0.275
0.481
0.278The tests PT1 and PT2 were investigated in two previous works [34, 35] from the point of view of the gain between the pressure and
power during the PP event. In the results, the gains for the tests PT1 and
PT2 are underestimated by the first work and overestimated for the second
one. This seems to indicate that the two operating points are somehow
critical for simulation.It is possible that the strong transient xenon-concentration change taken
place between test conditions PT1 and PT2 [8] could have “masked” the DR results for these
cases. The xenon concentration affects the stability and, in particular, the
DR value. The redistribution of the Xe concentration, following a large
scale power change, apparently may cause decrease in the DR [2].


#### 5.2.3. Case PT3

Also in this case, the process presents a fast decrease in the power
amplitude oscillation and, after about 30 seconds, oscillations are
terminated, as it can be seen in Figure9. The system presents good stability to the PP transient. After
the perturbation, power and the other parameters return to the steady state
values.DR and NF, from the power signal, were calculated using both algorithms ADRI
and DRAT considering a time window from 5.0 to 30.0 seconds. Table4 presents the results. The results
show a very good agreement between ADRI and DRAT for the DR calculations.
However, both algorithms tend to overestimate the experimental value of DR.
Calculated NF presents a value smaller than that for the experimental one,
as occurred for the case PT2.Table 4
DR and NF results in comparison with the experimental data for
PT3.
Data from
experimental results
Calculation (from 5 to
30 s)
ADRI
DRAT
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT3
0.344
0.437
0.454
0.289
0.467
0.275A sensitivity case for the test PT3 was performed, in which the core
nodalization was modified obtaining a new configuration. The number of
heated thermal-hydraulic channels in the core changed from 33, in the
original nodalisation, to 132 [29]. A
single pressure spike of 0.055 MPa was applied during 1 second in the
turbine. In this case, DR decreases slightly and NF increases with respect
to the base case. The values obtained are DR=0.438 and NF=0.312 Hz, and these results are closer to the
experimental DR and NF.


#### 5.2.4. Case PT4

The perturbation starts at the time zero. As it can be observed in Figure10, the pressure perturbation
brings the power to oscillate. In a short time (10 seconds), the
oscillations are damped and the power oscillates around the mean value with
amplitude of 0.61%. This “residual” power oscillation, ±0.6%, is very small and cannot be considered as an actual
indication of instability from the view point of the oscillation amplitude:
only amplitudes greater than 10% of the mean value are considered as an
indication of instability [37].Table5 presents calculated DR and
NF in comparison with the experimental data for PT4. The results are in
reasonable agreement with the experimental one though NF is a little smaller
in the prediction. The results given by the algorithms are in reasonable
agreement between each other for DR as well as for NF, considering the time
window from 4.2 to 30.0 seconds.Table 5
DR and NF results in comparison with the experimental data for
PT4.
Data from
experimental results
Calculation (from 4.2 to
30 s)
ADRI
DRAT
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT4
0.296
0.402
0.217
0.282
0.349
0.285The main conclusions from all the pressure-perturbation investigations are
the following.(i)
ADRI and DRAT are capable of predicting DR and NF, in most cases,
in reasonable agreement between them, in spite of the fact that
the algorithms are based on very different mathematical
assumptions.(ii)
The results obtained for DR showed that this value changes in
dependence of the time window considered in the analysis. It is,
therefore, very important to pay attention to selecting an
adequate signal time interval, representing the linear
dependencies of the system.(iii)
Performing DR calculations, using the algorithm ADRI, needs
having a sampled signal with a minimum of about 500 points. On
the other hand, DRAT is capable of evaluating DR for any number
of points sufficient to depict a reasonably complete swinging of
the considered parameter. Anyway, a higher number of points give
a more realistic result for DR.(iv)
The results of DR, from the LFST experimental data, were obtained
with a completely different calculation method with respect to
those adopted by ADRI and DRAT, and this fact can be a cause of
the discrepancies found between experimental and calculated
DRs.(v)
The calculated DR for the point PT4 is in good agreement with the
experimental one. However, the calculations showed a small
constant power oscillation that is observed after the
perturbation. This confirms that some nonlinear effects come
into play in the analyses.(vi)
For all the analyzed cases, the frequencies of the oscillations
varied between 0.3 and 0.5 Hz, which is a typical
frequency range for this kind of instability events.


### 5.3. Feedwater (FW) Temperature Decrease

A fault of FW preheaters (e.g., due to sudden depressurization in one preheater
on the heating side) and of FW pumps may cause an FW temperature decrease that
results in colder water at core inlet. This creates the potential for reducing
the volume occupied by steam in the core and a consequent increase in moderation
and fission power.The transient calculation was performed considering only the point PT3 at the
base of an operating condition. Considering that there are no experimental data
for comparison, this calculation represents only a sensitivity analysis.The event has been simulated in the Peach Bottom by the coupled RELAP5/PARCS. To
perform the transient calculation, the FW temperature value was reduced by 10
and 50 K (two separated cases) during five seconds. The feedwater
(element number 500 in the nodalization) corresponds to atmpdvol component in the RELAP5 input deck, in which it is
possible to vary the temperature as a function of time.It was observed that no significant variation in the power evolution has occurred
for the two analyzed cases. Power experienced an increase of 11% (caseΔTFW=−50 K) and, in few seconds, it returned to the steady-state
level as it can be seen in Figure 11.Figure 11
Relative power evolution considering two cases of FW temperature
decrease: 10 K and 50 K.


### 5.4. Recirculation Pump Trip (RPT)

In similarity with the previous case of FW temperature perturbation, there are no
experimental data available for comparing the results of the calculations
performed for this type of transient, which must be again regarded as a
sensitivity analysis.The stop of a recirculation pump causes a sharp decrease in the core flow, which
generates a significant negative reactivity insertion that tends to reduce power
and, consequently, the amount of steam generated.To simulate the event, the recirculation pump speed was brought to zero (in the
RELAP5 input deck) for one and five seconds, respectively, in two different
analyses. In the transient, the pump is shut down for a short time interval and
then it is switched on again. The relative power evolutions for the two cases
are shown in Figure12. One of the two
pumps is stopped at the time zero. As it can be seen, the variation in the pump
trip duration causes a small variation in the power oscillation amplitude, and
the oscillations are terminated at the same time.Figure 12
Relative power evolution considering one recirculation pump stopped for
two different time intervals: 5 seconds and 1 second.In addition, another case was considered in which both pumps were stopped, at the
same time, for one second. As it can be noted in Figure13, the amplitude of the power oscillation in this case
is higher, as it is expected to occur. The periods of oscillation, for both
cases, are practically the same; the reactor reaches stability nearly in the
same time for these two cases.Figure 13
Relative power evolution after recirculation pump trip.DR and NF calculations were performed for both cases using the algorithm ADRI.
The core power exhibits damped in-phase oscillations with a decay ratio value
less than 1.0, characterizing a stable system after the transient event. The
results of DR and NF, presented in Table6, are very similar to those found for the case of the pressure
perturbation.Table 6
DR and NF calculated by ADRI.
Calculation case
DR
NF (Hz)
(1) One pump stopped for 1 s
0.427
0.278
(2) Two pumps stopped for 1 s
0.345
0.280


### 5.5. Control Rod Bank Movement

To simulate this hypothetical transient, the control rod banks were withdrawn
beginning from the steady state positions for the case PT3. Figure14 shows the initial control rods
distribution. The position identified by “48” represents the bank
totally removed. The simulation was performed using the
MOV_BANK card in the PARCS input data. The rod banks
were continuously withdrawn, during the period from 20 to 100 seconds. Then the
calculation carried on until 200 seconds with the same rod banks configuration.
During the analysis, out-of-phase oscillations were observed. It must be
emphasized that the purpose of this study is not to simulate a realistic reactor
transient, but just to investigate the possible mode of oscillations that could
be observed in a hypothetical case in which the reactor was brought to unstable
conditions by raising its power.Figure 14
Initial control rods distribution (point PT3—steady state).Figure15 shows the power evolution
during the calculation. Power oscillations last until the end of the
calculation, but the amplitude falls drastically after 160 seconds of
calculation and the reactor, due to subsequent feedback effects, reaches less
unstable conditions. Only after about 90 seconds of calculation, strong
oscillations in the power are observed. The analyses showed that out-of-phase
oscillations are dominant after about 100 seconds in the transient
calculation.Total power evolution.
(a)
Overall transient
(b)
100 to 120 s time windowTwo TH core channels, 11 and 22, were taken as references to demonstrate the
out-of-phase phenomenon during the transient. In the TH nodalization, channels
11 and 22 are localized in two different parts of the core with respect to the
core center. Figure16 presents the
trend of mass-flow rate in the channels 11 and 22. Observing the mass-flow
evolution in a greater detail (Figure 16(b)), it is possible to check the out-of-phase behavior: when the
mass flow reaches a maximum value in a channel, in the other, it presents a
minimum value. In addition, Figure 17
represents the void fraction for the same channels in different quarters of the
core; also this parameter is found to behave out of phase in the two considered
positions.Inlet mass flow rate evolution for two selected channels located in
different quarters of the core.
(a)
Overall transient
(b)
100 to 120 s time windowVoid fraction evolution at mid height (axial level 12) in two selected
channels.
(a)
Overall transient
(b)
100 to 120 s time windowFigure18 shows the cladding temperature
evolution at the axial level 3 in six core TH channels. The variable was
observed for all 33 TH channels. Channel 11 reached the most elevated values of
temperature. The cladding temperature increases drastically in one extreme of
the fuel assembly (axial level 3) after the rod banks are removed. This
phenomenon is directly connected with the change in axial power distribution,
which is drastically affected by the rod banks withdrawn as it can be seen in
Figure 19. Since after rod withdrawal,
the coolant density is much higher at the bottom core inlet, the expected
bottom-peaked power profile is obtained. The temperature difference between the
cladding and coolant, in time, was verified and can reach for about 700°C testifying for the occurrence of dry out.Figure 18
Cladding temperature at the axial level 3 in six different core
channels.Figure 19
Evolution of the axial power profile during the removal of the control
rod banks.Figure20 illustrates five moments of
power evolution during one cycle of transient (from 121.101 to 122.921 seconds).
As it can be seen, 3D power plots provide clearer visualization of the time
evolution of the phenomenon.Relative power evolution during a period of oscillation of
1.82 seconds (from 121.101 to 122.921 seconds in the
transient calculation).
(a)
(b)
(c)
(d)
(e)The TH nodalization has been modified from 33 to 132 TH core channels. The main
aim of this change was to eliminate the existence of identical channels in
different halves of the core. In this way, the TH mapping is independent on each
half of the core permitting to verify the independent behavior of each one
during the transient. The obtained new results for the control rod banks
movement can be seen in the currently published work [38].


## 5.1. Thermal-Hydraulic and Neutronic Model

Peach Bottom Unit 2 is a direct-cycle BWR/4 of General Electric type that has
been subjected to stability testing. Three turbine trip tests and four series of
low flow stability tests (LFST) have been performed during the first quarter of
1977 at the end of cycle 2. The LFST were performed in the region of the
power-flow map where the highest probability of oscillations occurs.The Peach Bottom nodalisation for RELAP5 and PARCS was based on the benchmark
specification document for the turbine trip test (TT) [32], and on data in the related tests report [8]. Details of the adopted nodalization
methodology, developed by the University of Pisa, are described in [36]. The methodology was validated in
relation with the TT test [33, 36] and also for pressure perturbation
stability tests [29]. The Peach Bottom
NPP core was divided into 33 heated regions representing the 764 real core fuel
assemblies, modelled according to the RELAP5 code requirements; channels with
common characteristics were grouped together. In particular, each channel groups
a certain number of fuel assemblies; they were chosen according to their
thermal-hydraulic and kinetic properties, taking into account the lattice type,
the relative power, the inlet flow area, and the relative position within the
core.Figure3 shows the main elements of the
RELAP5 nodalization. Figure 4 represents
part of the nodalization corresponding to the reactor core; in the figure, the
identification number is related to the pipe component in the nodalization. The
core active zone was axially subdivided into 24 meshes with 15.24 cm
each. In a recent study, Ambrosini and Ferreri [39] investigated stability boundaries obtained from a RELAP5 model
for a boiling channel of 3.6 m with 48 and 24 meshes. The results showed
that the stability boundaries predicted with 48 and 24 nodes are very similar.
Therefore, the use of 24 meshes limits the complexity of the model reducing the
calculation time and conserving the accuracy of the results.Figure 3
Peach bottom BWR modelled by RELAP5.Figure 4
Detail of the plant nodalisation with the 33 TH channels in the reactor
core.To represent the reactor core neutronic behavior by the PARCS code, the reactor
core was discretized into parallelepipedal nodes, where the nuclear properties
are assumed to be constant. Radially, 18 fuel types and one reflector node were
defined, whereas, axially, the core was subdivided into 26 axial nodes; the
first and the last nodes represent the reflector zones. In total, 435
compositions or neutronic nodes were considered to represent the kinetic
behavior of the core.


## 5.2. Low Flow Stability Tests, Pressure Perturbation (PP)

The LFST were intended to measure the reactor core stability margins at the
limiting conditions used in design and safety analyses. Table1 gives the test condition for the four
stability points (PT1, PT2, PT3, and PT4), and Figure 5 shows the location of the respective test points in the
Peach Bottom power-flow map.Table 1
LFST-experimental conditions.
Tests
Power
Mass flow rate (Core Inlet)
Enthalpy (Core Inlet)
Pressure (Core Inlet)
DR
NF
%
MW
kg/s
(%)
(kJ/kg)
(MPa)
(Hz)
PT1
60.6
1995
6753.6
52.3
1184.6
7.06
0.121
0.441
PT2
51.7
1702
5657.4
43.8
1187.8
7.01
0.121
0.471
PT3
59.2
1948
5216.4
40.4
1184.6
7.10
0.344
0.437
PT4
43.5
1434
5203.8
40.3
1183.8
7.06
0.296
0.402Figure 5
Peach bottom-2 Low-Flow stability tests in the power/flow map.For all the four stability points, the steady-state simulations were firstly
performed using the RELAP5 code stand alone in order to estimate the
thermal-hydraulic operating conditions under the assumption of fixed and uniform
axial power distribution. These initial conditions are then used to perform the
coupled calculations. In the coupled steady-state calculation, results of the
axial power profile were obtained to the four cases and compared with the
available experimental curves with good agreement in all the four cases as it
can be verified in Figure6.PT1, PT2, PT3 and PT4—experimental and calculated mean axial power
profile.
(a)
(b)
(c)
(d)In the transient experiments, the magnitude of the pressure set point steps was
selected at approximately 8 psi (0.055 MPa), which gave a good
signal-to-noise ratio in the neutron flux response and did not cause operational
difficulties during the testing. Then the series of small pressure perturbation
tests conducted at each of the LFST conditions were composed of pseudorandom
binary switching of small step inputs to the pressure-regulator reference set
point.Typical reactor-core and vessel-pressure responses, and the average and local
neutron flux signals, were taken. The neutron flux to pressure transfer
functions was estimated from the data using the fast Fourier transform (FFT)
algorithm. From this transfer function, the stability margin of the core, in
terms of the decay ratio of the fundamental oscillatory mode of response, was
determined [8].Therefore, DR and NF calculated from the experimental power oscillations data are
available (Table1) and being compared,
in this work, with the data calculated from the simulations.The coupled transient calculations were performed considering boundary conditions
in which the reactor is disturbed with one pressure spike of 0.055 MPa in
the turbine. The turbine corresponds to atmpdvol
(time-dependent volume) component type in the RELAP5 input deck (element number
675 in the nodalization) that permits to impose pressure variation in time.The pressure perturbation propagates through the steam line (SL) and reaches the
core disturbing the mass flow rate and bringing the power to an oscillatory
behavior. Figures7, 8, 9 and 10 show the power behavior after the
pressure perturbation for the cases PT1, PT2, PT3, and PT4, respectively. In all
the cases, the transient begins at the time zero of calculation.Figure 7
Power evolution—PT1.Figure 8
Power evolution—PT2.Figure 9
Power evolution—PT3.Figure 10
Power evolution—PT4.


### 5.2.1. Case PT1

As it is shown in Figure7, power
oscillates after the pressure wave perturbation disturbs the core flow. The
power oscillation reaches the maximum value of 64% (7.4% higher
than the initial condition) and after about 40 seconds, the steady-state
conditions are re-established.The core mass flow rate oscillates reaching the maximum amplitude of
0.7%. The variation in the mass flow rate is about ten times smaller
than that of power, that is, a small perturbation of mass flow rate can
cause a great variation in the power.DR and NF from the power signal have been calculated using both algorithms
ADRI and DRAT and different time intervals. It has been observed that both
algorithms give results of DR and NF that can vary largely in accordance
with that of the time window considered in the calculation [29].The oscillation in the case PT1 presents a behavior very different in
comparison with the other three cases, PT2, PT3, and PT4, where the
oscillations involve few peaks and are quickly damped. Observing, in detail,
the final oscillation in power for PT1, it is more similar to the linear
oscillatory behavior of the cases PT2, PT3, and PT4. Therefore, DR and NF
were recalculated considering only the last power peaks for the case PT1
because the initial oscillations could be interpreted as a noise signal.
Table2 gives calculation and
experimental data of DR and NF; calculated values by ADRI and DRAT
correspond to different time windows. As it can be seen, DR predicted by
ADRI and DRAT tends to values closer to the experimental one when only the
last peaks are considered (time window from 37.5 to 70.0 s) in spite
of DRAT giving values of DR slightly higher than that from ADRI.Table 2
DR and NF results in comparison with the experimental data for
PT1.
Calculation
Data from experimental
results
ADRI
DRAT
Time window from 37.5 to
70.0 s
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT1
0.121
0.441
0.216
0.355
0.323
0.377
Time window from 0.0 to
55.0 s
0.729
0.416
1.001
0.383
Time window from 36.5 to
55.0 s
0.351
0.424
0.540
0.431


### 5.2.2. Case PT2

A fast decrease in the amplitude of power oscillations is observed and, after
approximately 30 seconds, oscillations are terminated (Figure8). In this case, the system has a very
stable behavior, the oscillations are completely damped and the power and
the other analysed parameters return to the steady state values.DR and NF, from the power signal, were calculated using both algorithms ADRI
and DRAT. Table3 presents the
numerical results. As it can be observed from the table, both algorithms
give values of DR and NF very close to each other. However, the calculation
results are not in agreement with the experimental data. In fact, the
experimental DR for the point PT2 is very small and it is exactly the same
as for the point PT1 (DBexp=0.121) though PT1 and PT2 are operating points relatively far
from each other in the power-flow map.Table 3
DR and NF results in comparison with the experimental data for the
PT2.
Data from
experimental results
Calculation (from 3.8 to
30 s)
ADRI
DRAT
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT2
0.121
0.471
0.498
0.275
0.481
0.278The tests PT1 and PT2 were investigated in two previous works [34, 35] from the point of view of the gain between the pressure and
power during the PP event. In the results, the gains for the tests PT1 and
PT2 are underestimated by the first work and overestimated for the second
one. This seems to indicate that the two operating points are somehow
critical for simulation.It is possible that the strong transient xenon-concentration change taken
place between test conditions PT1 and PT2 [8] could have “masked” the DR results for these
cases. The xenon concentration affects the stability and, in particular, the
DR value. The redistribution of the Xe concentration, following a large
scale power change, apparently may cause decrease in the DR [2].


### 5.2.3. Case PT3

Also in this case, the process presents a fast decrease in the power
amplitude oscillation and, after about 30 seconds, oscillations are
terminated, as it can be seen in Figure9. The system presents good stability to the PP transient. After
the perturbation, power and the other parameters return to the steady state
values.DR and NF, from the power signal, were calculated using both algorithms ADRI
and DRAT considering a time window from 5.0 to 30.0 seconds. Table4 presents the results. The results
show a very good agreement between ADRI and DRAT for the DR calculations.
However, both algorithms tend to overestimate the experimental value of DR.
Calculated NF presents a value smaller than that for the experimental one,
as occurred for the case PT2.Table 4
DR and NF results in comparison with the experimental data for
PT3.
Data from
experimental results
Calculation (from 5 to
30 s)
ADRI
DRAT
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT3
0.344
0.437
0.454
0.289
0.467
0.275A sensitivity case for the test PT3 was performed, in which the core
nodalization was modified obtaining a new configuration. The number of
heated thermal-hydraulic channels in the core changed from 33, in the
original nodalisation, to 132 [29]. A
single pressure spike of 0.055 MPa was applied during 1 second in the
turbine. In this case, DR decreases slightly and NF increases with respect
to the base case. The values obtained are DR=0.438 and NF=0.312 Hz, and these results are closer to the
experimental DR and NF.


### 5.2.4. Case PT4

The perturbation starts at the time zero. As it can be observed in Figure10, the pressure perturbation
brings the power to oscillate. In a short time (10 seconds), the
oscillations are damped and the power oscillates around the mean value with
amplitude of 0.61%. This “residual” power oscillation, ±0.6%, is very small and cannot be considered as an actual
indication of instability from the view point of the oscillation amplitude:
only amplitudes greater than 10% of the mean value are considered as an
indication of instability [37].Table5 presents calculated DR and
NF in comparison with the experimental data for PT4. The results are in
reasonable agreement with the experimental one though NF is a little smaller
in the prediction. The results given by the algorithms are in reasonable
agreement between each other for DR as well as for NF, considering the time
window from 4.2 to 30.0 seconds.Table 5
DR and NF results in comparison with the experimental data for
PT4.
Data from
experimental results
Calculation (from 4.2 to
30 s)
ADRI
DRAT
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT4
0.296
0.402
0.217
0.282
0.349
0.285The main conclusions from all the pressure-perturbation investigations are
the following.(i)
ADRI and DRAT are capable of predicting DR and NF, in most cases,
in reasonable agreement between them, in spite of the fact that
the algorithms are based on very different mathematical
assumptions.(ii)
The results obtained for DR showed that this value changes in
dependence of the time window considered in the analysis. It is,
therefore, very important to pay attention to selecting an
adequate signal time interval, representing the linear
dependencies of the system.(iii)
Performing DR calculations, using the algorithm ADRI, needs
having a sampled signal with a minimum of about 500 points. On
the other hand, DRAT is capable of evaluating DR for any number
of points sufficient to depict a reasonably complete swinging of
the considered parameter. Anyway, a higher number of points give
a more realistic result for DR.(iv)
The results of DR, from the LFST experimental data, were obtained
with a completely different calculation method with respect to
those adopted by ADRI and DRAT, and this fact can be a cause of
the discrepancies found between experimental and calculated
DRs.(v)
The calculated DR for the point PT4 is in good agreement with the
experimental one. However, the calculations showed a small
constant power oscillation that is observed after the
perturbation. This confirms that some nonlinear effects come
into play in the analyses.(vi)
For all the analyzed cases, the frequencies of the oscillations
varied between 0.3 and 0.5 Hz, which is a typical
frequency range for this kind of instability events.


## 5.2.1. Case PT1

As it is shown in Figure7, power
oscillates after the pressure wave perturbation disturbs the core flow. The
power oscillation reaches the maximum value of 64% (7.4% higher
than the initial condition) and after about 40 seconds, the steady-state
conditions are re-established.The core mass flow rate oscillates reaching the maximum amplitude of
0.7%. The variation in the mass flow rate is about ten times smaller
than that of power, that is, a small perturbation of mass flow rate can
cause a great variation in the power.DR and NF from the power signal have been calculated using both algorithms
ADRI and DRAT and different time intervals. It has been observed that both
algorithms give results of DR and NF that can vary largely in accordance
with that of the time window considered in the calculation [29].The oscillation in the case PT1 presents a behavior very different in
comparison with the other three cases, PT2, PT3, and PT4, where the
oscillations involve few peaks and are quickly damped. Observing, in detail,
the final oscillation in power for PT1, it is more similar to the linear
oscillatory behavior of the cases PT2, PT3, and PT4. Therefore, DR and NF
were recalculated considering only the last power peaks for the case PT1
because the initial oscillations could be interpreted as a noise signal.
Table2 gives calculation and
experimental data of DR and NF; calculated values by ADRI and DRAT
correspond to different time windows. As it can be seen, DR predicted by
ADRI and DRAT tends to values closer to the experimental one when only the
last peaks are considered (time window from 37.5 to 70.0 s) in spite
of DRAT giving values of DR slightly higher than that from ADRI.Table 2
DR and NF results in comparison with the experimental data for
PT1.
Calculation
Data from experimental
results
ADRI
DRAT
Time window from 37.5 to
70.0 s
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT1
0.121
0.441
0.216
0.355
0.323
0.377
Time window from 0.0 to
55.0 s
0.729
0.416
1.001
0.383
Time window from 36.5 to
55.0 s
0.351
0.424
0.540
0.431


## 5.2.2. Case PT2

A fast decrease in the amplitude of power oscillations is observed and, after
approximately 30 seconds, oscillations are terminated (Figure8). In this case, the system has a very
stable behavior, the oscillations are completely damped and the power and
the other analysed parameters return to the steady state values.DR and NF, from the power signal, were calculated using both algorithms ADRI
and DRAT. Table3 presents the
numerical results. As it can be observed from the table, both algorithms
give values of DR and NF very close to each other. However, the calculation
results are not in agreement with the experimental data. In fact, the
experimental DR for the point PT2 is very small and it is exactly the same
as for the point PT1 (DBexp=0.121) though PT1 and PT2 are operating points relatively far
from each other in the power-flow map.Table 3
DR and NF results in comparison with the experimental data for the
PT2.
Data from
experimental results
Calculation (from 3.8 to
30 s)
ADRI
DRAT
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT2
0.121
0.471
0.498
0.275
0.481
0.278The tests PT1 and PT2 were investigated in two previous works [34, 35] from the point of view of the gain between the pressure and
power during the PP event. In the results, the gains for the tests PT1 and
PT2 are underestimated by the first work and overestimated for the second
one. This seems to indicate that the two operating points are somehow
critical for simulation.It is possible that the strong transient xenon-concentration change taken
place between test conditions PT1 and PT2 [8] could have “masked” the DR results for these
cases. The xenon concentration affects the stability and, in particular, the
DR value. The redistribution of the Xe concentration, following a large
scale power change, apparently may cause decrease in the DR [2].


## 5.2.3. Case PT3

Also in this case, the process presents a fast decrease in the power
amplitude oscillation and, after about 30 seconds, oscillations are
terminated, as it can be seen in Figure9. The system presents good stability to the PP transient. After
the perturbation, power and the other parameters return to the steady state
values.DR and NF, from the power signal, were calculated using both algorithms ADRI
and DRAT considering a time window from 5.0 to 30.0 seconds. Table4 presents the results. The results
show a very good agreement between ADRI and DRAT for the DR calculations.
However, both algorithms tend to overestimate the experimental value of DR.
Calculated NF presents a value smaller than that for the experimental one,
as occurred for the case PT2.Table 4
DR and NF results in comparison with the experimental data for
PT3.
Data from
experimental results
Calculation (from 5 to
30 s)
ADRI
DRAT
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT3
0.344
0.437
0.454
0.289
0.467
0.275A sensitivity case for the test PT3 was performed, in which the core
nodalization was modified obtaining a new configuration. The number of
heated thermal-hydraulic channels in the core changed from 33, in the
original nodalisation, to 132 [29]. A
single pressure spike of 0.055 MPa was applied during 1 second in the
turbine. In this case, DR decreases slightly and NF increases with respect
to the base case. The values obtained are DR=0.438 and NF=0.312 Hz, and these results are closer to the
experimental DR and NF.


## 5.2.4. Case PT4

The perturbation starts at the time zero. As it can be observed in Figure10, the pressure perturbation
brings the power to oscillate. In a short time (10 seconds), the
oscillations are damped and the power oscillates around the mean value with
amplitude of 0.61%. This “residual” power oscillation, ±0.6%, is very small and cannot be considered as an actual
indication of instability from the view point of the oscillation amplitude:
only amplitudes greater than 10% of the mean value are considered as an
indication of instability [37].Table5 presents calculated DR and
NF in comparison with the experimental data for PT4. The results are in
reasonable agreement with the experimental one though NF is a little smaller
in the prediction. The results given by the algorithms are in reasonable
agreement between each other for DR as well as for NF, considering the time
window from 4.2 to 30.0 seconds.Table 5
DR and NF results in comparison with the experimental data for
PT4.
Data from
experimental results
Calculation (from 4.2 to
30 s)
ADRI
DRAT
LFST
DR
NF (Hz)
DR
NF (Hz)
DR
NF (Hz)
PT4
0.296
0.402
0.217
0.282
0.349
0.285The main conclusions from all the pressure-perturbation investigations are
the following.(i)
ADRI and DRAT are capable of predicting DR and NF, in most cases,
in reasonable agreement between them, in spite of the fact that
the algorithms are based on very different mathematical
assumptions.(ii)
The results obtained for DR showed that this value changes in
dependence of the time window considered in the analysis. It is,
therefore, very important to pay attention to selecting an
adequate signal time interval, representing the linear
dependencies of the system.(iii)
Performing DR calculations, using the algorithm ADRI, needs
having a sampled signal with a minimum of about 500 points. On
the other hand, DRAT is capable of evaluating DR for any number
of points sufficient to depict a reasonably complete swinging of
the considered parameter. Anyway, a higher number of points give
a more realistic result for DR.(iv)
The results of DR, from the LFST experimental data, were obtained
with a completely different calculation method with respect to
those adopted by ADRI and DRAT, and this fact can be a cause of
the discrepancies found between experimental and calculated
DRs.(v)
The calculated DR for the point PT4 is in good agreement with the
experimental one. However, the calculations showed a small
constant power oscillation that is observed after the
perturbation. This confirms that some nonlinear effects come
into play in the analyses.(vi)
For all the analyzed cases, the frequencies of the oscillations
varied between 0.3 and 0.5 Hz, which is a typical
frequency range for this kind of instability events.


## 5.3. Feedwater (FW) Temperature Decrease

A fault of FW preheaters (e.g., due to sudden depressurization in one preheater
on the heating side) and of FW pumps may cause an FW temperature decrease that
results in colder water at core inlet. This creates the potential for reducing
the volume occupied by steam in the core and a consequent increase in moderation
and fission power.The transient calculation was performed considering only the point PT3 at the
base of an operating condition. Considering that there are no experimental data
for comparison, this calculation represents only a sensitivity analysis.The event has been simulated in the Peach Bottom by the coupled RELAP5/PARCS. To
perform the transient calculation, the FW temperature value was reduced by 10
and 50 K (two separated cases) during five seconds. The feedwater
(element number 500 in the nodalization) corresponds to atmpdvol component in the RELAP5 input deck, in which it is
possible to vary the temperature as a function of time.It was observed that no significant variation in the power evolution has occurred
for the two analyzed cases. Power experienced an increase of 11% (caseΔTFW=−50 K) and, in few seconds, it returned to the steady-state
level as it can be seen in Figure 11.Figure 11
Relative power evolution considering two cases of FW temperature
decrease: 10 K and 50 K.


## 5.4. Recirculation Pump Trip (RPT)

In similarity with the previous case of FW temperature perturbation, there are no
experimental data available for comparing the results of the calculations
performed for this type of transient, which must be again regarded as a
sensitivity analysis.The stop of a recirculation pump causes a sharp decrease in the core flow, which
generates a significant negative reactivity insertion that tends to reduce power
and, consequently, the amount of steam generated.To simulate the event, the recirculation pump speed was brought to zero (in the
RELAP5 input deck) for one and five seconds, respectively, in two different
analyses. In the transient, the pump is shut down for a short time interval and
then it is switched on again. The relative power evolutions for the two cases
are shown in Figure12. One of the two
pumps is stopped at the time zero. As it can be seen, the variation in the pump
trip duration causes a small variation in the power oscillation amplitude, and
the oscillations are terminated at the same time.Figure 12
Relative power evolution considering one recirculation pump stopped for
two different time intervals: 5 seconds and 1 second.In addition, another case was considered in which both pumps were stopped, at the
same time, for one second. As it can be noted in Figure13, the amplitude of the power oscillation in this case
is higher, as it is expected to occur. The periods of oscillation, for both
cases, are practically the same; the reactor reaches stability nearly in the
same time for these two cases.Figure 13
Relative power evolution after recirculation pump trip.DR and NF calculations were performed for both cases using the algorithm ADRI.
The core power exhibits damped in-phase oscillations with a decay ratio value
less than 1.0, characterizing a stable system after the transient event. The
results of DR and NF, presented in Table6, are very similar to those found for the case of the pressure
perturbation.Table 6
DR and NF calculated by ADRI.
Calculation case
DR
NF (Hz)
(1) One pump stopped for 1 s
0.427
0.278
(2) Two pumps stopped for 1 s
0.345
0.280


## 5.5. Control Rod Bank Movement

To simulate this hypothetical transient, the control rod banks were withdrawn
beginning from the steady state positions for the case PT3. Figure14 shows the initial control rods
distribution. The position identified by “48” represents the bank
totally removed. The simulation was performed using the
MOV_BANK card in the PARCS input data. The rod banks
were continuously withdrawn, during the period from 20 to 100 seconds. Then the
calculation carried on until 200 seconds with the same rod banks configuration.
During the analysis, out-of-phase oscillations were observed. It must be
emphasized that the purpose of this study is not to simulate a realistic reactor
transient, but just to investigate the possible mode of oscillations that could
be observed in a hypothetical case in which the reactor was brought to unstable
conditions by raising its power.Figure 14
Initial control rods distribution (point PT3—steady state).Figure15 shows the power evolution
during the calculation. Power oscillations last until the end of the
calculation, but the amplitude falls drastically after 160 seconds of
calculation and the reactor, due to subsequent feedback effects, reaches less
unstable conditions. Only after about 90 seconds of calculation, strong
oscillations in the power are observed. The analyses showed that out-of-phase
oscillations are dominant after about 100 seconds in the transient
calculation.Total power evolution.
(a)
Overall transient
(b)
100 to 120 s time windowTwo TH core channels, 11 and 22, were taken as references to demonstrate the
out-of-phase phenomenon during the transient. In the TH nodalization, channels
11 and 22 are localized in two different parts of the core with respect to the
core center. Figure16 presents the
trend of mass-flow rate in the channels 11 and 22. Observing the mass-flow
evolution in a greater detail (Figure 16(b)), it is possible to check the out-of-phase behavior: when the
mass flow reaches a maximum value in a channel, in the other, it presents a
minimum value. In addition, Figure 17
represents the void fraction for the same channels in different quarters of the
core; also this parameter is found to behave out of phase in the two considered
positions.Inlet mass flow rate evolution for two selected channels located in
different quarters of the core.
(a)
Overall transient
(b)
100 to 120 s time windowVoid fraction evolution at mid height (axial level 12) in two selected
channels.
(a)
Overall transient
(b)
100 to 120 s time windowFigure18 shows the cladding temperature
evolution at the axial level 3 in six core TH channels. The variable was
observed for all 33 TH channels. Channel 11 reached the most elevated values of
temperature. The cladding temperature increases drastically in one extreme of
the fuel assembly (axial level 3) after the rod banks are removed. This
phenomenon is directly connected with the change in axial power distribution,
which is drastically affected by the rod banks withdrawn as it can be seen in
Figure 19. Since after rod withdrawal,
the coolant density is much higher at the bottom core inlet, the expected
bottom-peaked power profile is obtained. The temperature difference between the
cladding and coolant, in time, was verified and can reach for about 700°C testifying for the occurrence of dry out.Figure 18
Cladding temperature at the axial level 3 in six different core
channels.Figure 19
Evolution of the axial power profile during the removal of the control
rod banks.Figure20 illustrates five moments of
power evolution during one cycle of transient (from 121.101 to 122.921 seconds).
As it can be seen, 3D power plots provide clearer visualization of the time
evolution of the phenomenon.Relative power evolution during a period of oscillation of
1.82 seconds (from 121.101 to 122.921 seconds in the
transient calculation).
(a)
(b)
(c)
(d)
(e)The TH nodalization has been modified from 33 to 132 TH core channels. The main
aim of this change was to eliminate the existence of identical channels in
different halves of the core. In this way, the TH mapping is independent on each
half of the core permitting to verify the independent behavior of each one
during the transient. The obtained new results for the control rod banks
movement can be seen in the currently published work [38].


## 6. Conclusions

RELAP5/MOD3.3 thermal-hydraulic system code and the PARCS-2.4 3D neutron kinetic code
were coupled to simulate transients in a BWR. DR and NF calculated from the
experimental power oscillations data from the LFST were compared, in this work, with
the data calculated from the simulations.The coupled system was firstly validated for the test points PT1, PT2, PT3, and PT4
in steady-state conditions. Thereafter, simulations of PP events to the four cases
were shown and compared with the available experimental data. Two algorithms, ADRI
and DRAT, were used to perform the calculations. ADRI and DRAT predicted DR and NF
values in most of the cases in reasonable agreement between them in spite of the
fact that the algorithms are based on very different mathematical assumptions.Calculated DR is in reasonable agreement with experimental results for the tests PT3
and PT4, but it is overestimated for the tests PT1 and PT2. One possible cause for
the low value of DR found in the experiments for PT1 and PT2 is the change of the
xenon concentration taking place between those first test conditions.Other interesting transient cases, considered as sensitivity analyses (feedwater
temperature decrease, recirculation pump trip and control rod banks movement), have
been simulated using the same operating conditions of the LFST PT3. Their interest
is, anyway, related to the possibility of quantifying the margin to unstable
behavior bypassing the uncertainty introduced by different definitions adopted for
DR.In the simulations, the feedwater temperature decrease did not represent a
significant variation in the power evolution and the reactor seems to be very stable
in the analyzed case of FW temperature variation. In the case of the sudden pump
trip event, the core power exhibited damped in-phase oscillations with a decay-ratio
value less than 1.0, characterizing a stable system after the transient event. The
obtained results of DR and NF were close to those found in the case of the pressure
perturbation.The control rod bank movement transient was considered to study the out-of-phase
behavior occurring in the reactor as a consequence of raising power by the removal
of the control rod banks from the BWR core. The analyses taking into account two TH
channels, localized in two different quarters of the core, demonstrated that the
mass-flow rate and the void fraction are totally out of phase after 100 second of
calculation. Obviously, in the calculation, the scram intervention was not
considered because the main interest was to assess the core parameters evolution
during an out-of-phase oscillation. The simulation showed that, if there is not a
scram intervention, the cladding temperature can reach high values and dry-out can
occur.In summary, in-phase and out-of-phase modes of power oscillation in a BWR were
presented in this work. Though uncertainties still remain in relation to the
definition of DR to be adopted for comparison with experimental data and about the
effect of thermal-hydraulic radial discretization of the core on the obtained
results, the information obtained in this work will contribute for further more
detailed studies of such complex phenomena having considerable interest in the
safety of boiling-water nuclear reactors.



---

*Source: 423175-2008-01-03.xml*