Antonella Lombardi Costa,¹ Alessandro Petruzzi,² Francesco D'Auria,² and Walter Ambrosini²

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.

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.

Figure 1 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 for UO₂). 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) 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) in which $N$ is the number of the available data of $y_k$; the elements $x_1,x_2,x_3$ 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)

Making use of the minimum square technique, the following definitions apply: $\begin{array}{cc}& Q\left(x_1,x_2\right)={\displaystyle \sum _{k=2}^{N-1}{f}_k^2\left(x_1,x_2,x_3\right)};\\ & \frac{\partial Q}{\partial x_i}={\displaystyle \sum _{k=2}^{N-1}{f}_k\frac{\partial f_k}{\partial x_i}=0,\left(i=1,2,3\right)};\\ & \frac{\partial f_k}{\partial x_1}={\dot{y}}_k,\frac{\partial f_k}{\partial x_2}=y_k,\frac{\partial f_k}{\partial x_3}=1,\end{array}$(4) leading to the linear algebraic system (5) or (6) 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) it is found that (8) and then the general solution of the equation becomes (9)

As usual, two cases can be considered.

(i) $x_1^2/4-x_2>0$, that is, nonoscillatory behaviour, putting (10) the general solution is (11)

(ii) $x_1^2/4-x_2<0$, that is, oscillatory behaviour, putting (12) 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)

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-S¹ [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 Figure 2. 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.

Figure 3 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. Table 1 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.

Figure 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 Figure 6.

Figure 6: PT1, PT2, PT3 and PT4—experimental and calculated mean axial power profile.

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 (Table 1) 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 a tmpdvol (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. Figures 7, 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 Figure 7, 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. Table 2 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.

5.2.2. Case PT2

A fast decrease in the amplitude of power oscillations is observed and, after approximately 30 seconds, oscillations are terminated (Figure 8). 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. Table 3 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 (${\text{DB}}_{\text{exp}}=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.

The 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 Figure 9. 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. Table 4 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.

A 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 $\text{DR}=0.438$ and $\text{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 Figure 10, 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, $\pm 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].

Table 5 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.

The 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 a tmpdvol 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 $\mathrm{\Delta }{\text{T}}_{\text{FW}}=-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 Figure 12. 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 Figure 13, 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 Table 6, are very similar to those found for the case of the pressure perturbation.

Table 6: DR and NF calculated by ADRI.

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. Figure 14 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).

Figure 15 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.

Figure 15: Total power evolution.

Two 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. Figure 16 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.

Figure 16: Inlet mass flow rate evolution for two selected channels located in different quarters of the core.

Figure 17: Void fraction evolution at mid height (axial level 12) in two selected channels.

Figure 18 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.

Figure 20 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.

Figure 20: Relative power evolution during a period of oscillation of 1.82 seconds (from 121.101 to 122.921 seconds in the transient calculation).

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.

Notes

¹The acronym CRISSUE-S project stands for Critical Issues in Nuclear Reactor Technology, a state-of-the-art report.

Acknowledgments

The authors acknowledge CAPES (Brazil) and the University of Pisa (Italy) for the support. This work is carried out at the University of Pisa, Italy.

References

1. J. Dorning, “Models and stability analysis of boiling water reactors,” Tech. Rep. FG07-98ID13650, University of Virginia, Charlottesville, Va, USA, 2002.
2. F. D'Auria, W. Ambrosini, and T. Anegawa, et al., “State of the art report on boiling water reactor stability (SOAR ON BWRS),” 1997, OECD-CSNI Report OECD/GD (97) 13, Paris, France.
3. W. Ambrosini, F. D'Auria, and A. Giglioli, “Coupled thermal-hydraulic and neutronic instabilities in the LaSalle-2 BWR plant,” in Proceedings of the XIII Congresso Nazionale sulla Trasmissione del Calore, Bologna, Italy, June 1995.
4. F. D'Auria, V. Pellicoro, and O. Feldmann, “Use of RELAP5 code to evaluate the BOP influence following instability events in BWRs,” in Proceedings of the 4th ASME/JSME International Conference on Nuclear Engineering (ICONE '96), vol. 3, pp. 71–77, New Orleans, La, USA, March 1996.
5. J. Blazquez and J. Ruiz, “The Laguna verde BWR/5 instability event. Some lessons learnt,” Progress in Nuclear Energy, vol. 43, no. 1–4, pp. 195–200, 2003.
6. M. Kruners, “Analysis of instability event in Oskarshamn-3,” 1998, SKI Report, 98:42, ISSN 1104-1374.
7. M. Maqua, K. Kotthoff, and W. Pointner, “Neutron flux oscillations at German BWRs,” 2002-2003, Gesellschaft für Anlagen-und Reaktorsicherheit (GRS) mbH, Annual Report.
8. L. A. Carmichael and R. O. Niemi, “Transient and stability tests at peach bottom atomic power station unit 2 at end of cycle 2,” 1978, EPRI Report NP-564.
9. R. Oguma, “Investigation of power oscillation mechanisms based on noise analysis at Forsmark-1 BWR,” Annals of Nuclear Energy, vol. 23, no. 6, pp. 469–485, 1996.
10. OECD, “Neutronics/thermal-hydraulics coupling in LWR technology,” 2004, CRISSUE-S—WP2: state-of-the-art report, vol. 2, NEA 5436, ISBN 92-64-02084-5.
11. C. Aguirre, D. Caruge, and F. Castrillo, et al., “Natural circulation and stability performance of BWRs (NACUSP),” Nuclear Engineering and Design, vol. 235, no. 2–4, pp. 401–409, 2005.
12. J. Duderstadt and L. J. Hamilton, Nuclear Reactor Analysis, John Wiley & Sons, New York, NY, USA, 1976.
13. G. Verdú, D. Ginestar, V. Vidal, and R. Miró, “Modal decomposition method for BWR stability analysis,” Journal of Nuclear Science and Technology, vol. 35, no. 8, pp. 538–546, 1998.
14. R. Miró, D. Ginestar, G. Verdú, and D. Hennig, “A nodal modal method for the neutron diffusion equation. Application to BWR instabilities analysis,” Annals of Nuclear Energy, vol. 29, no. 10, pp. 1171–1194, 2002.
15. D. Ginestar, R. Miró, G. Verdú, and D. Hennig, “A transient modal analysis of a BWR instability event,” Journal of Nuclear Science and Technology, vol. 39, no. 5, pp. 554–563, 2002.
16. C. Pereira, G. Verdú, J. L. Muñoz-Cobo, and R. Sanchis, “BWR stability from dynamic reconstruction and autoregressive model analysis: application to Cofrentes Nuclear Power Plant,” Progress in Nuclear Energy, vol. 27, no. 1, pp. 51–68, 1992.
17. A. Dokhane, H. Ferroukhi, M. A. Zimmermann, and C. Aguirre, “Spatial and model-order based reactor signal analysis methodology for BWR core stability evaluation,” Annals of Nuclear Energy, vol. 33, no. 16, pp. 1329–1338, 2006.
18. A. Petruzzi and G. Galassi, “Developing of a software for valuation of DR and NF,” DIMNP Report NT 572(03), University of Pisa, Pisa, Italy, 2003.
19. OECD, “Neutronics/thermal-hydraulics coupling in LWR technology,” 2004, CRISSUE-S—WP3: achievements and recommendations report, vol. 3, NEA 5434, ISBN 92-64-02085-3.
20. OECD, “Neutronics/thermal-hydraulics coupling in LWR technology,” 2004, CRISSUE-S—WP1: data requirements and databases needed for transient simulations and qualification, vol. 1, NEA 4452-ISBN 92-64-02083-7.
21. A. Hotta, T. Anegawa, T. Hara, and H. Ninokata, “Simulation of BWR one-pump trip transient by TRAC-BF1/ENTRÉE,” Nuclear Technology, vol. 142, no. 3, pp. 205–229, 2003.
22. E. Uspuras, A. Kaliatka, and E. Bubelis, “Validation of coupled neutronic/thermal-hydraulic code RELAP5-3D for RBMK-1500 reactor analysis application,” Annals of Nuclear Energy, vol. 31, no. 15, pp. 1667–1708, 2004.
23. A. Lombardi Costa, A. Petruzzi, and F. D'Auria, “RELAP5-3D analysis of pressure perturbation at the peach bottom BWR during low-flow stability tests,” in Proceedings of the 14th International Conference on Nuclear Engineering (ICONE '06), Miami, Fla, USA, July 2006.
24. A. P. López and J. B. M. Sandoval, “A methodology for the coupling of RAMONA-3B neutron kinetics and TRAC-BF1 thermal-hydraulics,” Annals of Nuclear Energy, vol. 32, no. 6, pp. 621–634, 2005.
25. J.-J. Jeong, W. J. Lee, and B. D. Chung, “Simulation of a main steam line break accident using a coupled “system thermal-hydraulics, three-dimensional reactor kinetics, and hot channel analysis” code,” Annals of Nuclear Energy, vol. 33, no. 9, pp. 820–828, 2006.
26. W. Barten, P. Coddington, and H. Ferroukhi, “RETRAN-3D analysis of the base case and the four extreme cases of the OECD/NRC peach bottom 2 turbine trip benchmark,” Annals of Nuclear Energy, vol. 33, no. 2, pp. 99–118, 2006.
27. A. Hotta, M. Honma, H. Ninokata, and Y. Matsui, “BWR regional instability analysis by TRAC/BF1-ENTRÉE-I: application to density-wave oscillation,” Nuclear Technology, vol. 135, no. 1, pp. 1–16, 2001.
28. D. G. Cacuci, “Dimensionally adaptive dynamic switching and adjoint sensitivity analysis: new features of the RELAP5/PANBOX/COBRA code system for reactor safety transients,” Nuclear Engineering and Design, vol. 202, no. 2-3, pp. 325–338, 2000.
29. A. Lombardi Costa, BWR instability analysis by coupled 3D neutron-kinetic and thermal-hydraulic codes, Ph.D. thesis, Dipartimento di Ingegneria Meccanica, Nucleare e della Produzione, University of Pisa, Pisa, Italy, 2007, http://etd.adm.unipi.it/theses/available/etd-05132007-145646/.
30. US NRC, “RELAP5/MOD3.3 code manuals,” 2001, Idaho National Engineering Laboratory, NUREG/CR-5535.
31. H. G. Joo, D. Barber, G. Jiang, and T. J. Downar, “PARCS: a multi-dimensional two-group reactor kinetics code based on the non-linear analytic nodal method,” 1998, PU/NE-98-26, Purdue University.
32. J. Solis, K. Ivanov, B. Sarikaya, A. Olson, and K. W. Hunt, “Boiling water reactor turbine trip (TT) benchmark,” 2001, volume 1: final specifications, NEA/NSC/DOC.
33. A. Bousbia-Salah, J. Vedovi, F. D'Auria, K. Ivanov, and G. Galassi, “Analysis of the peach bottom turbine trip 2 experiment by coupled RELAP5-PARCS three-dimensional codes,” Nuclear Science and Engineering, vol. 148, no. 2, pp. 337–353, 2004.
34. U. S. Rohatgi, A. N. Mallen, H. S. Cheng, and W. Wulff, “Validation of the engineering plant analyzer methodology with peach bottom 2 stability tests,” Nuclear Engineering and Design, vol. 151, no. 1, pp. 145–156, 1994.
35. R. Taleyarkhan, R. T. Lahey Jr., A. F. McFarlane, and M. Z. Podowski, “Benchmarking and qualification of the NUFREQ-NPW code for best estimate prediction of multi-channel core stability margins,” in Proceedings of the Joint meeting of the European Nuclear Society and the American Nuclear Society, Washington, DC, USA, October-November 1988.
36. A. Bousbia-Salah, Overview of coupled system thermal-hidraulic 3D neutron kinetic code applications, Ph.D. thesis, Dipartimento di Ingegneria Meccanica, Nucleare e della Produzione, University of Pisa, Pisa, Italy, 2004.
37. P. K. Vijayan and A. K. Nayak, “Introduction to instabilities in natural circulation systems,” 2005, IAEA-TECDOC-1474—Natural Circulation in Water Cooled Nuclear Power Plants.
38. A. Lombardi Costa, C. Pereira, W. Ambrosini, and F. D'Auria, “Simulation of an hypothetical out-of-phase instability case in boiling water reactor by RELAP5/PARCS coupled codes,” Annals of Nuclear Energy, 2007.
39. W. Ambrosini and J. C. Ferreri, “Analysis of basic phenomena in boiling channel instabilities with different flow models and numerical schemes,” in Proceedings of the 14th International Conference on Nuclear Engineering (ICONE '06), Miami, Fla, USA, July 2006.