Abstract
While more than a decade ago reactor and thermal hydraulic calculations were tedious and often needed a lot of approximations and simplifications that forced the designers to take a very conservative approach, computational resources available nowadays allow engineers to cope with increasingly complex problems in a reasonable time. The use of best-estimate calculations provides tools to justify convenient engineering margins, reduces costs, and maximises economic benefits. In this direction, a suite of coupled best-estimate specific calculation codes was developed to analyse the behaviour of the Atucha II nuclear power plant in Argentina. The developed tool includes three-dimensional spatial neutron kinetics, a channel-level model of the core thermal hydraulics with subcooled boiling correlations, a one-dimensional model of the primary and secondary circuits including pumps, steam generators, heat exchangers, and the turbine with all their associated control loops, and a complete simulation of the reactor control, limitation, and protection system working in closed-loop conditions as a faithful representation of the real power plant. In the present paper, a description of the coupling scheme between the codes involved is given, and some examples of their application to Atucha II are shown.
1. Introduction
Atucha II is a pressurised heavy water reactor that is currently under construction in Argentina. It will provide a net power of 700 MW when connected to the grid. This reactor has some particular design features that obliges the computational model of the plant to take into account certain phenomena with a high degree of detail in order to retain the basic physics that determine both its steady-state and transient characteristics.
Considering that the reactor is a prototype, it is convenient to analyse the behaviour of the plant under certain operational transients and optimise its response by selecting proper gains and constants in the control system. In order to achieve this requirement, a calculation suite that couples different specific computational codes including a fine nodalisation of the core—in both neutronic and thermal-hydraulic aspects—was developed.
The engineering tool includes the ability to model three-dimensional neutronic kinetics coupled with a full thermal hydraulic representation of the core, up to an individual coolant channel description. The whole model is completed with a one-dimensional mathematical description of the plant, including the primary and secondary circuits with particular equations for pumps, steam generators, and the turbine. It also simulates the in-core and ex-core neutron detectors, whose signals are coupled to a software implementation of the reactor control, limitation, and protection systems. The tool is based on coupling specific codes that run in the same host computer and exchange information using shared memory objects. As an overall result, this suite of individual codes gives an autonomous tool that is able to compute transients near the operation design point with a high degree of detail using proved and validated nuclear calculation software within a reasonable execution time.
2. Plant Description
The plant is cooled and moderated by heavy water. Its core consists of 451 vertical fuel channels, contained in an in-vessel moderator tank. The coolant channels are surrounded by the moderating heavy water, which is enclosed in the moderator tank and is cooled by a separate system. For reactivity reasons, the moderator is maintained at a lower temperature than the reactor coolant. This is accomplished by extracting the moderating water from the moderator tank, cooling it in the moderator coolers, and feeding it back to the tank.
During full-load operation, 95% of the total thermal power is generated in the fuel, and the remaining 5% in the moderator, because of the neutron slow-down. Additionally, approximately 5% of the thermal power is transferred from the coolant to the moderator through the coolant channels walls, due to the temperature difference between the systems. The heat removed from the moderator is used for preheating the steam-generators feed water. The reactor coolant system and the moderator system are connected by pressure equalisation openings of the moderator tank closure head. Therefore, the pressure differences in the core between the primary coolant and moderator systems are comparatively small, which results in the thin walls for the coolant channels. This allows attaining a very high burn-up. Furthermore, the connection between the reactor coolant system and the moderator system permits the use of common auxiliary systems to maintain the necessary water quality.
Various methods are applied to control reactivity and thus the power output of the reactor. The reactor contains “black” (Hafnium absorbers) and “gray” (steel) control elements. These control elements are used to control the reactivity and the power distribution, to compensate the build-up of Xenon poisoning following a reactor power reduction, and to shut down the reactor. The reactivity value of all control elements is sufficient to shut the reactor safely down into a subcritical state. A second independent shutdown system is provided, which injects boric acid into the moderator tank. The reactivity can be also controlled by varying the moderator temperature within a certain range, which is advantageous for some operating modes.
The secondary loop circulates light water, and two U-tube steam generators feed a two-stage turbine coupled to a four-pole 21 kV generator. In normal operation, the thermal power is 2161 MW and the gross electrical power is 745 MW. Table 1 presents some of the main design parameters of the plant, Figure 1 depicts the primary circuit and Figures 2-3 show the nuclear reactor core. Coolant enters from the two cold legs to the reactor vessel and then flows through the downcomer before the two loops converge into the lower plenum. Mass flow is distributed differentially between the 451 channels by means of restriction nozzles that provide five different head losses. This way, five distinct hydraulic zones are configured (Figure 4) and peripheral channels, where the fission power is relatively low, are cooled by a smaller mass flow of heavy water than that corresponding to central channels. Consequently, the outlet temperature is almost the same for all the channels, resulting in a reasonably flat temperature distribution in the upper plenum.
3. Computational Codes
Atucha II has some unique characteristics that should be considered in detail when modelling the power plant for operational transient and accident analysis. One of such intrinsic features is that, on the one hand, the reactor has a slightly positive void reactivity coefficient and, on the other hand, subcooled boiling is expected to occur at some axial height in certain channels. To correctly model this situation, not only is a thermalhydraulic code capable of handling subcooled bubble generation needed, but also a detailed nodalisation of individual channels within the core. For if the channels were somehow averaged in order to decrease the complexity of the mathematical model, no subcooled void fraction would appear as a result. Another particular issue this reactor has is the fact that the control rods are slanted with different entrance angles. As the in-core instrumentation is arranged in lances inside the moderator tank, a control rod may locally disturb the flux near a neutron detector perturbing the signals that feed the reactor control system. A detailed three-dimensional neutronic model is needed in order to analyse the behaviour of the power distribution control system when a control rod is inserted near an in-core detector. Also, because of the heavy-water moderation, the reactor core is relatively big, and important xenon-instabilities may be induced during certain operational actions that the power regulation system ought to damp and control.
As a result, the analysis tool has to be able to handle spatial neutron kinetics, individual channel representation, and a full reactor power control, limitation, and protection systems simulation. The solution found is based on three specific calculation codes.(i)DYNETZ: an Atucha II ad hoc plant simulation code that includes both the primary and secondary circuits, and a simplified version of the reactor control and limitation logics [1]. (ii)PCE: an Atucha II ad hoc spatial kinetics core code [2]. (iii)RELAP: a well-known system thermalhydraulic code [3].
The developed tool is called DyPRA (DYNETZ-PCE-RELAP Acoplados) and is composed of some slight modifications of the listed codes above that allow them to exchange information and synchronise with other processes using shared memory segments and semaphores, respectively. The whole plant, with the exception of the coolant channels and the fission power distribution, is modelled by DYNETZ. The channels are simulated by 451 instances of RELAP, and the spatial neutron kinetic is computed by PCE. All these codes exchange information by means of shared memory, and they advance one integration step, one at a time, by operating on shared semaphores.
Each code by itself has been extensively validated by their developers as well by comparing code predictions with actual operational transients in other nuclear power plants. A thorough validation program is being specified to guarantee the validity of the individual code results when running under the proposed coupling scheme.
When called, DyPRA forks and executes the individual codes as daughters (as depicted in Figure 5), and afterwards it sets the shared semaphores to the proper red/green condition as the numerical integration proceeds (Figure 6). The information exchanged by the three codes is listed in Table 2. Once the simulation finishes, the daughters are politely terminated and the output is properly formatted for the analysis stage that may include key variables time plots, step-by step parameter spatial distribution, and logical signal plots or even geometrical representations of distributions of properties inside the core. Section 4 shows examples of the kind of information that can be extracted to analyse and optimise the behaviour of the reactor under certain transients.
A noteworthy feature of DyPRA and its related tools is that the whole project was designed so that the output analysis process may be done under any operating system using almost any text processing or data plotting program. Nevertheless, special effort was made to assure that there is at least one combination of these extra programs fully composed of free software. Some recoding of the calculation codes was needed in order to fit the standards and compile with the GNU compiler collection [4], and a small amount of output reformatting was introduced so that free software may be used as postprocessing tools. The development of the suite was done under GNU/Linux using only free software for the required tasks, that is, free compilers, debuggers, text editors, scripting systems, and so forth. Even the related documentation was prepared using a free typesetting system.
3.1. DyPRA
The main executable program is in charge of setting up the simulation, forking and executing the actual calculation codes (extended versions of DYNETZ, PCE, and RELAP as explained below), and advancing the integration steps one by one. After the simulation is over, it terminates the launched processes, cleans up allocated resources, and organises the outputs into subfolders. The standard output of the individual calculation codes is redirected to other files for later analysis. While running, DyPRA shows in the running terminal the reactor power and the instantaneous position of the control rods to give the user a prompt feedback about the evolution of the simulation.
The suite currently runs on GNU/Linux architecture, mainly because the individual codes consume a lot of resources (CPU, RAM, and hard disk) and the calculation cluster installed at the Atucha II site runs this operating system. On the other hand, the Linux kernel provides an efficient set of system calls for interprocess communication [5] that is the natural framework for the proposed kind of coupled calculation scheme. Nonetheless, with slight modifications it may run under other platforms as well.
It is possible to interrupt the execution by either pressing Ctrl-C or sending a TERM or INT signal to DyPRA. The program catches the interruption and then politely terminates its daughter processes, flushes buffers, removes shared locks, and frees allocated resources. This way, if the user detects a mistake in the simulation input, the execution may be rerun quickly and without the need for any manual intervention to clean up the mess that an interrupted simulation usually is.
Either manually interrupted or already finished simulations may be continued by means of the restart capabilities developed. Even though each individual code has its own restart mechanisms, a synchronisation scheme had to be developed to assure that at the times when a restart record is desired by the user, every code saves to the disk all the internal information needed for a successful continuation. Simulations may be restarted either from the last written record or from intermediate time steps. This last functionality may be used to insert disturbances at certain times and to compare the results with the original unperturbed case.
3.2. DYNETZ-451
DYNETZ [1] is a computer code developed in the 1980's by the original architect designer of the reactor, the former German KWU, to simulate the plant behaviour under certain design basis transients. It has a one-dimensional thermalhydraulic nodalisation of the whole plant that employs the homogeneous model to handle potential two-phase flow. The fission power is computed using the point kinetics equations with variable reactivity coefficients. The core is modelled as four average channels with six axial nodes. The code includes special models for the pressuriser, steam generators, turbine, pumps, and other particular equipment. The control and limitation logic is partially modelled because the point reactor equations do not provide certain information about changes in the power distribution or the ability to model differences in the in-core neutron detectors that form the input to the reactor control logic. Moreover, individual control rod movement is not implemented, and the eighteen control rods are represented by six groups of three rods that move continuously inside the core and provide a net reactivity that depends only on the bank axial position.
Nowadays, with the availability of spatial kinetics, it is possible to incorporate into the code a full description of the control and limitation systems, taking actions depending on the neutron detectors readings and coping not only with overall power changes but with its distribution as well. With this respect, it is necessary to allow individual control rods to be moved. Moreover, it is desired to model the actual rod movement using the standard three-coil displacement mechanism that results in discrete steps both in space and in time.
On the other hand, neither the homogeneous model nor the average-channel technique allow one to take into account the effects of subcooled boiling in some of the central channels. In order to have a second code to model the core, a method of combining the calculation of the thermalhydraulic conditions of the primary circuit by different codes had to be designed.
All these modifications, plus the ones needed to exchange information using shared memory segments and the ability to synchronise with shared semaphores, resulted in a code called DYNETZ-451. This program is the one that is launched by DyPRA when run. It waits until the parent sets the proper semaphore to green, reads information from the other codes (power distribution and initial thermalhydraulic core conditions), and computes the steady state. It sets another semaphore to green to indicate DyPRA that the steady-state conditions were computed and then waits again before advancing one integration time step. DYNETZ lets RELAP compute the hydrodynamic conditions in the core by setting the pressure in the lower and upper plenums. It then modifies its internal core head loss to match the mass flow calculated by RELAP and thus is able to compute the whole primary system as a closed loop.
As DYNETZ-451 is a standalone program that may be called independently (although it makes no sense without the corresponding kinetics and core thermalhydraulics codes), DyPRA may be told not to launch DYNETZ-451 so it may be run from a debugger. This way, internal variables may be monitored and particular control actions can be tracked back to their trigger conditions. This debugging capability is enhanced by the fact that the whole suite usually runs in a dedicated calculation machine while the debug can be done remotely from the user's workstation within a friendly graphical interface running under any modern operating system, that does not need to be GNU/Linux.
3.3. PUMITACPL
The neutronic characteristics of Atucha II present some challenges to the state-of-the-art spatial kinetics code. On the one hand, the separation of the moderator from the coolant in this type of heavy water reactors makes the utilisation of light water reactor codes useless for coupled calculations. On the other hand, the fact that the control rods are inclined with respect to the vertical poses an extra nodalisation problem that cannot be handled by most of the industry standard codes. Both existing Argentinian power reactors (Atucha I and Embalse) use PUMA [6] as the basic code for neutronic calculations. This program is able to deal with very general reactor problems with a great variety of nodalisations and characteristics, and it has also been validated against experimental measurements in different reactors. It can also handle LOGO instructions to further refine the reactor computations or to add particular features to the code. However, PUMA is a very big program that presents some difficulties if one desires to couple it with external codes using shared memory. With the objective of having a calculation tool that is simple but at the same time retains the basic physics underlying the core of Atucha II, the code PCE [2] was developed by CNEA. This software, also known as PUMITA, was coded from scratch in Fortran 90 (the original PUMA language is PL-1) and even though it has some flexibility, it was designed ad hoc to solve only Atucha II particular transients. PCE was validated against PUMA for several kinds of calculations for Atucha II [7].
Although PCE does not have external coupling capabilities, and the input only allows the user to enter homogeneous thermalhydraulic parameters, the internal matrices memory design included the ability to interpolate the macroscopic nuclear cross-sections of the materials using distributions of properties (temperatures, densities, boron concentration, etc.). Thus, the modification of the software to allow shared memory coupling and semaphore synchronisation with external processes did not require major changes. Besides the coupling capabilities, improved in-core detector simulation, control rod dilution effects corrections, some memory optimisations, and parallelisation using OpenMP [8] were introduced. As a result, a new neutronic tool called PUMITACPL [9] was created. It may be run standalone as an enhanced version of PCE, called from within the DyPRA suite or executed with another particular coupled calculation scheme. In fact, PUMITACPL was coupled to a few ad hoc codes to study some neutronic effects. For example, the numerical dilution of the control rods incremental cross-sections in the finite difference formulation of the diffusion equation was studied with a coupled scheme in such a way that the control rods were moved externally. In fact, after this analysis, two correction methods were proposed [9]. The PUMITACPL distribution package also includes examples to make frequency response studies by coupling the main code with a simple program that oscillates the position of certain control rods in time.
If a coupled simulation is desired, a coupling file must be provided. This file should be a plain text where the user defines by means of some keywords the names of the shared memory object where PUMITACPL should read, for example, the instantaneous positions of the control rods and the temperature distributions inside the core, and the objects where to write the power distribution and the neutron detectors readings. Also, the names of the semaphores to wait for and to set to green after advancing one time step should be given in the coupling file. For a DyPRA coupled calculation, this information is set up by a collection of scripts consistently with the shared objects defined in the DYNETZ-451 source code and in the coupling files of RELAP5CPL, as explained below.
3.4. RELAP5CPL
RELAP5CPL [10] is an extension of RELAP5/Mod3.3 [3] designed to allow the code to import and export internal variables from and to an arbitrary number of shared memory objects. One or more semaphores may be used to synchronise the calculation with the external coupled processes. The basic idea is that, if provided from the command line, RELAP5CPL may parse a coupling file after the problem input file. The coupling file should contain keywords that define imports and exports. An import or export is a combination of a shared memory identifier and a list of RELAP's internal variables, plus the names of the associated shared semaphores to wait for or to set into green. The extension is designed to allow certain flexibility in the coupling scheme, that is, variable number and size of shared segments, zero or more semaphores, synchronisation only after either a fixed number of steps or a fixed time increment, ability to terminate the simulation via shared memory, and so forth.
In principle, any variable stored in the fa array [11] may be read from or written to a shared memory segment. This way, boundary conditions, power distributions, control actions, and even heat transfer coefficients may be set externally whereas any of the results obtained by RELAP (usually temperatures, qualities, densities, pressures, mass flows, etc.) may be made available for other codes to utilise them. Imports and exports are processed in the order specified in the coupling file after the call to hydro and before the call to rkin in RELAP's tran routine. Not only were some of the original routines slightly modified, but also extra functions were coded (in C instead of Fortran 77) such as the coupling file parser and the shared object access routines.
The extension was designed to couple RELAP to a wide variety of external processes, such as neutronic codes, control codes, particular models, and even other instances of RELAP. Consider for example the case where a fine nodalisation of a whole plant does not fit into a single RELAP input file, but both the primary and secondary circuits as standalone loops do. Then, two coupled RELAPs may be run exchanging the conditions in each side of the steam generators via shared memory. RELAP5CPL gives a straightforward solution to implement this scheme by constructing two simple coupling files. Another example of the usefulness of coupling RELAP with external codes is illustrated in reference [12].
For the computation of DyPRA transients, a bash script generates the 451 distinct inputs from five templates (one for each hydraulic zone), and 451 coupling files that utilise shared segments and semaphore names consistently with the ones defined in DYNETZ-451 and in the coupling file of PUMITACPL.
4. Sample Results
As an engineering simulator designed to analyse the behaviour of a particular nuclear power plant, DyPRA allows a wide variety of results to be obtained. A very important part of the suite is the development of powerful postanalysis tools. To illustrate what is DyPRA capable of, a brief survey of example results is shown in this section.
Every transient simulation must start from a restart record, either one generated by a previous run or one containing the steady-state conditions of the plant. The steady-state restart registers are constructed from some initial guesses by separately coupling the calculation codes until they all converge to a solution that may be used as the initial guess for the full-coupled calculation without introducing numerical instabilities. The transient simulation may be a run to study the design point conditions or to include either operative changes or external disturbances.
Once a simulation has finished, the outputs generated by the calculation codes may be parsed and analysed to study different aspects of the calculated evolution. Figure 7 shows the thermal power distribution in the reactor core as a three-dimensional contour maps in the 451 channels. Figures 7(a) and 7(b) illustrate the typical distribution in normal operation at full power with the control rods in its nominal position. Individual channel or control rods may be turned on and off, and the whole model may be cut, zoomed, or rotated as desired by the user, giving the ability to study qualitatively the influence of the control rods in the power distribution. This instantaneous representation of the core may be advanced in time as desired to analyse transients, and the graphic tool may generate video files. For example, Figures 7(c) and 7(d) show how the distribution changes under a reactor trip because of the control rods inclination angles. It can be seen that even though the total power is smaller, the residual power distribution is very different from the total power at normal operation.
(a)
(b)
(c)
(d)
Asymmetries in the flux distribution caused by control rods movements or xenon poisoning can be graphically studied with the aid of the representation shown in Figure 8. There, two three-dimensional views of the thermal flux distribution in the plane located at the half of the active zone can be seen. Again, the views may be rotated or scaled, and transient videos may be constructed.
(a)
(b)
The fine thermalhydraulic nodalisation allows to carefully study not only the steady-state distributions but also the expected conditions under certain transients. For example, taking advantage of the reliability of RELAP's two-phase models, the overall pressure may be reduced a few bars giving rise to more void fraction in the core and a whole new spectrum of plant dynamics. Also, the outlet temperature profile may be easily investigated (Figure 9) and compared to the predicted original hydraulic zones design. Moreover, having a detailed coolant temperature distribution (Figure 10) makes it possible to compute more accurately the heat transferred to the moderator by radial conduction through the coolant channels walls.
One of the main reasons to launch the project that resulted in the development of DyPRA was the necessity of having a way to verify the original design of the reactor power control system. Even though originally DYNETZ had routines that modelled most of the control systems, the physical models employed did not allow a full scope simulation of the control logics. For example, consider Figure 11 that presents the evolution of the indication of the four ex-core ionisation chambers that are used by the reactor power control system to form the actual power value during a 100% to 80% power ramp. The effect of the continuous increase of the xenon concentration that gives rise to a wavy form of the neutron flux can be seen. The control rods are being extracted to overcome the negative reactivity caused by the overall power reduction. Using point kinetics equations, it is impossible to take into account the potential differences that may rise between the different chambers due to spatial kinetics effects or thermalhydraulic distribution changes during this kind of transients. The incorporation of three-dimensional neutronic models and best-estimate thermalhydraulic nodalisation of the core allows the analysis of the behaviour of the power regulation system under more realistic conditions.
Another example of control system optimisation is shown in Figure 12, where a sequence of results obtained during the design of an improved logic to reduce frequency of rods movements in normal operation can be seen. Because of the neutronic characteristics of the reactor, the core is never strictly critical. It goes periodically from subcritical to supercritical according to the insertion or withdrawal of the regulation control rods. Besides optimising the set of parameters of the control system (gains, dead-band sizes, instrument calibrations, etc.) DyPRA allows to propose different regulation logics. Figure 12(a) shows the difference of the reading reported by one of the ex-core ionisation chambers (the one closest to the corrected thermal power) with respect to the power set point versus time for two thousand seconds of normal operation at full power using the original control rod movement logic. The total power is bounded inside a given dead band. It is interesting to note that because of the positive feedback effects, sometimes the reactivity inserted by one rod step is not enough to reach the other limit of the dead band. Figure 12(b) shows the same reading but inserting an extra first-order lag in the generation of the rods set point giving as a result fewer movements. However, this case presented an asymmetry between power defect and power excess because of the azimuthal power distribution control. Additionally, as there are now less rod insertions than before, the effect of the thermalhydraulic feedbacks is bigger, and sometimes sequential movements are needed in order to change the sign of the overall reactivity. Figure 12(c) was obtained by suppressing the azimuthal distribution control, and, finally, Figure 12(d) presents both the extra phase lag in the set point generation and an improved version of the power distribution control. The overall result from (a) to (d) was to reduce the frequency of control rod movements.
(a)
(b)
(c)
(d)
In order to exploit one of the main advantages that the usage of DyPRA gives, that is a detailed information of all the signals generated by the instrumentation and control systems, a postprocessing spreadsheet was developed. Screenshots of this data analysis tool are shown in Figures 13, 14, and 15. When a simulation finishes, the data analysis spreadsheet can import the output files and display time plots of the main control variables of Atucha II, such as thermal and generator power, control rods set points and actual positions, ex-core and in-core detector readings, and so forth. The information arranged in this way simplifies the comprehension of the instantaneous plant conditions and helps the engineer to analyse the actual response back to the causes that produced it. It also helps to detect inconsistencies, either in the code programming or in the expected results.
All the figures presented in this sections were produced by processing the outputs generated by DyPRA with free graphical postprocessing and plotting tools. Eventually, some small programs and scripts were developed to adapt the output of the calculation codes to fit the format of the free available data analysis applications.
5. Conclusions
DyPRA is a novel tool that allows the computation of Atucha II operational transients taking into account several distinctive features of the reactor by coupling particular calculation codes. It consists of an engineering simulator that allows the study and analysis of conditions in the whole plant, with special focus on the reactor core. It was developed by extending an industry standard thermalhydraulic program, an ad hoc plant simulator and an indigenous three-dimensional kinetic code. A fourth code that synchronises the integration of the governing equations along with some scripts to automate the process of input generation and a collection of postprocessing scripts and tools were written from scratch. The suite is able to represent several particular features, including the effect of subcooled boiling on the transient power distribution and vice versa, disturbances in the in-core neutronic instrumentation, nonuniform xenon-poisoning and its control mechanisms, among others. Moreover, DyPRA provides a framework to study and analyse improvements to the reactor control logics, and a very convenient way to optimise its adjustable parameters such as gains, dead-band sizes, and time constants. In fact, a great deal of knowledge has been gained not only by analysing DyPRA's outputs but also during the development itself.
The proposed coupling method based on information interchange using shared memory segments is very efficient, as the access speed is several orders of magnitude greater than that of files for example. Thus, the overhead required by the coupling scheme is negligible, and the running time is defined only by the complexity of the mathematical equations to be solved. As the core thermalhydraulic nodalisation is quite fine and the spatial kinetics require a lot of computational effort, the overall running time is in the order of ten to twenty seconds of calculation time per real problem second, depending mostly on the nature of the neutronic transient being studied.
Even though still under development and improvement, one year after the development started, the proposed ideas converged to an usable version that already threw some interesting results about the plant’s inherent dynamics and its regulation. The results are constantly being compared and validated with existing engineering results, obtaining reasonable coincidences in common cases. The definite validation will take place during the plant start-up tests. Actual transient calculation results with their related discussion and comparisons of numerical predictions with measurements taken from the real plant will be the subject of future papers.