- Email: [email protected]
Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment
Annals of Nuclear Energy xxx (xxxx) xxx
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment R.A. Busquim e Silva a,b,c, K. Shirvan b, J.J. Cruz a, R.P. Marques a, A.L.F. Marques a, J.R.C. Piqueira a,⇑ a
Polytechnic School of University of Sao Paulo, Av. Prof. Luciano Gualberto, 380, Sao Paulo, Brazil Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA, USA c Navy Technological Center in Sao Paulo, Av. Prof. Lineu Prestes, 2468 Sao Paulo, Brazil b
a r t i c l e
i n f o
Article history: Received 18 May 2019 Received in revised form 22 August 2019 Accepted 2 October 2019 Available online xxxx Keywords: NPP control PWR Digital instrumentation MATLAB PARCS RELAP
a b s t r a c t The operation of current nuclear power plants (NPP) and the deployment of advanced nuclear reactors rely on digital instrumentation and control (I&C) and information technology systems. The growing use of digital systems implies considerable security and safety challenges, as the usual standard codes used for designing, licensing, and analyzing accidents over the last few decades have been limited in terms of I&C features. This study presents an innovative method to couple PARCS/RELAP5 with the high-performance language MATLAB/Simulink. As a proof-of-concept proof, a 2772 MWt pressurized water reactor modeled using the coupled package PARCS/RELAP5 is simulated. A supervisory system is implemented that allows MATLAB/Simulink to oversee the nuclear codes for exchanging online data. In addition, I&C high-level computing tools are added to the NPP model, endorsing the application of ‘‘on-the-fly” decision logics. The coding strategy does not change the nuclear binary source codes, as the PARCS/RELAP5 codes have been validated for design and licensing. By running open- and closedloop cases for control rod assembly movement, we found that PARCS/RELAP5 capabilities for digital I&C modeling and simulations can be extended by adding MATLAB/Simulink flexibility. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction In the last four decades, computer nuclear codes have been extensively used by industry licensees and regulatory boards to simulate nuclear power plants (NPPs) and fuel cycle facility (FCF) systems. Simulations include design assessment and safety/ accident analysis. In addition, the nuclear industry foresees a massive increase in the use of digital instrumentation and control (I&C) systems, while also being wary of the emergence of new and sophisticated cyber threats. These factors have motivated the development of advanced coupled computational tools for improved management, cyber security assessment, higher process performance, and cost reduction [1,2]. The growing use of digital I&C systems raised one major security and safety issue related to the possibility that unauthorized codes may change the NPP behavior as a result of a cyber attack. The development of tools that allow the assessment of the impact on the facility, from the thermo-hydraulic (TH) and neutronic (NK) point of view, provides significant opportunities for conducting cyber defense experiments. By using MATLAB/Simulink, real ⇑ Corresponding author. E-mail address: [email protected] (J.R.C. Piqueira).
control equipment can be interfaced with the supervisory system to determine the consequences of sabotage by exploitation of security vulnerabilities that may result in loss of confidentiality, integrity and availability. In this study, an innovative computational method to interconnect two widely used nuclear codes is presented. The thermalhydraulic system code RELAP5 (R5) and the neutronic nuclear code PARCS 3D (P3D) are connected by using a well-known language for technical computing, namely, MATLAB/Simulink (MATLAB). The supervisory system that is developed extends the capabilities and features of the parallel virtual machine (PVM) package and P3D/ R5 I&C by coupling them with MATLAB. In addition, by accessing and managing the R5 and P3D input, output, and restart files during code execution, basic and advanced MATLAB toolboxes and nested functions may be applied to NPP models. The supervisory system was developed only for I&C modeling and simulation and not for actual I&C systems on NPP. For research purpose, the supervisory system was coupled to a 2772MWt pressurized water reactor (PWR) Babcock & Wilcox Co. (B&W) as modeled in [1,3]. Moreover, this simulation strategy enhances the coupling features without changing the nuclear binary source codes, which is a major innovation given that nuclear codes have been validated
https://doi.org/10.1016/j.anucene.2019.107098 0306-4549/Ó 2019 Elsevier Ltd. All rights reserved.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
2
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
during the last decades based on multiple benchmark activities for designing and licensing purposes. With the source code remaining unchanged, reliability of the simulation results is guaranteed. The implementation described herein includes: - Use of high-level I&C MATLAB modeling capabilities to add flexibility to the simulation, considering that P3D/R5 has limited supervision tools compared to MATLAB; - Application of on-the-fly decision logic based on P3D/R5 simulation results (i.e., within each MATLAB time step); - Ability to assess separately the decision logic (MATLAB) from the NPP behavior and from the TH and NK points of view; - Handling of the input files by writing new information before calling P3D/R5; and - Saving exchanged data and output information for further analysis. The next section presents some information about the simulation codes. After these preliminaries, two application sections describe how to use the developed tools; one section presents the supervisory system and the other the simulation strategies. Results and discussion conclude this work.
2. Simulation codes Since the 1970s, nuclear codes have been developed to simulate TH and NK NPP conditions. The analysis of NPP behavior usually begins with a TH system code that may include point kinetics (PK) or a more complex model to estimate the thermal power generation within the core [4]. The NK codes are mostly based on nodal expansion methods derived from neutron diffusion theory, where the macroscopic cross sections depend on feedback parameters such as fuel temperature, moderator density, and moderator temperature [5]. The use of nuclear TH and NK coupled codes allows a more complete NPP system analysis than when using standalone codes, which includes the feedback effect, core NK behavior, and TH dynamics. In addition, coupling tools allow more realistic simulation of normal operation (steady-state) and accident (transient) scenarios. This work uses the P3D [3] and R5 [6], which are advanced computer tools for TH and NK simulation, and the computing environmental software MATLAB [7] as a programming platform for implementing the supervisory system that drives the coupled simulations.
2.1. PARCS neutronic code The Purdue Advanced Reactor Core Simulator (PARCS) is a 3D reactor core simulator spatial kinetics code that solves the steady-state and time-dependent multi-group neutron diffusion and low-order transport equations in three-dimensional Cartesian geometry [8,9]. Furthermore, it has an embedded simplified TH solution, which accounts for the feedback effect when running in standalone mode, as PARCS uses an internal temperature calculation scheme based on input power. P3D is computer code that performs standalone calculations for light water reactors (LWRs), or it can be coupled to R5 using the PVM interface or directly to the TH system code TRACE-M [9]. P3D has been used for NK calculations by the United States Nuclear Regulatory Commission (NRC) [8]. Although P3D has few graphical user interfaces (GUIs), it works primarily by means of input and restart files for scenario configuration and output files for presenting results. Our simulations used PARCS version 2.71d.
2.2. RELAP5 system thermal-hydraulic code The Reactor Excursion and Leak Analysis Program (RELAP) is a tool for analyzing loss of coolant accidents (LOCA) and system transients [10]. R5 (i.e., version 5 of RELAP) has been used for licensing and regulation, evaluation of operator guidelines, and plant accident analysis. The code was developed to provide the firstorder effect needed to accurately predict system transients and to be sufficiently cost effective to enable sensitivity and parametric studies. Although it was initially designed for nuclear applications, the latest versions can be used for nuclear and non-nuclear applications involving vapor, liquid, non-condensable gases, and nonvolatile solutions [4]. R5 calculates the expected evolution of TH properties during steady state or transients, including the interaction among system components by using an internal temperature calculation scheme based on power input. It implements a two-fluid and sixequation model to simulate the two-phase coolant behavior and, consequently, the system’s TH behavior. However, the PK model inherent to R5 cannot consider multidimensional effects caused by spatial variations of feedback parameters. R5 also has few GUIs. However, it runs by means of input, output, and restart files. Our simulations used RELAP5/Mod3.2 version. 2.3. MATLAB computing environmental code The Matrix Laboratory (MATLAB) is a high-level numerical computing environmental software developed by The MathWorks Inc. for technical computing and programming not only for control engineering but also for data acquisition, data analysis, mathematical modeling, algorithm development, parallel computing, and web development. It has been in operation since 1984 [7]. MATLAB uses specific languages and codes controlled by its command window. MATLAB also allows interfacing with codes written in other languages including C, C++, Java, and FORTRAN. Hardware connections with field-programmable gate arrays or programmable logic controllers are available. It includes numerous toolboxes and packages such as Simulink, which adds graphical multi-domain simulation and model-based designs for dynamic systems, and MuPAD, which allows symbolic computing capabilities. 3. Supervisory system This work implements a supervisory system that couples P3D, R5 and MATLAB. This system oversees the nuclear codes online exchange of data; adds I&C high level computing tools within the NPP model; and endorses the application of ‘‘on the fly” decision logics. Therefore, the designed three codes coupling strategy is the heart of the supervisory system. The coupling strategy applied in this study extends the PVM coupling features without changing the P3D/R5 binary source codes. This enables the simulation results obtained using the developed supervisory code to be used for designing and licensing purposes. Through the access and management of R5 and P3D input, output, and restart files, basic and advanced MATLAB toolboxes and nested functions can be applied to NPP models specially developed for P3D and R5. The coupling strategy among the three codes includes: 1) The development of R5 and P3D standalone input and mapping files, which allows for proper TH and NK coupling; 2) The development of routines enabling the exchange of information between the coupled P3D/R5 and MATLAB in such a manner that MATLAB oversees the P3D/R5 exchange of data
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
within a predefined time step established for each simulation scenario; and 3) The development of functions to write and change R5 and P3D input files cards during code execution based on I&C decision rules.
3
- a control decision logic; - predefined R5 and P3D input file templates; and - a change of R5 and P3D input file cards during the simulation. Fig. 1 presents the time precedence among the codes. Basically, the simulation runs as follows:
3.1. P3D/R5 coupling The P3D and R5 codes work together to exchange information through a general interface [3], which manages the data flow between both codes by using the PVM. To guarantee the proper exchange of information between the codes, the NK-node and TH channel structures must match each other. This is accomplished by generating an input file named Maptab. The adequate generation of the Maptab file and correct determination of the mapping weight factors are necessary for obtaining accurate coupled code solutions. Moreover, the core nodalization in the radial and axial directions and their correlation have a considerable effect on determining the local core TH parameters and power distribution during transients. The mapping determines where to set NK node power in the TH channel (for both direct moderator heating and a heat structure fuel source) and where to assign TH properties in the NK node. Therefore, the coupling of the TH channels with the NK nodes in radial and axial directions depends on the adequate assignment of mapping weights between the two meshes. In this study, the mapping weight is based on the volumetric geometric fraction that correlates the TH and NK mesh in the radial and axial directions. The user is responsible for the assignments of the weighting factors. Then, a MATLAB script generates the Maptab file following the steps suggested in reference [11]. The TH conditions provided by the external systems code R5 at every time step in the NK spatial kinetics code P3D for the simulated case are the moderator temperature, vapor and liquid densities, void fraction, and averaged, centerline, and surface nuclear fuel temperatures [11]. These parameters are required at each NK node for the feedback calculation. The nodal/space dependent power distribution is then calculated by P3D and sent to R5, which calculates the heat conduction in the core heat structures. The TH code R5 addresses all the TH issues, and the NK code P3D spatial kinetics solves the NK challenges [1]. 3.2. P3D/R5/MATLAB coupling The PVM coupled P3D/R5 package and the MATLAB codes work together to exchange information through a leading-edge developed routine, implemented by means of a MATLAB system function (S-function) that manages the data flow between them. An Sfunction is a user coded dynamically linked subroutine that MATLAB can automatically load and execute. S-functions enable the code to work as built-in Simulink blocks [12]. The supervisory system uses the developed routine/S-function to stop/start and manage the input/output P3D/R5 files, and it uses the output from the same routine/s-function to implement the decision logic, which is applied to P3D/R5 closing the feedback loop during code execution. The supervisory system code allows the use of three nuclear codes in a single piece of software by means of: - an R5 and P3D spatial mesh that couples TH and NK using the PVM; - a predefined MAPTAB file that allows the correct matching between the NK node structure and TH channels (axially and radially); - a script that synchronizes the exchange of information; - on-the-fly access to the R5 and P3D restart and output files;
- The supervisory system code drives the steady-state run to generate the initial conditions for the transient simulations (restart files). Therefore, P3D/R5 runs within a Dt time step up to the predefined MATLAB DTSS; - During the transient operation, the supervisory system code drives the simulation based on a MATLAB DTTR time step that depends on the scenario being evaluated; - At each DTTR (DTTR Dt), the supervisory code calls P3D/R5 throughout the PVM and reads/uses the nuclear code restart information from the previous DTTR (or DTSS if it is the first transient supervisory time step); - The supervisory code determines the values and cards that must be written/changed in the input files based on the control decision logic; and - The supervisory code writes all changes in the input files and calls the nuclear P3D/R5 code again. The P3D/R5/MATLAB coupling strategy enables a more efficient and flexible I&C modeling and simulation for design, protection, and use of advanced control techniques. By using an S-function to access the coupled P3D/R5 package, the supervisory Simulink interface enables the incorporation of MATLAB algorithms and native toolboxes for instrumentation, control and simulation of dynamic systems. 3.3. NPP main modeling data The evaluation of the supervisory system performance is based on a set of NPP simulations and benchmarks. The NPP design is based on the 2772-MWt PWR B&W, which is well known from several studies published after the TMI Unit 1 accident. The main TMI NPP technical specifications can be found in the 1999 benchmark study released by the Nuclear Energy Agency from the Organization for Economic Cooperation and Development for mainstream line break [13] and references, such as [1] and [3]. Detailed information on P3D/R5 modeling can be found in reference [9,11]. The NPP core has a 17 17 assembly configuration. Sixty control rod assemblies (CRA) grouped into 7 banks (full-length control rods) are also employed. The core active height is 357.12 cm and
Fig. 1. Time step precedence among the codes.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
4
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
21.81cm reflectors on the top and bottom of the core are used, contributing to a total height of 400.74 cm. The position of the CRA insertion is given from the bottom of the lower reflector, where heights of 378.93 and 36.23 cm indicate completely withdrawn and completely inserted CRAs, respectively. The insertion is measured in steps of 0.35 cm (zero and 971 steps are equivalent to completely inserted and completely withdrawn CRAs, respectively). For NK modeling purposes, each radial node has the same flow area (FA) width of 21.81 cm, and each assembly is considered to be neutronically homogeneous. The FA has different U235 enrichments and burn-ups as well as a different number of burnable absorbers. The reactor is assumed to be in end of cycle (EOC), that is, 650 effective full power days (EFPD) or 24.58 GWd/MT. Axially, different compositions are used for each assembly, where each composition has a complete set of diffusion coefficients and macroscopic cross sections. The core is also divided axially in 28 NK regions. For TH purposes, the active core is axially divided into six fuel volumes and two reflectors. Radially, 18 TH channels and one regular channel represent the reflectors (channel 19).
4. Simulation strategy The main goal of supervisory system development is to enhance nuclear instrumentation and control modeling as well as simulation capabilities for future testing of digital I&C system reliability and redundancy for both safety and cyber security. For the purpose of this study, the performance of the supervisory system was assessed by assigning the control rod drive mechanism (CRDM) CRA positions by simulating the cases presented in Table 1. For these simulations, open- and closed-loop feedback supervisory setups were implemented. Moreover, the simulation results were benchmarked against P3D/R5 standalone runs. For comparison purposes, the nuclear code cards used the supervisory system post-simulation data. Prior to the steady-state runs, a set of P3D/R5 simulations were performed to identify the initial conditions of the control rod assembly that lead to criticality under steady-state conditions. For the same CRA axial positions (Banks 1 to 7), it was found that 50% power corresponds to 411 steps withdrawn and keff = 1.00083. In addition, 80% corresponds to 817 steps withdrawn and keff = 1.00079. During simulation, the supervisory code writes the adequate P3D cards into the input files, accesses the restart files, and manages the time step and output files. The R5 cards may also be changed on the fly depending on the simulation scenarios. Fig. 2 and Fig. 3 present the high-level Simulink block diagrams used for simulation. The plant block (blue) uses the S-function mechanism for extending the capabilities of the Simulink environment by means of compiled code dynamically linked to subroutines – all subroutines to access the P3D/R5 files are called under the Plant block. The Decision Logic diagram (green) contains decision logic for each one of the cases under study (e.g. the PI controller). The block diagrams in orange allow the use of data within the MATLAB workspace (e.g. for plotting).
Table 1 Cases under investigation. Case
Case Description
A B C D
CRA total insertion from 411 to 0 steps Power set point from 50 to 80% relative power CRA withdrawn from 411 to 817 steps Power set point from 50% to hot zero power (HZP)
4.1. Open-loop set up In the open-loop setup as presented in Fig. 2, the output of the plant (i.e. the relative reactor power, that is, variable ‘‘rktpow” in the plant block) has no effect on the input of the plant (i.e., the CRA position, that is, variable ‘‘bank pos cmd” in the Decision Logic block) that is applied to the core CRDM. No feedback information is used to adjust or control the reactor power. For simulation purposes, all CRA were assumed to be in the same axial position at the beginning of the simulation. 4.2. Closed-loop feedback setup Closed-loop systems compare the output information to the desired set point (i.e., the reference) and generate an error. This error is then used as input for a controller designed to have the output more closely from the reference. The design and implementation of the controller depends on the system being controlled. For the closed-loop feedback set up, Fig. 3 no CRA axial position is set in advance. During the simulation, the supervisory system measures the nuclear power ‘‘rktpow” (process variable) at each supervisory DTTR time step and sends its value to the Decision Logic block. The Decision Logic block receives the reference ‘‘rktpow_ref” (which may be set by the reactor operator, for instance), compares both values, and generates an error. Then, the controller calculates the CRA bank position ‘‘bank_pos_cmd” to be applied to the core CDRM. The controller implemented here was a proportional-integral (PI) type. A PI controller is a control loop mechanism widely used in the industry for minimizing the error e(t) over time between a desired set point and a measured process variable. It is named PI controller due to a proportional term (P ¼ K p eðt Þ) and and integral term Rt (I ¼ K i 0 eðsÞ dt). The P term takes into account a gain factor Kp that is proportional to the difference between the desired set point and the measured variable. The I term accounts for accumulated past values of the error. By integrating the residual error between the desired set point and the measured variable, it minimizes its residual due to cumulative values. A PI controller is commonly used in various practical applications even though it has many limitations, including an inability to provide optimal control (it may be tuned without direct knowledge of the process) and the use of constant controller parameters. However, it can be easily applied by knowing the dynamics of the plant. The PI input is the error e(t) and the output is the control signal u(t) applied to the plant according to the following equation:
Z uðt Þ ¼ K p eðtÞ þ K i 0
t
eðsÞ dt
Here, Kp allows the controller to make a change to the output proportional to the error. An improper tuning of the proportional term may lead to an unstable (too high) or slow sensitive (too small) response. In addition, Ki is the integral term that mitigates the steady state error as derived from the proportional term as it simultaneously speeds up the output response toward the set point. The integral parameter was set assuming an integral time Ti to minimize the output/set point offset and making K i ¼ 1=T i , and Kp was set to avoid an unstable output response. For all cases under investigation, Kp was taken as 5 and Ti was assumed to be 6. The proportional and derivative terms were determined by applying the Ziegler-Nichols method [14] followed by manual fine tuning. (It is worth noting that the PI implementation works in a continuous time framework, and the plant output is taken at each DTTR).
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
5
Fig. 2. MATLAB/Simulink open-loop interface.
Moreover, the NPP ran in steady state mode for 1210 s (all plots started in t = 1200 s for better viewing). Then, from 1210 to 1380 s, it ran in transient conditions. During the steady state, the coupled TH and NK solution converged within the 1200 s such that all its values remained almost constant. During the transient mode, the supervisory system operated as shown in Fig. 1. Moreover, the MATLAB DTTR time step was set to 1 s for all cases. P3D/R5 ran on a 103-s Dt time step. For comparison, a closed-loop simulation with DTTR = 0.2 s was also performed.
5. Simulation results and discussion The supervisory system adequately managed to access and change cards in the P3D/R5 input files based on the decision logic output. It also accessed data in the output and restart files, as shown Fig. 1. However, the supervisory system was able to drive the reactor dynamics. Therefore, its results were very similar to the P3D/R5 standalone runs. This supervisory system was designed to apply to more complex on-the-fly I&C tools as well as to offer a
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
6
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
Fig. 3. MATLAB closed-loop interface.
cyber security strategy for future applications by using MATLAB real digital system connectivity. While the cases under investigation were simple from I&C point-of-view, the model details are realistic enough to verify the computational methodology. For simulation purposes, the MATLAB DTTR time step was set at 1 s. Thus, the decision logic action was delayed by 1 s when compared to the P3D/R5 standalone runs. DTTR was defined by the user at the beginning of the simulation. The nuclear power behavior is explained by the CRA insertions/withdraws and the correspondent CRA reactivity trend. Table 4 summarizes the simulation results by presenting the standard deviation (SD) and mean between the P3D/R5 and supervisory reactor power.
5.1. CRA total insertion from 411 to 0 steps (Case A) The open-loop (see Fig. 2) axial CRA insertion from 411 to 0 steps withdrawn as shown in Fig. 4 was set in the Simulink decision logic block in which no feedback information was given. The power response after the total CRA insertion is shown in Fig. 5. As expected, the nuclear reactor power decreased as the control rod assemblies were inserted, and the CRA reactivity was negative (see Fig. 6). It can be seen that the reactor power remained under shutdown conditions after 1,360 s, while the CRA reactivity was expected to return gradually to zero.
Fig. 4. CRA bank position: Case A.
In addition, Fig. 7 and Fig. 8 show that the 3D core power distribution decreased as the CRA was inserted, which was expected, and it remained relatively uniform during the transient simulation.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
7
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
Fig. 5. Total power: Case A. Fig. 8. 3D power distribution t = 1320 s.
Fig. 6. CRA reactivity: Case A.
This behavior was similar for all cases under investigation. The 3D power distribution represents the assembly peaking factor (i.e., the assembly power divided by the average assembly power in the core, which is the reactor total power divided by the number of fuel assemblies). In addition, to compare results, the same pre-defined axial position was set in the coupled P3D/R5 package ‘‘move_bank” card. However, a modeling limitation existed in the card when the standalone simulation was run: the pair (time, bank position) could be repeated up to 10 times. However, the supervisory system runs P3D/R5 and writes new ‘‘move_bank” values within the MATLAB DTTR. Thus, the supervisory system can move the CRA 10 times per time step, which enables better control. Table 2 shows the bank position assumed for the P3D/R5 benchmark (up to 1200 s, whereby all CRA are 411 steps withdrawn and, after 1320 s, remain fully inserted). Fig. 9 and Fig. 10 show that no significant difference was found between the open-loop supervisory and P3D/R5 package benchmark runs. Fig. 11 and Fig. 12 present the power and CRA reactivity differences between the P3D/R5 standalone and supervisory runs. It can be seen that the differences were very small. The main differences were derived from the 1 s delay between the results and the user-defined MATLAB DTTR time step.
5.2. Power set point from 50 to 80% relative power (Case B) This case was evaluated by a closed-loop setup (see Fig. 3) assuming a 50 to 80% relative power set point. Fig. 13 shows that the application of the PI controller led, as expected, to an increase in nuclear power toward the reference. The control signal calculated by the PI controller, which was applied by the supervisory system to the CRDM, is shown in Fig. 14.
Table 2 CRA position: Case A.
Fig. 7. 3D power distribution t = 1210 s.
Time (s)
CRA Position (steps)
Time (s)
CRA Position (steps)
1200 1210 1220 1240 1250
411 411 372 300 261
1260 1280 1300 1320 1380
224 148 73 0 0
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
8
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
Fig. 9. Total power: Case A and P3D/R5.
Fig. 12. CRA reactivity difference.
Fig. 10. CRA reactivity: Case A and P3D/R5.
Fig. 13. Total power: Case B.
Fig. 11. Power difference.
Fig. 14. CRA position: Case B.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
5.3. CRA withdrawn from 411 to 817 steps (Case C) For Case C, the CRA bank position obtained (calculated) from Case B (see Fig. 14) was applied directly to the CRDM using the open-loop setup. The expected nuclear power trend followed a ramp from 50 to 80% nominal power. In addition, a P3D/R5 package standalone case was simulated by setting the P3D ‘‘move_bank” card according to Table 3. (This was a post-simulation data setup obtained from Fig. 14 by direct observation.) Then, Cases B and C as well as the P3D/R5 standalone run and reference were all compared. Fig. 15 and Fig. 16 show that, in addition to the differences derived from the MATLAB DTTR and Table 3 CRA position: Case C. Time (s)
CRA Position (steps)
Time (s)
CRA Position (steps)
1200 1210 1220 1240 1260
411 411 424 467 525
1280 1300 1320 1380 1400
590 660 732 806 817
9
P3D/R5 time steps as well as the discrete P3D/R5 package application of the ‘‘move_bank” card, all simulation results were quite similar. 5.4. Power set point from 50% to HZP (Case D) This case was also evaluated by a closed-loop reference from 50% power to HZP setup. Then, the axial position (postsimulation data) from the supervisory results was set in the P3D/ R5 package ‘‘move_bank” card, and a standalone simulation was performed. The results are presented in Fig. 17 and Fig. 18. As it can be seen, the results were very similar. The differences as presented in Fig. 19 and Fig. 20 were due to the supervisory system time step and P3D ‘‘bank_move” card that was limited in terms of the number of pairs (time, bank position). 5.5. Standard deviation, mean, maximum and minimum benchmark data To benchmark supervisory results against the P3D/R5 package standalone runs, Table 4 presents the SD, mean, and maximum and minimum values calculated over the difference (DP) between
Fig. 17. Total power: Case D. Fig. 15. Total Power: Case B.
Fig. 16. CRA reactivity: Case B and P3D/R5.
Fig. 18. CRA reactivity: Case D and P3D/R5.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
10
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
Fig. 19. Power difference.
Fig. 21. Total power: Case B.
Fig. 20. CRA difference. Fig. 22. CRA bank position: Case B.
Table 4 DP simulation data. Case
SD (%)
Mean (%)
Max (%)
Min (%)
A B C D
0.301 0.302 0.290 0.308
0.226 0.211 0.175 0.250
0.0 1.036 0.398 0.324
1.123 0.671 0.001 1.617
the reactor power obtained with the P3D/R5 package (open- and closed-loop) and supervisory system. It can be seen that the SD was close to 0.3% for all cases under simulation and that the mean remained close to zero (in the range 0:25 DP 0:2). These were acceptable results considering the dynamics of the plant. The maximum and minimum values were also suitable results considering the delay due to the MATLAB DTTR time step. 5.6. Simulation results for distinct MATLAB time step The application of smaller values of MATLAB DTTR may lead to better results. However, the computational time will increase.
For simulation purposes, a DTTR = 0.2 s discrete PI (parameters Kp = 3 and Ti = 5) was implemented and both results were compared. Fig. 21 and Fig. 22 show that no significant difference was found between the results for a total power increase from 50 to 80%. However, it must be noted that the appropriate MATLAB DTTR must be chosen based on the expected dynamics of the plant under the simulations performed. 5.7. Computation running time The coupled package P3D/R5 is a time-consuming computational code. Nevertheless, most computational time is consumed by the neutronic and thermo hydraulic calculation. The supervisory system coupling task includes the opening of different files for binary read access as well as formatting and writing data, which are also time-consuming. In addition, it uses feedback information as input for the decision logic. For further improvements to the supervisory system regarding its computational demand, the elapsed time can be seen in Table 5 and Table 6 for the cases under investigation (all runs were performed on a Dell Laptop with an Intel Core i7 7500U CPU @ 2.7 GHz running Windows 10 at 64 bits). The simulation time
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
11
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx Table 5 Time step DTTR = 1 s elapsed time (s). Case
Supervisory System (P3D/R5/MATLAB)
Nuclear Power from 50 to 0% Nuclear Power from 50 to 80%
Open-loop
Closed-loop
1.713 103 1.762 103
1.740 103 1.770 103
P3D/R5 Standalone
Ratio Supervisory System (P3D/R5/MATLAB) Open-loop to P3D/R5
Ratio Supervisory System (P3D/R5/MATLAB) Closed-loop to P3D/R5
1.490 103 1.491 103
1.15 1.18
1.17 1.19
Table 6 Time step DTTR = 0.2 s elapsed time (s). Case
Supervisory System (P3D/R5/MATLAB) Closed-loop
P3D/R5 Standalone
Ratio Supervisory System (P3D/R5/MATLAB) to P3D/R5
Nuclear Power from 50 to 0% Nuclear Power from 50 to 80%
3.081 x104 3.086 x104
1.490 103 1.491 103
20.68 20.70
Table 7 EKF cyber attack scenario. Core Conditions
CRA at t = 0
Attack Threshold
Cyber-Attack
Steady-State 100% Power
Fully/Partially withdraw
DPower>±10% nominal power D reactivity >±0.1 $ D moderator temperature >±10 °C
CRA goes from fully inserted to fully withdraw HMI compromised
was approximately 1.15 to 1.19 times higher for the supervisory system at DTTR = 1 s, and there was a greater computational demand of approximately 20.7 times for the smaller MATLAB time step of 0.2 s. This was expected, as the supervisory executed a greater number of open/write actions for a small time step. (For data consistency, the elapsed time did not consider the time when the simulation was not running or paused.) No significant difference was found between open- and closed-loop DTTR = 1 s simulations. The running time depended on the MATLAB DTTR time step, the number of restart block files desired, and the dynamics of the simulation. The additional time as a result of the MATLAB coupling may be reduced by retaining information in the CPU memory and by reducing the number of R5 restart files to one for each MATLAB DTTR. Although this was not the goal of this study, further implementation (see reference [1]) intends to minimize the computational time. 5.8. EKF as cyber security point defense option: a simplified application The use of nuclear codes, validated with well-established criteria and methods, provide nuclear safety analysis advantages over conventional non nuclear codes used for I&C modeling and simulation. On the other hand, specially designed I&C modeling and simulation software may add flexibility to potential users that may not invest resources to qualify new software for their reactor control system. Moreover, the use of the supervisory system provides significant opportunities for conducting cyber defense investigations due to the capability of assessing both digital system function and facility function impacts after a cyber-attack. This section presents a practical application of the stochastic-forecasting extended Kalman filter (EKF) algorithm developed in reference [1,3] integrated with the supervisory system, that provides a point defense option against cyberattacks. The EKF model is based on the point kinetics reactor equations, on the temperature feedback effects and on simplified core-coolant mass and energy balance equations. The EKF is an extension of the
Kalman Filter (KF) for nonlinear systems, which works on a stream of noisy input data to generate a statistically optimal estimate of the variable of interest. The KF is optimal linear estimator that minimizes the mean square error of the estimated variable if the noise has a Gaussian distribution. The EKF estimation is based entirely on the one-directional neutron flux (or reactor power). For simulation purposes, noise is added to the measured signal. The EKF algorithm produces estimates of the reactor power, moderator temperature and reactivity and compares these estimated values to their measured counterparts. Based on the detection criteria presented in Table 7, a control signal is sent to the NPP protection system. The scenario analyzed under this study considers: - The control rod drive mechanism is compromised during the NPP normal operation (100% power). - The control rods are take to ‘‘all rod out” position increasing reactivity. - The human-machine interface (HMI) is compromised by hiding the aforementioned action. - The total power, reactivity and moderator temperature HMI indications are not updated. - The EKF, using an one-directional/direct neutron flux signal, estimates the reactivity and the moderator temperature. - Detection criteria, designed to minimize false positives and missing alarms; - The detection criteria assume low and high thresholds for detection purposes. Fig. 23 presents the EKF implementation. The Decision Logic block implements the attack threshold shown in Table 7. The simulation was run long enough to reach steady-state. Then, at t = 400.1 s, a control rod unexpectedly moves. The reactivity increases from 0 to 0.15$ and the total power increases from 100 to 122% of the nominal power. By using an EKF time step of 0.4 s (for the P3D/R5 model, which represents the real NPP, time step is 0.2 s), the filter estimates the power, reactivity and moderator temperature based on one-directional neutron flux
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
12
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
Fig. 23. MATLAB EKF interface.
measurements. Based on the attack thresholds presented in Table 7, the CRA with the higher rod worth [1] is then automatically fully inserted: this action takes the system running at a safer and reduced level of 89% of the nominal power. The simulation allows the user, by mean of nuclear codes coupled with an I&C modelling and simulation tool (MATLAB), to assess the facility impact after an attack on a specific system function (control drive rod mechanism), and this information may eventually guide safety and security considerations regarding digital systems used for controlling important systems in NPP. 6. Conclusion The operation of current NPP and the deployment of new advanced reactors rely on digital I&C and IT systems. Digital I&C and IT systems are responsible for acquisition, transmission, analysis, delivery, and storage of data. By providing systems integration (actuators, sensors, displays, databases), digital systems improve performance and lead to reduced cost and better management. The pervasive use of digital systems motivates the development of new modeling and simulation tools to ensure minimum probability of reactor failure due to safety and cyber security threats. In this study, we implemented an innovative computational method to interconnect the nuclear-coupled P3D/R5 package and high-performance computational language MATLAB by developing a supervisory system. A P3D/R5 model simulated a 2772-MWt PWR B&W NPP during in both steady-state and transient, and MATLAB was used to drive the simulation by accessing (read and write) input, output, and restart files during code execution based on user-defined decision logic. This supervisory system aimed to extend the nuclear I&C modeling and simulation capabilities. Therefore, this implementation included:
a) Use of high level I&C MATLAB modeling capabilities to add flexibility to the simulation, as P3D/R5 has extremely limited supervision and control tools compared to MATLAB; b) Application of on-the-fly decision logic based on P3D/R5 simulation results (i.e., within each MATLAB time step, any I&C MATLAB decision rule has precedence over P3D/R5); c) Ability to assess separately the decision logic (MATLAB) from the NPP TH and NK points of view; d) Handling of input files by writing new information during code execution before calling P3D/R5 at each MATLAB time step; e) Use of a coding strategy that enhances the coupling features without changing the nuclear binary source codes. This is a critical feature once all risk is removed of invalidating results for designing, licensing, or accident analysis; and f) Saving exchanged data and output information for further analysis. In this study, the supervisory system performance assessment was based on a set of NPP simulations (code-to-code comparison) by accessing the relative nuclear power and control rod assembly reactivity. For simulation purposes, a direct CRA insertion from 411 (equivalent for the cases under simulation to 50% power) to 0 steps withdrawn and a reference set point from 50 to 80% and 50 to 0% were applied. The decision logic was based on a proportional integral controller. It is worth noting that different control methods (i.e., the decision logic) may produce different results from the PI implemented here. The main purpose of the PI controller and open-loop setup is to allow a benchmark among the codes to verify the supervisory system’s behavior. Because of MATLAB’s high versatility, the user can apply distinct decision logic based on manifold high-level computing tools. The supervisory system presented similar results when compared to the P3D/R5 package standalone runs. The direct
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
application of a decision rule by an open-loop approach gave us the same results as when editing the P3D input file and running the P3D/R5 package. This was only possible because we had access to post-simulation supervisory data, which were used to set the P3D/R5 package cards. Moreover, the application of the openloop CRA bank position dynamic (Cases c and d, see Table 1) using the values obtained from a closed-loop reactor power reference set point illustrated one of the many applications of the supervisory code. The flexibility of changing the MATLAB time step DTTR allowed us to meet the target computational times. In addition, a simplified practical application of the EKF algorithm developed in [1,3] and the supervisory system under study suggest that the EKF may be a powerful method for accurate detection and indication of abnormal plant operations due to cyberattacks. More specifically, the simulation presented here showed that the detection system was able to assess the facility impact (power change) after an attack on a system function (control drive rod mechanism). To summarize, the simulation of four distinct cases using open and closed-loop approaches showed that the strategy for coupling P3D/R5 and MATLAB was an improved testing and simulation method for different I&C scenarios. After the development of a strategy to access and write on-the-fly P3D/R5 cards (within each MATLAB time step) by a set of scripts and an s-function, the application of Simulink and MATLAB toolboxes was straightforward and added versatility to the simulation. The computational time was not a constraint when weighed against the benefits added by the environment technical language MATLAB. As a conclusion, this study suggested that the supervisory system may extend the nuclear capabilities of digital I&C modeling and simulations for safety and future cyber security assessment by incorporating the MATLAB flexibility into two globally used nuclear codes. Acknowledgments This work was supported by the Sao Paulo Research Foundation (FAPESP), grant 2016/046000; the Brazilian Research Council (CNPq), grant 302883/2018-5; and the International Atomic Energy Agency (IAEA), research project J02008, contract 22527.
13
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.anucene.2019.107098. References R. Busquim e Silva, Implications of advanced computational methods for reactivity initiated accidents in nuclear reactors. University of Sao Paulo, PhD Thesis, 2015. Available in: http://www.teses.usp.br/teses/disponiveis/3/3139/tde20072016-142605/pt-br.php. Kim, M.G., Lee, J.I., ]. Implication of LOCA characteristics of large PWR and SMR for future development of intelligent nuclear power plant control system. Ann. Nucl. Energy 127, 237–247. Busquim e Silva, R., Ferreira Marques, A.L., Cruz, J.J., Shirvan, K., Kazimi, M.S., ]. Reactivity estimation during a reactivity-initiated accident using the extended Kalman filter. Ann. Nucl. Energy 85, 753–762. Cheng, K., Meng, T., Zhao, F., Tan, S., ]. Development and validation of a thermal hydraulic transient analysis code for offshore floating nuclear reactor based on RELAP5/SCDAPSIM/MOD3.4. Ann. Nucl. Energy 127, 215–226. Janbazi, F., Rezaie, M.R., ]. ATWS severe accident analysis in the loss of flow scenario using the MELCOR code in Bushehr nuclear Power Plant. Kerntechnik 84, 123– 130. ed. 2. Tabadarl, Z. et al., ]. Simulation of a control rod ejection accident in a VVER-100/ V446 using RELAP5/Mod3.2. Ann. Nucl. Energy 45, 106–114. Hanselman, D., Littlefield, B., ]. Mastering MATLAB. Prentice-Hall, New Jersey. U.S. Nuclear Regulatory Commission (NRC). Computer Codes. Available in: https:// www.nrc.gov/about-nrc/regulatory/research/safetycodes.html. Accessed on May 8th, 2019. T. Kozlowski, D. Lee, D. Barber, C.H. Lee, J.Y. Cho, T. Downar, Purdue Advanced Reactor Core Simulator, PHYSOR 2002, pp. 1–7, Seoul, Korea, 2002. Yang, J., Yang, Y., Deng, C., Ishii, M., 0]. Best Estimate Plus Uncertainty analysis of a large break LOCA on Generation III reactor with RELAP5. Ann. Nucl. Energy 127, 326–340. T. Downar, Y. Xu, V. Seker, PARCS v2.7. U.S. NRC Core Neutronics Simulator, School of Nuclear Engineering, Purdue University, W. Lafayette, Indiana, 2006. MathWorks, Accelerating the Pace of Engineering and Science; Available in: https:// www.mathworks.com/videos/embracing-technical-computing-trends-withmatlab-accelerating-the-pace-of-engineering-and-science-90383.html, Accessed on May 12th, 2019. K.N. Ivanov et al. Pressurized Water Reactor Main Steam Line Break (MSLB) Benchmark. Volume I: Final Specifications. Organization for Economic CoOperation and Development (OECD), Nuclear Energy Agency (NEA) (1999). Ziegler, J.B., Nichols, N.B., 4]. Optimum settings for automatic controllers. ASME Trans. 64, 759–768.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
Contents lists available at ScienceDirect
Annals of Nuclear Energy journal homepage: www.elsevier.com/locate/anucene
Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment R.A. Busquim e Silva a,b,c, K. Shirvan b, J.J. Cruz a, R.P. Marques a, A.L.F. Marques a, J.R.C. Piqueira a,⇑ a
Polytechnic School of University of Sao Paulo, Av. Prof. Luciano Gualberto, 380, Sao Paulo, Brazil Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA, USA c Navy Technological Center in Sao Paulo, Av. Prof. Lineu Prestes, 2468 Sao Paulo, Brazil b
a r t i c l e
i n f o
Article history: Received 18 May 2019 Received in revised form 22 August 2019 Accepted 2 October 2019 Available online xxxx Keywords: NPP control PWR Digital instrumentation MATLAB PARCS RELAP
a b s t r a c t The operation of current nuclear power plants (NPP) and the deployment of advanced nuclear reactors rely on digital instrumentation and control (I&C) and information technology systems. The growing use of digital systems implies considerable security and safety challenges, as the usual standard codes used for designing, licensing, and analyzing accidents over the last few decades have been limited in terms of I&C features. This study presents an innovative method to couple PARCS/RELAP5 with the high-performance language MATLAB/Simulink. As a proof-of-concept proof, a 2772 MWt pressurized water reactor modeled using the coupled package PARCS/RELAP5 is simulated. A supervisory system is implemented that allows MATLAB/Simulink to oversee the nuclear codes for exchanging online data. In addition, I&C high-level computing tools are added to the NPP model, endorsing the application of ‘‘on-the-fly” decision logics. The coding strategy does not change the nuclear binary source codes, as the PARCS/RELAP5 codes have been validated for design and licensing. By running open- and closedloop cases for control rod assembly movement, we found that PARCS/RELAP5 capabilities for digital I&C modeling and simulations can be extended by adding MATLAB/Simulink flexibility. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction In the last four decades, computer nuclear codes have been extensively used by industry licensees and regulatory boards to simulate nuclear power plants (NPPs) and fuel cycle facility (FCF) systems. Simulations include design assessment and safety/ accident analysis. In addition, the nuclear industry foresees a massive increase in the use of digital instrumentation and control (I&C) systems, while also being wary of the emergence of new and sophisticated cyber threats. These factors have motivated the development of advanced coupled computational tools for improved management, cyber security assessment, higher process performance, and cost reduction [1,2]. The growing use of digital I&C systems raised one major security and safety issue related to the possibility that unauthorized codes may change the NPP behavior as a result of a cyber attack. The development of tools that allow the assessment of the impact on the facility, from the thermo-hydraulic (TH) and neutronic (NK) point of view, provides significant opportunities for conducting cyber defense experiments. By using MATLAB/Simulink, real ⇑ Corresponding author. E-mail address: [email protected] (J.R.C. Piqueira).
control equipment can be interfaced with the supervisory system to determine the consequences of sabotage by exploitation of security vulnerabilities that may result in loss of confidentiality, integrity and availability. In this study, an innovative computational method to interconnect two widely used nuclear codes is presented. The thermalhydraulic system code RELAP5 (R5) and the neutronic nuclear code PARCS 3D (P3D) are connected by using a well-known language for technical computing, namely, MATLAB/Simulink (MATLAB). The supervisory system that is developed extends the capabilities and features of the parallel virtual machine (PVM) package and P3D/ R5 I&C by coupling them with MATLAB. In addition, by accessing and managing the R5 and P3D input, output, and restart files during code execution, basic and advanced MATLAB toolboxes and nested functions may be applied to NPP models. The supervisory system was developed only for I&C modeling and simulation and not for actual I&C systems on NPP. For research purpose, the supervisory system was coupled to a 2772MWt pressurized water reactor (PWR) Babcock & Wilcox Co. (B&W) as modeled in [1,3]. Moreover, this simulation strategy enhances the coupling features without changing the nuclear binary source codes, which is a major innovation given that nuclear codes have been validated
https://doi.org/10.1016/j.anucene.2019.107098 0306-4549/Ó 2019 Elsevier Ltd. All rights reserved.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
2
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
during the last decades based on multiple benchmark activities for designing and licensing purposes. With the source code remaining unchanged, reliability of the simulation results is guaranteed. The implementation described herein includes: - Use of high-level I&C MATLAB modeling capabilities to add flexibility to the simulation, considering that P3D/R5 has limited supervision tools compared to MATLAB; - Application of on-the-fly decision logic based on P3D/R5 simulation results (i.e., within each MATLAB time step); - Ability to assess separately the decision logic (MATLAB) from the NPP behavior and from the TH and NK points of view; - Handling of the input files by writing new information before calling P3D/R5; and - Saving exchanged data and output information for further analysis. The next section presents some information about the simulation codes. After these preliminaries, two application sections describe how to use the developed tools; one section presents the supervisory system and the other the simulation strategies. Results and discussion conclude this work.
2. Simulation codes Since the 1970s, nuclear codes have been developed to simulate TH and NK NPP conditions. The analysis of NPP behavior usually begins with a TH system code that may include point kinetics (PK) or a more complex model to estimate the thermal power generation within the core [4]. The NK codes are mostly based on nodal expansion methods derived from neutron diffusion theory, where the macroscopic cross sections depend on feedback parameters such as fuel temperature, moderator density, and moderator temperature [5]. The use of nuclear TH and NK coupled codes allows a more complete NPP system analysis than when using standalone codes, which includes the feedback effect, core NK behavior, and TH dynamics. In addition, coupling tools allow more realistic simulation of normal operation (steady-state) and accident (transient) scenarios. This work uses the P3D [3] and R5 [6], which are advanced computer tools for TH and NK simulation, and the computing environmental software MATLAB [7] as a programming platform for implementing the supervisory system that drives the coupled simulations.
2.1. PARCS neutronic code The Purdue Advanced Reactor Core Simulator (PARCS) is a 3D reactor core simulator spatial kinetics code that solves the steady-state and time-dependent multi-group neutron diffusion and low-order transport equations in three-dimensional Cartesian geometry [8,9]. Furthermore, it has an embedded simplified TH solution, which accounts for the feedback effect when running in standalone mode, as PARCS uses an internal temperature calculation scheme based on input power. P3D is computer code that performs standalone calculations for light water reactors (LWRs), or it can be coupled to R5 using the PVM interface or directly to the TH system code TRACE-M [9]. P3D has been used for NK calculations by the United States Nuclear Regulatory Commission (NRC) [8]. Although P3D has few graphical user interfaces (GUIs), it works primarily by means of input and restart files for scenario configuration and output files for presenting results. Our simulations used PARCS version 2.71d.
2.2. RELAP5 system thermal-hydraulic code The Reactor Excursion and Leak Analysis Program (RELAP) is a tool for analyzing loss of coolant accidents (LOCA) and system transients [10]. R5 (i.e., version 5 of RELAP) has been used for licensing and regulation, evaluation of operator guidelines, and plant accident analysis. The code was developed to provide the firstorder effect needed to accurately predict system transients and to be sufficiently cost effective to enable sensitivity and parametric studies. Although it was initially designed for nuclear applications, the latest versions can be used for nuclear and non-nuclear applications involving vapor, liquid, non-condensable gases, and nonvolatile solutions [4]. R5 calculates the expected evolution of TH properties during steady state or transients, including the interaction among system components by using an internal temperature calculation scheme based on power input. It implements a two-fluid and sixequation model to simulate the two-phase coolant behavior and, consequently, the system’s TH behavior. However, the PK model inherent to R5 cannot consider multidimensional effects caused by spatial variations of feedback parameters. R5 also has few GUIs. However, it runs by means of input, output, and restart files. Our simulations used RELAP5/Mod3.2 version. 2.3. MATLAB computing environmental code The Matrix Laboratory (MATLAB) is a high-level numerical computing environmental software developed by The MathWorks Inc. for technical computing and programming not only for control engineering but also for data acquisition, data analysis, mathematical modeling, algorithm development, parallel computing, and web development. It has been in operation since 1984 [7]. MATLAB uses specific languages and codes controlled by its command window. MATLAB also allows interfacing with codes written in other languages including C, C++, Java, and FORTRAN. Hardware connections with field-programmable gate arrays or programmable logic controllers are available. It includes numerous toolboxes and packages such as Simulink, which adds graphical multi-domain simulation and model-based designs for dynamic systems, and MuPAD, which allows symbolic computing capabilities. 3. Supervisory system This work implements a supervisory system that couples P3D, R5 and MATLAB. This system oversees the nuclear codes online exchange of data; adds I&C high level computing tools within the NPP model; and endorses the application of ‘‘on the fly” decision logics. Therefore, the designed three codes coupling strategy is the heart of the supervisory system. The coupling strategy applied in this study extends the PVM coupling features without changing the P3D/R5 binary source codes. This enables the simulation results obtained using the developed supervisory code to be used for designing and licensing purposes. Through the access and management of R5 and P3D input, output, and restart files, basic and advanced MATLAB toolboxes and nested functions can be applied to NPP models specially developed for P3D and R5. The coupling strategy among the three codes includes: 1) The development of R5 and P3D standalone input and mapping files, which allows for proper TH and NK coupling; 2) The development of routines enabling the exchange of information between the coupled P3D/R5 and MATLAB in such a manner that MATLAB oversees the P3D/R5 exchange of data
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
within a predefined time step established for each simulation scenario; and 3) The development of functions to write and change R5 and P3D input files cards during code execution based on I&C decision rules.
3
- a control decision logic; - predefined R5 and P3D input file templates; and - a change of R5 and P3D input file cards during the simulation. Fig. 1 presents the time precedence among the codes. Basically, the simulation runs as follows:
3.1. P3D/R5 coupling The P3D and R5 codes work together to exchange information through a general interface [3], which manages the data flow between both codes by using the PVM. To guarantee the proper exchange of information between the codes, the NK-node and TH channel structures must match each other. This is accomplished by generating an input file named Maptab. The adequate generation of the Maptab file and correct determination of the mapping weight factors are necessary for obtaining accurate coupled code solutions. Moreover, the core nodalization in the radial and axial directions and their correlation have a considerable effect on determining the local core TH parameters and power distribution during transients. The mapping determines where to set NK node power in the TH channel (for both direct moderator heating and a heat structure fuel source) and where to assign TH properties in the NK node. Therefore, the coupling of the TH channels with the NK nodes in radial and axial directions depends on the adequate assignment of mapping weights between the two meshes. In this study, the mapping weight is based on the volumetric geometric fraction that correlates the TH and NK mesh in the radial and axial directions. The user is responsible for the assignments of the weighting factors. Then, a MATLAB script generates the Maptab file following the steps suggested in reference [11]. The TH conditions provided by the external systems code R5 at every time step in the NK spatial kinetics code P3D for the simulated case are the moderator temperature, vapor and liquid densities, void fraction, and averaged, centerline, and surface nuclear fuel temperatures [11]. These parameters are required at each NK node for the feedback calculation. The nodal/space dependent power distribution is then calculated by P3D and sent to R5, which calculates the heat conduction in the core heat structures. The TH code R5 addresses all the TH issues, and the NK code P3D spatial kinetics solves the NK challenges [1]. 3.2. P3D/R5/MATLAB coupling The PVM coupled P3D/R5 package and the MATLAB codes work together to exchange information through a leading-edge developed routine, implemented by means of a MATLAB system function (S-function) that manages the data flow between them. An Sfunction is a user coded dynamically linked subroutine that MATLAB can automatically load and execute. S-functions enable the code to work as built-in Simulink blocks [12]. The supervisory system uses the developed routine/S-function to stop/start and manage the input/output P3D/R5 files, and it uses the output from the same routine/s-function to implement the decision logic, which is applied to P3D/R5 closing the feedback loop during code execution. The supervisory system code allows the use of three nuclear codes in a single piece of software by means of: - an R5 and P3D spatial mesh that couples TH and NK using the PVM; - a predefined MAPTAB file that allows the correct matching between the NK node structure and TH channels (axially and radially); - a script that synchronizes the exchange of information; - on-the-fly access to the R5 and P3D restart and output files;
- The supervisory system code drives the steady-state run to generate the initial conditions for the transient simulations (restart files). Therefore, P3D/R5 runs within a Dt time step up to the predefined MATLAB DTSS; - During the transient operation, the supervisory system code drives the simulation based on a MATLAB DTTR time step that depends on the scenario being evaluated; - At each DTTR (DTTR Dt), the supervisory code calls P3D/R5 throughout the PVM and reads/uses the nuclear code restart information from the previous DTTR (or DTSS if it is the first transient supervisory time step); - The supervisory code determines the values and cards that must be written/changed in the input files based on the control decision logic; and - The supervisory code writes all changes in the input files and calls the nuclear P3D/R5 code again. The P3D/R5/MATLAB coupling strategy enables a more efficient and flexible I&C modeling and simulation for design, protection, and use of advanced control techniques. By using an S-function to access the coupled P3D/R5 package, the supervisory Simulink interface enables the incorporation of MATLAB algorithms and native toolboxes for instrumentation, control and simulation of dynamic systems. 3.3. NPP main modeling data The evaluation of the supervisory system performance is based on a set of NPP simulations and benchmarks. The NPP design is based on the 2772-MWt PWR B&W, which is well known from several studies published after the TMI Unit 1 accident. The main TMI NPP technical specifications can be found in the 1999 benchmark study released by the Nuclear Energy Agency from the Organization for Economic Cooperation and Development for mainstream line break [13] and references, such as [1] and [3]. Detailed information on P3D/R5 modeling can be found in reference [9,11]. The NPP core has a 17 17 assembly configuration. Sixty control rod assemblies (CRA) grouped into 7 banks (full-length control rods) are also employed. The core active height is 357.12 cm and
Fig. 1. Time step precedence among the codes.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
4
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
21.81cm reflectors on the top and bottom of the core are used, contributing to a total height of 400.74 cm. The position of the CRA insertion is given from the bottom of the lower reflector, where heights of 378.93 and 36.23 cm indicate completely withdrawn and completely inserted CRAs, respectively. The insertion is measured in steps of 0.35 cm (zero and 971 steps are equivalent to completely inserted and completely withdrawn CRAs, respectively). For NK modeling purposes, each radial node has the same flow area (FA) width of 21.81 cm, and each assembly is considered to be neutronically homogeneous. The FA has different U235 enrichments and burn-ups as well as a different number of burnable absorbers. The reactor is assumed to be in end of cycle (EOC), that is, 650 effective full power days (EFPD) or 24.58 GWd/MT. Axially, different compositions are used for each assembly, where each composition has a complete set of diffusion coefficients and macroscopic cross sections. The core is also divided axially in 28 NK regions. For TH purposes, the active core is axially divided into six fuel volumes and two reflectors. Radially, 18 TH channels and one regular channel represent the reflectors (channel 19).
4. Simulation strategy The main goal of supervisory system development is to enhance nuclear instrumentation and control modeling as well as simulation capabilities for future testing of digital I&C system reliability and redundancy for both safety and cyber security. For the purpose of this study, the performance of the supervisory system was assessed by assigning the control rod drive mechanism (CRDM) CRA positions by simulating the cases presented in Table 1. For these simulations, open- and closed-loop feedback supervisory setups were implemented. Moreover, the simulation results were benchmarked against P3D/R5 standalone runs. For comparison purposes, the nuclear code cards used the supervisory system post-simulation data. Prior to the steady-state runs, a set of P3D/R5 simulations were performed to identify the initial conditions of the control rod assembly that lead to criticality under steady-state conditions. For the same CRA axial positions (Banks 1 to 7), it was found that 50% power corresponds to 411 steps withdrawn and keff = 1.00083. In addition, 80% corresponds to 817 steps withdrawn and keff = 1.00079. During simulation, the supervisory code writes the adequate P3D cards into the input files, accesses the restart files, and manages the time step and output files. The R5 cards may also be changed on the fly depending on the simulation scenarios. Fig. 2 and Fig. 3 present the high-level Simulink block diagrams used for simulation. The plant block (blue) uses the S-function mechanism for extending the capabilities of the Simulink environment by means of compiled code dynamically linked to subroutines – all subroutines to access the P3D/R5 files are called under the Plant block. The Decision Logic diagram (green) contains decision logic for each one of the cases under study (e.g. the PI controller). The block diagrams in orange allow the use of data within the MATLAB workspace (e.g. for plotting).
Table 1 Cases under investigation. Case
Case Description
A B C D
CRA total insertion from 411 to 0 steps Power set point from 50 to 80% relative power CRA withdrawn from 411 to 817 steps Power set point from 50% to hot zero power (HZP)
4.1. Open-loop set up In the open-loop setup as presented in Fig. 2, the output of the plant (i.e. the relative reactor power, that is, variable ‘‘rktpow” in the plant block) has no effect on the input of the plant (i.e., the CRA position, that is, variable ‘‘bank pos cmd” in the Decision Logic block) that is applied to the core CRDM. No feedback information is used to adjust or control the reactor power. For simulation purposes, all CRA were assumed to be in the same axial position at the beginning of the simulation. 4.2. Closed-loop feedback setup Closed-loop systems compare the output information to the desired set point (i.e., the reference) and generate an error. This error is then used as input for a controller designed to have the output more closely from the reference. The design and implementation of the controller depends on the system being controlled. For the closed-loop feedback set up, Fig. 3 no CRA axial position is set in advance. During the simulation, the supervisory system measures the nuclear power ‘‘rktpow” (process variable) at each supervisory DTTR time step and sends its value to the Decision Logic block. The Decision Logic block receives the reference ‘‘rktpow_ref” (which may be set by the reactor operator, for instance), compares both values, and generates an error. Then, the controller calculates the CRA bank position ‘‘bank_pos_cmd” to be applied to the core CDRM. The controller implemented here was a proportional-integral (PI) type. A PI controller is a control loop mechanism widely used in the industry for minimizing the error e(t) over time between a desired set point and a measured process variable. It is named PI controller due to a proportional term (P ¼ K p eðt Þ) and and integral term Rt (I ¼ K i 0 eðsÞ dt). The P term takes into account a gain factor Kp that is proportional to the difference between the desired set point and the measured variable. The I term accounts for accumulated past values of the error. By integrating the residual error between the desired set point and the measured variable, it minimizes its residual due to cumulative values. A PI controller is commonly used in various practical applications even though it has many limitations, including an inability to provide optimal control (it may be tuned without direct knowledge of the process) and the use of constant controller parameters. However, it can be easily applied by knowing the dynamics of the plant. The PI input is the error e(t) and the output is the control signal u(t) applied to the plant according to the following equation:
Z uðt Þ ¼ K p eðtÞ þ K i 0
t
eðsÞ dt
Here, Kp allows the controller to make a change to the output proportional to the error. An improper tuning of the proportional term may lead to an unstable (too high) or slow sensitive (too small) response. In addition, Ki is the integral term that mitigates the steady state error as derived from the proportional term as it simultaneously speeds up the output response toward the set point. The integral parameter was set assuming an integral time Ti to minimize the output/set point offset and making K i ¼ 1=T i , and Kp was set to avoid an unstable output response. For all cases under investigation, Kp was taken as 5 and Ti was assumed to be 6. The proportional and derivative terms were determined by applying the Ziegler-Nichols method [14] followed by manual fine tuning. (It is worth noting that the PI implementation works in a continuous time framework, and the plant output is taken at each DTTR).
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
5
Fig. 2. MATLAB/Simulink open-loop interface.
Moreover, the NPP ran in steady state mode for 1210 s (all plots started in t = 1200 s for better viewing). Then, from 1210 to 1380 s, it ran in transient conditions. During the steady state, the coupled TH and NK solution converged within the 1200 s such that all its values remained almost constant. During the transient mode, the supervisory system operated as shown in Fig. 1. Moreover, the MATLAB DTTR time step was set to 1 s for all cases. P3D/R5 ran on a 103-s Dt time step. For comparison, a closed-loop simulation with DTTR = 0.2 s was also performed.
5. Simulation results and discussion The supervisory system adequately managed to access and change cards in the P3D/R5 input files based on the decision logic output. It also accessed data in the output and restart files, as shown Fig. 1. However, the supervisory system was able to drive the reactor dynamics. Therefore, its results were very similar to the P3D/R5 standalone runs. This supervisory system was designed to apply to more complex on-the-fly I&C tools as well as to offer a
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
6
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
Fig. 3. MATLAB closed-loop interface.
cyber security strategy for future applications by using MATLAB real digital system connectivity. While the cases under investigation were simple from I&C point-of-view, the model details are realistic enough to verify the computational methodology. For simulation purposes, the MATLAB DTTR time step was set at 1 s. Thus, the decision logic action was delayed by 1 s when compared to the P3D/R5 standalone runs. DTTR was defined by the user at the beginning of the simulation. The nuclear power behavior is explained by the CRA insertions/withdraws and the correspondent CRA reactivity trend. Table 4 summarizes the simulation results by presenting the standard deviation (SD) and mean between the P3D/R5 and supervisory reactor power.
5.1. CRA total insertion from 411 to 0 steps (Case A) The open-loop (see Fig. 2) axial CRA insertion from 411 to 0 steps withdrawn as shown in Fig. 4 was set in the Simulink decision logic block in which no feedback information was given. The power response after the total CRA insertion is shown in Fig. 5. As expected, the nuclear reactor power decreased as the control rod assemblies were inserted, and the CRA reactivity was negative (see Fig. 6). It can be seen that the reactor power remained under shutdown conditions after 1,360 s, while the CRA reactivity was expected to return gradually to zero.
Fig. 4. CRA bank position: Case A.
In addition, Fig. 7 and Fig. 8 show that the 3D core power distribution decreased as the CRA was inserted, which was expected, and it remained relatively uniform during the transient simulation.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
7
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
Fig. 5. Total power: Case A. Fig. 8. 3D power distribution t = 1320 s.
Fig. 6. CRA reactivity: Case A.
This behavior was similar for all cases under investigation. The 3D power distribution represents the assembly peaking factor (i.e., the assembly power divided by the average assembly power in the core, which is the reactor total power divided by the number of fuel assemblies). In addition, to compare results, the same pre-defined axial position was set in the coupled P3D/R5 package ‘‘move_bank” card. However, a modeling limitation existed in the card when the standalone simulation was run: the pair (time, bank position) could be repeated up to 10 times. However, the supervisory system runs P3D/R5 and writes new ‘‘move_bank” values within the MATLAB DTTR. Thus, the supervisory system can move the CRA 10 times per time step, which enables better control. Table 2 shows the bank position assumed for the P3D/R5 benchmark (up to 1200 s, whereby all CRA are 411 steps withdrawn and, after 1320 s, remain fully inserted). Fig. 9 and Fig. 10 show that no significant difference was found between the open-loop supervisory and P3D/R5 package benchmark runs. Fig. 11 and Fig. 12 present the power and CRA reactivity differences between the P3D/R5 standalone and supervisory runs. It can be seen that the differences were very small. The main differences were derived from the 1 s delay between the results and the user-defined MATLAB DTTR time step.
5.2. Power set point from 50 to 80% relative power (Case B) This case was evaluated by a closed-loop setup (see Fig. 3) assuming a 50 to 80% relative power set point. Fig. 13 shows that the application of the PI controller led, as expected, to an increase in nuclear power toward the reference. The control signal calculated by the PI controller, which was applied by the supervisory system to the CRDM, is shown in Fig. 14.
Table 2 CRA position: Case A.
Fig. 7. 3D power distribution t = 1210 s.
Time (s)
CRA Position (steps)
Time (s)
CRA Position (steps)
1200 1210 1220 1240 1250
411 411 372 300 261
1260 1280 1300 1320 1380
224 148 73 0 0
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
8
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
Fig. 9. Total power: Case A and P3D/R5.
Fig. 12. CRA reactivity difference.
Fig. 10. CRA reactivity: Case A and P3D/R5.
Fig. 13. Total power: Case B.
Fig. 11. Power difference.
Fig. 14. CRA position: Case B.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
5.3. CRA withdrawn from 411 to 817 steps (Case C) For Case C, the CRA bank position obtained (calculated) from Case B (see Fig. 14) was applied directly to the CRDM using the open-loop setup. The expected nuclear power trend followed a ramp from 50 to 80% nominal power. In addition, a P3D/R5 package standalone case was simulated by setting the P3D ‘‘move_bank” card according to Table 3. (This was a post-simulation data setup obtained from Fig. 14 by direct observation.) Then, Cases B and C as well as the P3D/R5 standalone run and reference were all compared. Fig. 15 and Fig. 16 show that, in addition to the differences derived from the MATLAB DTTR and Table 3 CRA position: Case C. Time (s)
CRA Position (steps)
Time (s)
CRA Position (steps)
1200 1210 1220 1240 1260
411 411 424 467 525
1280 1300 1320 1380 1400
590 660 732 806 817
9
P3D/R5 time steps as well as the discrete P3D/R5 package application of the ‘‘move_bank” card, all simulation results were quite similar. 5.4. Power set point from 50% to HZP (Case D) This case was also evaluated by a closed-loop reference from 50% power to HZP setup. Then, the axial position (postsimulation data) from the supervisory results was set in the P3D/ R5 package ‘‘move_bank” card, and a standalone simulation was performed. The results are presented in Fig. 17 and Fig. 18. As it can be seen, the results were very similar. The differences as presented in Fig. 19 and Fig. 20 were due to the supervisory system time step and P3D ‘‘bank_move” card that was limited in terms of the number of pairs (time, bank position). 5.5. Standard deviation, mean, maximum and minimum benchmark data To benchmark supervisory results against the P3D/R5 package standalone runs, Table 4 presents the SD, mean, and maximum and minimum values calculated over the difference (DP) between
Fig. 17. Total power: Case D. Fig. 15. Total Power: Case B.
Fig. 16. CRA reactivity: Case B and P3D/R5.
Fig. 18. CRA reactivity: Case D and P3D/R5.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
10
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
Fig. 19. Power difference.
Fig. 21. Total power: Case B.
Fig. 20. CRA difference. Fig. 22. CRA bank position: Case B.
Table 4 DP simulation data. Case
SD (%)
Mean (%)
Max (%)
Min (%)
A B C D
0.301 0.302 0.290 0.308
0.226 0.211 0.175 0.250
0.0 1.036 0.398 0.324
1.123 0.671 0.001 1.617
the reactor power obtained with the P3D/R5 package (open- and closed-loop) and supervisory system. It can be seen that the SD was close to 0.3% for all cases under simulation and that the mean remained close to zero (in the range 0:25 DP 0:2). These were acceptable results considering the dynamics of the plant. The maximum and minimum values were also suitable results considering the delay due to the MATLAB DTTR time step. 5.6. Simulation results for distinct MATLAB time step The application of smaller values of MATLAB DTTR may lead to better results. However, the computational time will increase.
For simulation purposes, a DTTR = 0.2 s discrete PI (parameters Kp = 3 and Ti = 5) was implemented and both results were compared. Fig. 21 and Fig. 22 show that no significant difference was found between the results for a total power increase from 50 to 80%. However, it must be noted that the appropriate MATLAB DTTR must be chosen based on the expected dynamics of the plant under the simulations performed. 5.7. Computation running time The coupled package P3D/R5 is a time-consuming computational code. Nevertheless, most computational time is consumed by the neutronic and thermo hydraulic calculation. The supervisory system coupling task includes the opening of different files for binary read access as well as formatting and writing data, which are also time-consuming. In addition, it uses feedback information as input for the decision logic. For further improvements to the supervisory system regarding its computational demand, the elapsed time can be seen in Table 5 and Table 6 for the cases under investigation (all runs were performed on a Dell Laptop with an Intel Core i7 7500U CPU @ 2.7 GHz running Windows 10 at 64 bits). The simulation time
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
11
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx Table 5 Time step DTTR = 1 s elapsed time (s). Case
Supervisory System (P3D/R5/MATLAB)
Nuclear Power from 50 to 0% Nuclear Power from 50 to 80%
Open-loop
Closed-loop
1.713 103 1.762 103
1.740 103 1.770 103
P3D/R5 Standalone
Ratio Supervisory System (P3D/R5/MATLAB) Open-loop to P3D/R5
Ratio Supervisory System (P3D/R5/MATLAB) Closed-loop to P3D/R5
1.490 103 1.491 103
1.15 1.18
1.17 1.19
Table 6 Time step DTTR = 0.2 s elapsed time (s). Case
Supervisory System (P3D/R5/MATLAB) Closed-loop
P3D/R5 Standalone
Ratio Supervisory System (P3D/R5/MATLAB) to P3D/R5
Nuclear Power from 50 to 0% Nuclear Power from 50 to 80%
3.081 x104 3.086 x104
1.490 103 1.491 103
20.68 20.70
Table 7 EKF cyber attack scenario. Core Conditions
CRA at t = 0
Attack Threshold
Cyber-Attack
Steady-State 100% Power
Fully/Partially withdraw
DPower>±10% nominal power D reactivity >±0.1 $ D moderator temperature >±10 °C
CRA goes from fully inserted to fully withdraw HMI compromised
was approximately 1.15 to 1.19 times higher for the supervisory system at DTTR = 1 s, and there was a greater computational demand of approximately 20.7 times for the smaller MATLAB time step of 0.2 s. This was expected, as the supervisory executed a greater number of open/write actions for a small time step. (For data consistency, the elapsed time did not consider the time when the simulation was not running or paused.) No significant difference was found between open- and closed-loop DTTR = 1 s simulations. The running time depended on the MATLAB DTTR time step, the number of restart block files desired, and the dynamics of the simulation. The additional time as a result of the MATLAB coupling may be reduced by retaining information in the CPU memory and by reducing the number of R5 restart files to one for each MATLAB DTTR. Although this was not the goal of this study, further implementation (see reference [1]) intends to minimize the computational time. 5.8. EKF as cyber security point defense option: a simplified application The use of nuclear codes, validated with well-established criteria and methods, provide nuclear safety analysis advantages over conventional non nuclear codes used for I&C modeling and simulation. On the other hand, specially designed I&C modeling and simulation software may add flexibility to potential users that may not invest resources to qualify new software for their reactor control system. Moreover, the use of the supervisory system provides significant opportunities for conducting cyber defense investigations due to the capability of assessing both digital system function and facility function impacts after a cyber-attack. This section presents a practical application of the stochastic-forecasting extended Kalman filter (EKF) algorithm developed in reference [1,3] integrated with the supervisory system, that provides a point defense option against cyberattacks. The EKF model is based on the point kinetics reactor equations, on the temperature feedback effects and on simplified core-coolant mass and energy balance equations. The EKF is an extension of the
Kalman Filter (KF) for nonlinear systems, which works on a stream of noisy input data to generate a statistically optimal estimate of the variable of interest. The KF is optimal linear estimator that minimizes the mean square error of the estimated variable if the noise has a Gaussian distribution. The EKF estimation is based entirely on the one-directional neutron flux (or reactor power). For simulation purposes, noise is added to the measured signal. The EKF algorithm produces estimates of the reactor power, moderator temperature and reactivity and compares these estimated values to their measured counterparts. Based on the detection criteria presented in Table 7, a control signal is sent to the NPP protection system. The scenario analyzed under this study considers: - The control rod drive mechanism is compromised during the NPP normal operation (100% power). - The control rods are take to ‘‘all rod out” position increasing reactivity. - The human-machine interface (HMI) is compromised by hiding the aforementioned action. - The total power, reactivity and moderator temperature HMI indications are not updated. - The EKF, using an one-directional/direct neutron flux signal, estimates the reactivity and the moderator temperature. - Detection criteria, designed to minimize false positives and missing alarms; - The detection criteria assume low and high thresholds for detection purposes. Fig. 23 presents the EKF implementation. The Decision Logic block implements the attack threshold shown in Table 7. The simulation was run long enough to reach steady-state. Then, at t = 400.1 s, a control rod unexpectedly moves. The reactivity increases from 0 to 0.15$ and the total power increases from 100 to 122% of the nominal power. By using an EKF time step of 0.4 s (for the P3D/R5 model, which represents the real NPP, time step is 0.2 s), the filter estimates the power, reactivity and moderator temperature based on one-directional neutron flux
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
12
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
Fig. 23. MATLAB EKF interface.
measurements. Based on the attack thresholds presented in Table 7, the CRA with the higher rod worth [1] is then automatically fully inserted: this action takes the system running at a safer and reduced level of 89% of the nominal power. The simulation allows the user, by mean of nuclear codes coupled with an I&C modelling and simulation tool (MATLAB), to assess the facility impact after an attack on a specific system function (control drive rod mechanism), and this information may eventually guide safety and security considerations regarding digital systems used for controlling important systems in NPP. 6. Conclusion The operation of current NPP and the deployment of new advanced reactors rely on digital I&C and IT systems. Digital I&C and IT systems are responsible for acquisition, transmission, analysis, delivery, and storage of data. By providing systems integration (actuators, sensors, displays, databases), digital systems improve performance and lead to reduced cost and better management. The pervasive use of digital systems motivates the development of new modeling and simulation tools to ensure minimum probability of reactor failure due to safety and cyber security threats. In this study, we implemented an innovative computational method to interconnect the nuclear-coupled P3D/R5 package and high-performance computational language MATLAB by developing a supervisory system. A P3D/R5 model simulated a 2772-MWt PWR B&W NPP during in both steady-state and transient, and MATLAB was used to drive the simulation by accessing (read and write) input, output, and restart files during code execution based on user-defined decision logic. This supervisory system aimed to extend the nuclear I&C modeling and simulation capabilities. Therefore, this implementation included:
a) Use of high level I&C MATLAB modeling capabilities to add flexibility to the simulation, as P3D/R5 has extremely limited supervision and control tools compared to MATLAB; b) Application of on-the-fly decision logic based on P3D/R5 simulation results (i.e., within each MATLAB time step, any I&C MATLAB decision rule has precedence over P3D/R5); c) Ability to assess separately the decision logic (MATLAB) from the NPP TH and NK points of view; d) Handling of input files by writing new information during code execution before calling P3D/R5 at each MATLAB time step; e) Use of a coding strategy that enhances the coupling features without changing the nuclear binary source codes. This is a critical feature once all risk is removed of invalidating results for designing, licensing, or accident analysis; and f) Saving exchanged data and output information for further analysis. In this study, the supervisory system performance assessment was based on a set of NPP simulations (code-to-code comparison) by accessing the relative nuclear power and control rod assembly reactivity. For simulation purposes, a direct CRA insertion from 411 (equivalent for the cases under simulation to 50% power) to 0 steps withdrawn and a reference set point from 50 to 80% and 50 to 0% were applied. The decision logic was based on a proportional integral controller. It is worth noting that different control methods (i.e., the decision logic) may produce different results from the PI implemented here. The main purpose of the PI controller and open-loop setup is to allow a benchmark among the codes to verify the supervisory system’s behavior. Because of MATLAB’s high versatility, the user can apply distinct decision logic based on manifold high-level computing tools. The supervisory system presented similar results when compared to the P3D/R5 package standalone runs. The direct
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098
R.A. Busquim e Silva et al. / Annals of Nuclear Energy xxx (xxxx) xxx
application of a decision rule by an open-loop approach gave us the same results as when editing the P3D input file and running the P3D/R5 package. This was only possible because we had access to post-simulation supervisory data, which were used to set the P3D/R5 package cards. Moreover, the application of the openloop CRA bank position dynamic (Cases c and d, see Table 1) using the values obtained from a closed-loop reactor power reference set point illustrated one of the many applications of the supervisory code. The flexibility of changing the MATLAB time step DTTR allowed us to meet the target computational times. In addition, a simplified practical application of the EKF algorithm developed in [1,3] and the supervisory system under study suggest that the EKF may be a powerful method for accurate detection and indication of abnormal plant operations due to cyberattacks. More specifically, the simulation presented here showed that the detection system was able to assess the facility impact (power change) after an attack on a system function (control drive rod mechanism). To summarize, the simulation of four distinct cases using open and closed-loop approaches showed that the strategy for coupling P3D/R5 and MATLAB was an improved testing and simulation method for different I&C scenarios. After the development of a strategy to access and write on-the-fly P3D/R5 cards (within each MATLAB time step) by a set of scripts and an s-function, the application of Simulink and MATLAB toolboxes was straightforward and added versatility to the simulation. The computational time was not a constraint when weighed against the benefits added by the environment technical language MATLAB. As a conclusion, this study suggested that the supervisory system may extend the nuclear capabilities of digital I&C modeling and simulations for safety and future cyber security assessment by incorporating the MATLAB flexibility into two globally used nuclear codes. Acknowledgments This work was supported by the Sao Paulo Research Foundation (FAPESP), grant 2016/046000; the Brazilian Research Council (CNPq), grant 302883/2018-5; and the International Atomic Energy Agency (IAEA), research project J02008, contract 22527.
13
Appendix A. Supplementary data Supplementary data to this article can be found online at https://doi.org/10.1016/j.anucene.2019.107098. References R. Busquim e Silva, Implications of advanced computational methods for reactivity initiated accidents in nuclear reactors. University of Sao Paulo, PhD Thesis, 2015. Available in: http://www.teses.usp.br/teses/disponiveis/3/3139/tde20072016-142605/pt-br.php. Kim, M.G., Lee, J.I., ]. Implication of LOCA characteristics of large PWR and SMR for future development of intelligent nuclear power plant control system. Ann. Nucl. Energy 127, 237–247. Busquim e Silva, R., Ferreira Marques, A.L., Cruz, J.J., Shirvan, K., Kazimi, M.S., ]. Reactivity estimation during a reactivity-initiated accident using the extended Kalman filter. Ann. Nucl. Energy 85, 753–762. Cheng, K., Meng, T., Zhao, F., Tan, S., ]. Development and validation of a thermal hydraulic transient analysis code for offshore floating nuclear reactor based on RELAP5/SCDAPSIM/MOD3.4. Ann. Nucl. Energy 127, 215–226. Janbazi, F., Rezaie, M.R., ]. ATWS severe accident analysis in the loss of flow scenario using the MELCOR code in Bushehr nuclear Power Plant. Kerntechnik 84, 123– 130. ed. 2. Tabadarl, Z. et al., ]. Simulation of a control rod ejection accident in a VVER-100/ V446 using RELAP5/Mod3.2. Ann. Nucl. Energy 45, 106–114. Hanselman, D., Littlefield, B., ]. Mastering MATLAB. Prentice-Hall, New Jersey. U.S. Nuclear Regulatory Commission (NRC). Computer Codes. Available in: https:// www.nrc.gov/about-nrc/regulatory/research/safetycodes.html. Accessed on May 8th, 2019. T. Kozlowski, D. Lee, D. Barber, C.H. Lee, J.Y. Cho, T. Downar, Purdue Advanced Reactor Core Simulator, PHYSOR 2002, pp. 1–7, Seoul, Korea, 2002. Yang, J., Yang, Y., Deng, C., Ishii, M., 0]. Best Estimate Plus Uncertainty analysis of a large break LOCA on Generation III reactor with RELAP5. Ann. Nucl. Energy 127, 326–340. T. Downar, Y. Xu, V. Seker, PARCS v2.7. U.S. NRC Core Neutronics Simulator, School of Nuclear Engineering, Purdue University, W. Lafayette, Indiana, 2006. MathWorks, Accelerating the Pace of Engineering and Science; Available in: https:// www.mathworks.com/videos/embracing-technical-computing-trends-withmatlab-accelerating-the-pace-of-engineering-and-science-90383.html, Accessed on May 12th, 2019. K.N. Ivanov et al. Pressurized Water Reactor Main Steam Line Break (MSLB) Benchmark. Volume I: Final Specifications. Organization for Economic CoOperation and Development (OECD), Nuclear Energy Agency (NEA) (1999). Ziegler, J.B., Nichols, N.B., 4]. Optimum settings for automatic controllers. ASME Trans. 64, 759–768.
Please cite this article as: R. A. Busquim e Silva, K. Shirvan, J. J. Cruz et al., Advanced method for neutronics and system code coupling RELAP, PARCS, and MATLAB for instrumentation and control assessment, Annals of Nuclear Energy, https://doi.org/10.1016/j.anucene.2019.107098