Development and Testing of TRACE/PARCS ECI Capability for Modelling CANDU Reactors with Reactor Regulating System Response

The use of the USNRC codes TRACE and PARCS has been considered for the coupled safety analysis of CANDU reactors. A key element of CANDU simulations is theinteractions betweenthermal-hydraulic and physic phenomena with theCANDU reactor regulating system (RRS).Todate,noorlimiteddevelopmenthas takenplaceinTRACE-PARCSin thisarea.Inthiswork,thesystemthermal-hydrauliccode TRACE_Mac1.0 is natively coupled with the core physic code PARCS_Mac1.0, and RRS control is implemented via the exterior communications interface (ECI) in TRACE. ECI is used for coupling the external codes to TRACE, including additional physical models and control system models. In this work, a Python interface to the TRACE ECI library is developed, along with an RRS model written in Python. This coupling was tested using a CANDU-6 IAEA code coupling benchmark and a 900MW CANDU model for various transients. For the CANDU-6 benchmark, the transients did not include RRS response, however, the TRACE_Mac1.0/PARCS_Mac1.0 coupling and ECI script functionality was compared to the previous benchmark simulations, which utilized external coupling. For the 900MW CANDU simulations, all aspects of the ECI module and RRS were included. The results from the CANDU-6 benchmark when using the built-in coupling are comparable to those previously achieved using external coupling between the two codes with coupled simulations taking 2 x to 3 x less execution time. The 900MW CANDU simulations successfully demonstrate the RRS functionality for the loss of ﬂow events, and the coupled solutions demonstrate adequate performance for ﬁgure-of-eight ﬂow instability modeling.


Introduction
In nuclear power plants, including CANDU stations, many phenomena arise because of coupled interactions between the reactor physic (nuclear chain reaction) and system thermalhydraulic phenomena. ese range from simple reactivity feedbacks to 3D spatial power changes to coupled flow instabilities. For this reason, many core physic and system thermalhydraulic codes include coupling capabilities. is study focuses on the USNRC-developed codes, PARCS (for core physics), and TRACE (for system thermal-hydraulics) [1,2]. While these codes are designed for the safety analysis of light water reactors, they have capabilities for the analysis of other reactor types, including CANDU reactors. For example, TRACE and PARCS were used to model the IAEA ICSP benchmark problem, including an uncertainty analysis [3].
In addition to coupled phenomena, it is also important to consider the effect of the reactor's control systems on transient reactor response. In TRACE, these are typically modeled using signal variables and control blocks. Signal variables can read process parameters, while control blocks perform analog and logical calculations. Both can be used to actuate devices, such as valves and pumps. e coupled TRACE-PARCS code then allows the TRACE control the blocks to alter reactivity device configurations in PARCS and mimic the response of a reactor control system. is reactivity device coupling was introduced in PARCS 2.7 [4] and used in a previous study for pressurized water reactor (PWR) analysis [5,6] but is not officially documented in the code's user manuals.
Currently, most work in CANDU safety analysis is performed using codes developed specifically for CANDU reactors. While these codes operate under the same principles as more general-purpose safety analysis codes, they include models and correlations that are better suited for CANDU reactors specifically. e examples of codes developed and used for CANDU design and analysis include RFSP for core physics and reactor regulation [7], along with CATHENA and TUF for system thermal-hydraulics [8,9]. e Canadian Nuclear Safety Commission (CNSC) has been evaluating the use of various independent codes for CANDU safety analysis, particularly TRACE, motivating various studies on its applicability [10], including this study.
Most previous studies involving the implementation of CANDU-specific control systems embed the reactor regulating system (RRS) capabilities into the source code. One study of a CANDU loss of flow transient using RELAP5 added station-specific controller routines to the source code [11]. Alternatively, RRS emulators can be developed within an outside environment and included as part of external coupling methodologies. It was done in a prior study using TRACE and PARCS to simulate multiple CANDU transients [3]. In such an arrangement, an external script is developed, which executes each coupled code independently over short periods of time and exchanges information at the end of each time step. During this information exchange, the response of the reactor regulating system can be determined and used to alter the associated input files to include control device changes prior to initiating the next time step. While the former methodology allows for fast execution times and avoids the I/O bottlenecks in the second method, it is inflexible compared to an external RRS emulator. Ideally, a solution involving the close coupling of the thermal-hydraulic and physic phenomena with an external RRS emulator would be the ideal solution. It would provide the flexibility of RRS development within its own shell while still maintaining the computational advantages. Using TRACE-PARCS with the exterior communications interface (ECI) module provides one avenue to achieve these goals.
TRACE includes a general coupling capability with the exterior communications interface, which allows TRACE to be coupled with any code that uses the ECI library [12]. It can be used for additional physical models (e.g., fuel performance models, subchannel models, CFD models) and for detailed control systems that would be difficult to implement using only TRACE control blocks. is study focuses on the implementation of RRS response via the ECI interface.
is study develops a model for much of the CANDU RRS and tests the coupling framework and RRS model on two existing coupled PARCS-TRACE models (one CANDU-6 and one 900 MW class CANDU).
is work converts the models to utilize the built-in PARCS-TRACE coupling and integrate the other systems using ECI capability with the goal of streamlining the execution of the models while still obtaining accurate results.
While one of the goals of ECI is to avoid the source code changes to TRACE, some changes were nevertheless required. is work utilizes the modified versions of TRACE V5.1262 and PARCS V3.31, hereafter referred to as TRACE_Mac1.0 and PARCS_Mac1.0, respectively, wherever the results and discussion are specific to the modified codes.
e changes are relatively minor, primarily to facilitate code coupling for this work's applications, and they do not add or modify any physical models.

PARCS-TRACE Coupling Background.
e PARCS core physic code is included within the distribution of TRACE. e two codes are compiled into a single executable. e PARCS and TRACE data structures are not linked directly but are linked by an internal "general interface." PARCS and TRACE access shared data transfer tables rather than directly accessing each other's data structures. Upon initialization, PARCS sets up the data transfer tables based on the data in the MAPTAB file. e MAPTAB file is a usercreated file that specifies how the nodes of a PARCS model are to be coupled with TRACE components [13].
TRACE and PARCS are coupled explicitly as the two solvers run independently. At the end of each TRACE time step, thermal-hydraulic data is passed to PARCS, and a step is performed in PARCS-converting the flux in steady state mode or advancing time in transient mode. e updated core power data is passed back to TRACE for the next time step. e accuracy of this coupling is equivalent to externally controlled coupling for a given data transfer frequency.
is accuracy has been evaluated in a previous study using the external coupling methodology [3]. However, there is a performance advantage in performing the data transfers in-memory, keeping both codes initialized in memory, and running in lockstep, making it practical to perform data transfer on every TRACE time step. In addition, PARCS has the capability to skip TRACE time steps based on user input as oftentimes it is unnecessary to update the neutronic model on every thermalhydraulic time step.
To model CANDU reactivity devices, this work utilized an undocumented ability to couple PARCS control rod banks to TRACE signal variables. Rather than physical control rods in a light water reactor (LWR), the models in this work mimic the liquid zone devices and adjusters of a CANDU [14]. It is done using the %CRSIG card [4]  is research focuses on method 3, which is the only method that can dynamically react to the rest of the system in a generalized fashion.

Exterior Communication Interface (ECI) Overview.
e exterior communications interface is a library of Fortran and C subroutines that permits the coupling of other codes to TRACE without modifying the TRACE source code or executable [12]. It is also capable of coupling multiple TRACE submodels together to parallelize the simulation of a larger model. Using ECI, through TRACE, allows one to directly access and control reactivity device configurations within PARCS. e two main components of the ECI are as follows [12]: (1) e ECI library, which is embedded within TRACE.
A separate copy of the library is included so that programmers may embed it in their own programs. (2) e ECI driver, which is a standalone Java program that must run in the background, is responsible for starting child "satellite" processes and setting up the actual interprocess communication through sockets or shared memory. e ECI coupling model is "request-driven" [12]. ECI defines 18 synchronization points at which the parallel processes exchange data, corresponding to different points in the TRACE program flow, as shown in Figure 1. Each program specifies its requests during the initialization of the simulation. Each request identifies the variable to be coupled, along with a synchronization point and the direction of data transfer. Multiple requests may exist for the same variable. ECI then locates all of the requested variables and constructs the transfer tables. During execution, the central process (usually TRACE) is responsible for time step control and status monitoring, while satellite processes can request a smaller time-step size or report on their convergence status.
PARCS itself does not have ECI support. us, ECI programs can only interact with PARCS indirectly through TRACE. It is done by transferring data from signal variables, fluid components, or heat structures in TRACE. It limits the manipulation of PARCS to control banks by the manipulation of TRACE control blocks and signal variables. In Figure 1, PARCS initialization and PARCS execution occur after Input and before EndStep, respectively. e ECI library, distributed with TRACE, is the set of codes required to create an ECI-compatible program. It is written primarily in Fortran 90 and designed to work with Fortran 90 programs by including the library in the program's source code (rather than as a dynamic or static library). While most of the modules can be used as-is, two modules are templates that must be completed by the programmer, SpecExTrans and TimeEvolve [12]. e former allows the program to locate variables requested by other processes, while the latter contains the required program flow and synchronization points.
An ECI satellite program will typically retrieve initial conditions at Init and old-time values at OldTime, as shown in Figure 1. Values can be sent and new-time values retrieved wherever appropriate.
A complete ECI program, therefore, requires the following: (1) e ECI library distributed with TRACE (2) A completed version of the SpecExTrans module if its variables are to be visible to other processes (3) e TimeEvolve subroutine that includes the synchronization points and program flow from Figure 1 and all of the program's own computations (4) A subroutine that sets up all of the data requests that the program requires (5) A "main" function or a subroutine that performs all of the necessary ECI subroutines, including reading command-line arguments, preparing the data transfers, and calling TimeEvolve Specific implementation details are omitted for brevity. ey may be found in the ECI manual [12]. e program structure is summarized in Figure 2.

e Python ECI Package and RRS Module.
is work modified and adapted the Fortran ECI library to produce a Python interface, compiling the modified ECI library using the F2PY program [15], creating a compiled Fortran library that can be imported by a Python program or module. It allows for the creation of ECI-compatible Python programs.
is package contains items 1 and 2 as outlined above and allows for the development of items 3-5 as a Python program.
e actual package consists of the following components: (1) e distributed ECI library with modifications to improve compatibility with Python, such as improved handling of command line arguments being passed through Python. (2) A Fortran module named PyInterface, which includes the subroutines to be exposed to Python using F2PY. ese subroutines call the ECI library's internal subroutines. is interface also includes an allocatable array, which is used to facilitate data transfers.
(3) An object-oriented Python module is named __init__.py. It introduces an object named Varia-bleData, which contains and manages ECI-linked variables. is object manages the Fortran allocatable array along with the data transfer requests.
A Python program using this package will follow the same program flow as an ECI-enabled Fortran program described in the previous section and in Figure 2. In summary, such a program will set up all of its data requests using a VariableData object. en, it executes a Python Science and Technology of Nuclear Installations version of TimeEvolve, following the same time step structure shown in Figure 1. Each requested component variable, corresponding to a TRACE component, is accessed using a Variable object. is object-oriented approach adds a layer of abstraction between the program and the underlying Fortran array as the program can simply call "get" and "set" methods to retrieve or change the data. e full implementation is summarized in Figure 3, showing how each component of the program is linked and how data is passed through ECI and through the Python ECI package. e ECI library also includes a number of global variables used for time step control, program flow control, and error reporting. Unlike component variables, these global variables are directly exposed by the interface created by F2PY and can be accessed directly by the program or through the ECI package. e ECI package is used to couple a reactor regulating system (RRS) model to the reactor model, along with a model for shutdown system 1 (SDS1) trips. e program is written in a modular fashion so that each RRS module (e.g., the liquid zone control) and each set of physical devices (e.g., the actual liquid zone compartments) are represented by individual Python objects. ese Python objects are linked to ECI variables along with each other. Figure 4 outlines the RRS module, showing the implemented components in blue and the ECI-linked input variables in orange.
During each time step, each control system function is modeled. Physical devices advance their state on each time step, while control modules update on a fixed interval (0.5 seconds for most RRS modules, and 2.0 seconds for liquid zone spatial control) to match the discrete time step nature of the digital RRS update frequency in CANDU reactors. High-level structure of ECI program. From the entry point, the program must call ECI library subroutines to initialize ECI coupling, including scheduling data transfers. e program must then call the TimeEvolve subroutine, which is responsible for synchronization and all program-specific computations.  SDS1 trips are also checked on every time step (because the SDS systems can actuate at any instant and are not subjected to the 0.5 s discrete RRS response limits). Much of the implementation is based on the description in [16]. Both backup and restart capabilities are implemented. If TRACE requests a time step back-up, the states of the RRS modules are reverted to their old-time values. Restart capability is implemented by writing the RRS state to a JSONformatted file at the end of the simulation and reading the state at the start of a restarted run. e following simplifications were made in the current version of the model: (i) No flux mapping routine (FLU) or fully instrumented channel (FINCH) mapping is included for spatial flux reconstruction. Instead, zonal powers are determined directly from the relevant PARCS nodes, equivalent to every channel being a FINCH with no measurement lag time. (ii) CANDU Setback is an RRS function that gradually reduces reactor powers and can be triggered from multiple signals, including from neutronic or process parameters. In this work, only neutronic setbacks are considered. Furthermore, setback on high local flux is not implemented as its true implementation relies on FLU [16]. erefore, setback is only implemented on high zonal flux and high flux tilt. (iii) Stepback in a CANDU is an RRS response to a large detected perturbation and causes a large reduction in rector power over a short interval of time.
Stepback triggers were implemented for neutronics parameters only, along with a manual stepback to 60% (to simulate a turbine trip Stepback). (iv) SDS1 trip triggers were implemented only for high neutron power, high log rate, low inlet feeder flow, and high heat transport system (HTS) pressure. e low flow trip uses the total flow rather than individual channel flows since the TRACE model aggregates the channels in groups of 60. (v) e platinum detectors are treated as perfect detectors that only respond to the thermal neutron flux with no delayed components.    Science and Technology of Nuclear Installations (vi) e ion chambers are not modeled, and the reactor power and log-neutron-rate signals from PARCS are used in its place. e PARCS reactor power is also used for the SDS1 high power trip.
Several additional Python programs are coupled with ECI, which are as follows: (1) A power calculator that reads the nodal powers and calculates channel and zone powers before passing them back to TRACE. It makes these powers visible to ECI, SNAP, and AptPlot. ECI is utilized here as the summation blocks built into TRACE currently use an inefficient implementation for data retrieval. (2) A data collection process that reads specified ECI variables every time step and outputs them to a tabdelimited data file for analysis. (3) A profiler that serves no purpose, except to measure the time spent waiting at each synchronization point. is data can be used to optimize the model's runtime by revealing which steps take up the most CPU time.

Coupling of PARCS and ECI to TRACE.
For modeling a CANDU reactor, the PARCS and TRACE models were coupled using a MAPTAB file created to map the PARCS neutronic nodes with the TRACE fluid and heat structure cells. Each fuel channel is mapped to the corresponding representative thermal-hydraulic channel. Channel powers are summed to get the total average channel power for each group of 60 channels, while fuel and coolant properties are shared among all channels for a given representative channel. Both models use the same axial division (one bundle length per node/cell) so the axial mapping is 1 : 1.
Coupling can be specified using either of the two sets of cards: (1) VOLRMAP cards: used for the automatic mapping of the PARCS mesh to the TRACE VESSEL or CHAN component. It is applicable to LWR analysis. (2)  As the TABLE cards must specify the mapping for each individual PARCS mesh cell, a script was written to automate the writing of these cards based on the mapping in the specification, accounting for the flow direction of each flow pass.
ECI programs are coupled by adding a "task list" file. is file contains a list of all the processes to run along with command-line arguments. When TRACE is run, it will read the task list and start all of the necessary processes.
While the use of ECI coupling reduces the need to modify the original source codes, certain changes were required in this work. As mentioned in the introduction, the codes with these modifications are referred to as TRACE_Mac1.0 and PARCS_Mac1.0, respectively.
(i) One signal variable parameter had two conflicting data validation checks when used with the %CRSIG card. As the parameter was otherwise unused for this type of signal variable, the %CRSIG-specific validation check was removed. (ii) Signal variable type 119 (PARCS cell power) was found to not function. Additional code was implemented to make the signal variable functional. (iii) A new signal variable type was added to retrieve the PARCS thermal flux. (iv) e POWER component was modified so that, in a coupled PARCS_Mac1.0-TRACE_Mac1.0 model, noncoupled POWER components (i.e., those not mapped to PARCS) will function as they would in a standalone TRACE model. (v) e SpecExTrans ECI module in TRACE was modified to permit access to signal variable data rather than requiring indirect access through control blocks to reduce data propagation delays in the coupling scheme. (vi) e size of the ECI transfer buffer was increased to accommodate the transfers of larger arrays, such as individual bundle powers. (vii) An error that could periodically occur when calculating D 2 O properties was modified to force a time step backup instead of terminating the simulation.

e CANDU-6 IAEA Code Coupling
Benchmark. e IAEA "Numerical Benchmarks for Multiphysics Simulation of Pressurized Heavy Water Reactor Transients" [17] is a set of standardized tests for the coupled simulation of postulated pressurized heavy water reactor (PHWR) transients. e purpose of the benchmarks is to provide test problems that may be implemented in different physic and thermalhydraulic codes. e following events are simulated on a stylized CANDU-6 model as part of the benchmark: (1) Steady state: to establish initial conditions prior to a transient, the coupled model is run until convergence. e converged model can then be used as a starting point for performing transient analysis. A null transient is used to test if the initial conditions are sufficiently converged. e fully coupled TRACE-PARCS model is compared with other codes using internal or external coupling to determine the agreement in the converged steady-state solution.
(2) Adjuster absorber rod withdrawal: in this loss of regulation (LOR) transient, adjuster rods 7 and 14 are withdrawn at a rate of 10 cm/s. e simulation is run for 25 seconds. ere is no credit for SDS. e key result is the set of channels that exceeds its respective critical channel power [17]. with a specified pump speed profile. e simulation is run for 25 seconds. ere is no credit for SDS. e key result is the set of channels that exceeds its respective critical channel power [17]. (4) Inlet header break: in this loss of coolant accident (LOCA), there is a break in the inlet header 2 of 0.0645 m 2 over the first 0.1 seconds of the transient, with an external pressure of 1 atm. SDS-1 is triggered at 120% of nominal core power, and the rod drop follows a specified profile. e simulation is run for 5 seconds.
e key result is the maximum bundle enthalpy [17]. e reactor regulating system (RRS) is not credited for any of these transients, and thus, it does not need to be simulated in this model. e adjuster rod withdrawal can be carried out using the MOVE_BANK card in PARCS. e benchmark is not described in full detail here, however, a summary of key features follows. e CANDU-6 design has 380 fuel channels and a two-loop heat transport system, with each loop serving half the core (left and right halves). Each loop has two flow passes in opposite directions. Adjacent channels have opposite flow directions, forming a checkerboard pattern as shown in Figure 5. Each flow pass has its own inlet and outlet headers. In the thermal-hydraulic model, only the main circuit, pressurizer, and some loop and pass interconnects are modelled. Each flow pass is represented by seven thermal-hydraulic channels, for a total of 28 thermal-hydraulic channels for the core, as shown in Figure 6. e boilers have representative flow resistance on the primary side, however, the secondary side is modeled as a boundary condition with a constant heat transfer coefficient and temperature. e pressurizer is also simplified, being represented by a large pipe. A boundary condition of saturated liquid at 10 MPa is connected to the pressurizer to set initial conditions in the steady state model.
In both neutronic and thermal-hydraulic models, each fuel channel is axially divided into 12 nodes of one bundle length each. A radial reflector is present around the core with a width of slightly more than 2 channel pitches. ere is no axial reflection.
For the core neutronic model, the nominal thermal power is 2000 MW. e model is highly simplified, and most structural materials are not modeled. e burnup distribution provided in the benchmark includes a significant flux tilt. e control devices modeled include 21 adjuster rods, 28 shutoff rods, and 14 liquid zone controllers. Adjuster rods start inserted at their nominal position, liquid zone controllers start 50% filled, and shutoff rods start outside the core. e adjuster rods and liquid zone controller positions are not changed during steady state or transient analysis, except if specified by the specific transient. e benchmark specifications included a full set of branch structures and a reduced set, with most participants electing to use the reduced set. For this work, the set of cross-sections generated with SCALE, available in Ref. [3], was primarily used. e transients were also modeled with the benchmark crosssections that utilize a reduced branch structure for comparison purposes. e initial model for this work was developed in Ref. [3] based on the specification for the code coupling benchmark. In Ref. [3], external coupling via scripts was used for TRACE and PARCS coupling. For the transients in Ref. [3], each code is run for 0.1 seconds of simulation time per step, until the end of the transient simulation is reached. It is a leapfrogging calculation with data exchange occurring at a coarse time step (0.1 seconds), while each individual code can divide the coarse time step into multiple fine time steps as needed.
In this work, the model was modified to fully utilize TRACE-PARCS coupling and ECI coupling, and while RRS was not utilized, the ECI scripts were still tested for functionality.
is coupled mode permits a tighter degree of coupling, with data exchange occurring at every TRACE time step. Since data exchange occurs entirely in memory, file management is greatly simplified, as files do not have to be generated for every step of data exchange.

2.6.
e 900 MW CANDU Model and Simulations. e 900 MW CANDU model is based on work performed in Ref. [18] related to station blackout transients and was converted from RELAP5 to TRACE. In addition to the 480 fuel channels, different reactor powers, and loop flows, several other notable changes were made compared to the CANDU 6 model discussed above. Compared to the CANDU-6 benchmark, this model includes pressurizer level control using feed and bleed, pressure control using pressurizer heaters, and heat transfer by the pressure tube and calandria tube, and it models the secondary side of the steam generators, including pressure and level control. However, only eight representative fuel channels are modeled, two per core pass. Unlike the CANDU-6 model, the 900 MW CANDU model is not based on a specific benchmark.
(1) e time step size for RRS may be decoupled from the time step size of TRACE. As PARCS converges the flux each time it runs the eigenvalue calculation, the flux responds instantaneously to control device movements, regardless of the time step size taken by TRACE. erefore, the RRS model steps forward by a fixed value of 0.5 seconds for each time step taken by TRACE-PARCS. (2) A modified version of the power error calculation (CEP) module is used, which gives a power error proportional to the reactivity. erefore, reactivity devices will respond to converge the core reactivity to zero. Under transients, the RRS control goes back to ensure convergence on power.
(3) e setback, stepback, and reactor trip modules are disabled to avoid spurious triggers while the steady state is being converged. e tests that were performed with this model are summarized as follows: (1) Two functionality tests designed to ensure that the coupled model with RRS behaves as expected (2)  Prior to testing the coupled system against data, some simple tests were performed on the functionality of the coupled system. ese include a zone control failure (fill valve fails closed for one zone) and an adjuster rod pull (rods 1 and 9). ere are no reference results for these integrated tests. Hence, they are evaluated based on whether the observed behavior is consistent with the expected control system's behavior and CANDU phenomena.
e first integral tests examine the figure-of-eight flow oscillation phenomenon, where data is available from a scaled CANDU integral test facility, described in the next section, along with the theory and simulation data on a CANDU-6 facility [19]. e figure-of-eight flow oscillation is an instability characterized by low-frequency density waves (T �14 s) propagating through the system [19]. Under  e second case begins with nominal conditions (at 100%FP) and decreases the system pressure setpoint to 9300 kPa along with blocking the interconnects. e results of these simulations are compared with prior theory and simulation [19]. It is not a benchmark comparison. Hence, an exact match is not expected, however, it is expected that the phenomenon is reproduced under the conditions where it is expected and has the expected properties. e final transient that was simulated was a coupled loss of flow (LOF) event resulting from a loss of class IV power event that occurred at an operating CANDU station. To mimic the loss of power during the event, HTS circulation pumps, feed pumps, and steam generator feedwater pumps lose power at zero time, coincident with a turbine trip. e auxiliary feedwater pump is started two seconds later based on available information from the station. e turbine trip triggers a stepback to 60%FP [18].
It was necessary to add the condenser steam discharge valves (CSDVs), atmospheric steam discharge valves (ASDVs), and main steam safety valves (MSSVs) to the LOF model. e CSDVs are available for the first 13.5 seconds of the transient until the condenser vacuum is lost, while the ASDVs and MSSVs are available for the entire transient. e pressure thresholds for these valves are listed in Table 1. SDS1 is available so that the reactor can trip on high neutron power, high neutron log rate, low inlet feeder flow, or high HTS pressure. e following model simplifications and assumptions are present in this model: (i) e feedwater pumps are not explicitly modeled.
Hence, the feedwater flow rate is modeled as a boundary condition instead. e trip of the main (3) (1) (2) (9) seconds (equal to the interval of the stepback algorithm during a stepback), and the deceleration of the absorbers at the end of the stepback is treated as being instantaneous. (iv) Control absorbers and shutdown rods are modeled as being dropped by gravity with damping, with the damping set to match the insertion rate in [18]. (v) ASDV and CSDV flow areas were set to achieve the nominal flow rates in [18]. MSSVs are set to the flow area provided in [18]. (vi) e pressurizer level setpoint is fixed at 6.5 m to match the station data and reference simulations.

e RD-14M
Model. e RD-14M facility is an experimental facility designed to be a full-length but scaled model of the CANDU primary heat transport system, including 10 fuel channel simulators, which are scaled channels with seven fuel elements each (compared to 37 fuel elements for a typical CANDU bundle). e facility is used to model various phenomena related to the CANDU heat transport system and provides experimental data on these phenomena. e ability of computer codes to model these phenomena can be evaluated by modeling the RD-14M facility and comparing the simulation's results to the experimental data [20].
In this work, previously presented in Ref. [20], an RD-14M TRACE model was adapted to perform flow instability tests, simulating the same phenomenon as for the 900 MW CANDU model. e primary and secondary sides of the model are shown in Figures 11 and 12, respectively. e pressure boundary condition on Header 8, along with the emergency coolant injection system, both shown in Figure 11, were used in a previous study, simulating a header break experiment [10], however, they are isolated from the rest of the system in this work. e surge tank functions as a pressure boundary condition using the TRACE pressurizer component. is component automatically adds or removes energy at a specified rate to maneuver its pressure to the setpoint. On the secondary side, the feedwater flow rate and temperature along with the steam line outlet pressure are specified as boundary conditions, as shown in Figure 12.
ese remain fixed throughout the simulations performed in this work. ese simulations were performed to reproduce similar conditions as those in the experimental tests [21]. In these experiments, pump speed reductions and system pressure reductions were used to induce oscillations. Two different interconnect designs based on different scaling methodologies were tested, along with tests using no header interconnect. One interconnect, labeled "geometric similarity," was scaled based on conserving the momentum equation, while the other, labeled "dynamic similarity," was scaled based on conserving the ratio of the interconnect flow resistance and the heat transport system flow resistance. ese labels are shown in Figure 11. It should be noted that at most one of these interconnects is connected at any time.
In these experiments, flow oscillations were produced when no header interconnect was used. e "geometric similarity" interconnect stabilizes the system against a pressure reduction but not a pump speed reduction, while the "dynamic similarity" interconnect stabilizes the system against both reductions. In the experiments where oscillations occur, the oscillation period is roughly 19 seconds [21], which is longer than the 14 seconds expected for a CANDU-6 heat transport loop [19]. As RD-14M is a full-length loop, comparable oscillation periods are expected.
Simulations were performed following the experimental procedures based on the available details [21], and the results were compared with the experimental results. If the results differed significantly, then the sensitivity of the results to the test procedure and test parameters was evaluated. RD-14M TRACE_Mac1.0 simulations do not include PARCS coupling, though the perturbations were performed using an ECI-coupled script.

CANDU-6 IAEA Code Coupling Benchmark.
e CANDU-6 benchmark model described in the methodology section was simulated using the coupled PARCS_Mac1.0-TRACE_Mac1.0 for the following cases:  Successful functionality of the coupled model has been demonstrated, and the general behavior of the models is consistent with that of the original version of the model that used an external coupling script. e model and cases were run for the full SCALE-generated branch structure and the reduced "benchmark specifications" branch structure from [3]. e quantitative results, shown in Table 2, agree well with the referenced externally coupled results [3], with the differences of at most a few percent. e difference in the results between the two different sets of cross-section data agrees with the prior study, with the reduced branch structure cases reaching higher maximum powers.
Possible sources of differences between the two studies include the increased data transfer frequency from using the built-in coupling (especially for the LOCA scenario), along with a difference between code versions, preventing a direct comparison.
e referenced study included a sensitivity study to the information exchange frequency. Increasing this frequency from every 0.1 seconds to every 0.01 seconds for the LOCA simulation increased the power peak by roughly 45 MW, along with making it occur 0.04 seconds earlier [3]. An information exchange frequency of 0.05 seconds results in a very similar peak power every 0.01 seconds [3]. Figure 13 shows the evolution of the core power during LOCA. Voiding in the broken flow loop results in positive reactivity in half of the core, resulting in a power excursion with a power tilt toward the voided half of the core. As the shutdown rods take a couple of seconds to drop into the core, peak channel power occurs in the lower left quadrant of the core, as the power initially decreases at the top of the core before decreasing in the rest of the core. Figure 14 shows the end state of the broken flow pass for the LOCA transient. e outlet header and fuel channels are almost fully voided. Cladding temperatures are elevated in most of this flow pass, with Figure 15 showing their evolution in time. e clad temperatures in CHAN24 exceed 1000 K. Figure 16 shows the end state of the adjuster rod withdrawal transient. e reactor power is elevated throughout the entire core but with localized peaking where the adjuster rods have been withdrawn.
e original externally coupled model and the internally coupled model were run to compare the performance of the two models. e greatest performance benefit is found when  Science and Technology of Nuclear Installations running transient analysis, with a run time of 160 seconds for the internally coupled case, compared to 505 seconds for the externally coupled model-approximately a factor of 3 difference. is benefit arises from avoiding the need for frequent restarts as is needed for the externally coupled case. ere is less of a performance benefit for steady-state analysis. PARCS can be configured to skip TRACE time steps, though there is currently an arbitrary limit of 20 TRACE time steps per PARCS execution. Each time PARCS is executed by TRACE, it runs until convergence, adding significant computation time compared to a standalone TRACE run, as well as when compared to the externally coupled model, where only a small number of restarts and PARCS executions are required.

900 MW CANDU Model Tests.
As mentioned in the methodology section, three sets of tests were performed, which are as follows: (1) Two functionality tests to evaluate that the coupled model is working as expected (2) Flow instability tests that induce flow oscillations in an off-normal state (3) A loss-of-flow transient initiated by a loss of Class IV power As mentioned in the methodology section, two functionality tests were performed to ensure robust RRS behavior. e first test was performed to simulate an abnormal operating occurrence (AOO), where the liquid zone control value fails, resulting in the draining of the liquid zone. e transient begins with the fill rate for Zone 5 being set to zero. It causes the liquid zone to rapidly empty, as shown in Figure 17. e result is increasing the reactor power, shown in Figure 18, which the RRS responds to by increasing the fill rate of the other zones, particularly adjacent zones (Zones 6 and 12), for which the spatial power error is the greatest. When the zone power in Zone 5 exceeds 110% of its nominal value, a setback is triggered to bring the zone power below 105% as expected. e result is an average reactor power of approximately 95%. At this point, the reactor continues to operate in a new steady state, though the average liquid zone controller level gradually decreases because of the build-up of xenon-135. After 50 minutes (not shown), the average zone water level is approximately 20%. e adjuster driven transient begins with adjuster rods 1 and 9 being withdrawn at the maximum rate possible based on the adjuster drive motor design. e adjuster and absorber rod movements are shown in Figure 19. It causes a small bulk power excursion by roughly 1% along with a flux tilt, as shown in Figure 20. e RRS setback logic is triggered  first, with liquid zone levels increasing rapidly as shown in Figure 21, however, high zone powers then trigger the stepback logic multiple times, resulting in partial absorber rod drops and the reduction of bulk power to 80%. Afterwards, there is further setback to 60% full power because of zone tilt.
In the RRS simulator produced by this work, the RRS attempts to maintain the current reactor power after a stepback. As the reactor power is still decreasing, this results in a brief negative power error that triggers the withdrawal of the first two adjuster banks, as seen in Figure 19. e banks stop at 25% of the maximum withdrawal as the power error briefly exceeds positive 3%, triggering RRS to stop withdrawing adjusters and start reinserting one bank. As bank C includes one of the failed rods that is detected as being withdrawn, it tries to reinsert bank C, however, nothing happens as one rod is stuck out and the others are already fully inserted.
After the initial stage of the transient, RRS responds to the xenon transient by withdrawing the control absorbers (up to a 75% average liquid zone level) and then reducing the average zone level. e liquid zone levels and reactor power remain highly asymmetric because of the asymmetric adjuster withdrawal.
Overall, the functionality tests show behavior consistent with what is expected based on the design of the RRS model in this work. e flow instability tests using TRACE_Mac1.0 and PARCS_Mac1.0 were carried out by running a steady-state  Science and Technology of Nuclear Installations calculation with the header interconnects unblocked, and the reactor power set to the initial power for the respective case, either 108%FP or 100%FP. en, the transient run is performed using the steady state calculation's conditions as the initial conditions for the transient. For the transient runs, the header interconnects are blocked at time t � 0. en, for the system pressure reduction case, the pressure setpoint is set to 9300 kPa at time t � 0. In both cases, flow oscillations were produced. Figure 22 shows the base case for the 108%FP scenario. At this power, an average void fraction of 14% corresponds to an average flow quality of 2%. Flow oscillations grow until 340 seconds to a large amplitude, with these oscillations also appearing in the reactor outlet header (ROH) pressures. When the oscillations grow large enough, the coolant voiding in the core rapidly changes and moves from one end of the core to the other as the voiding at any given time is the greatest toward the outlets of the fuel channels of the lowerpressure flow pass. e result is oscillations in both the endto-end flux tilt and total reactor power at rates which RRS is incapable of compensating for. At this point, reactor setback, stepback, or trip would be expected, however, for this simulation, these are disabled from acting so that the continuation of the transient can be observed. e average system pressure increases as the coolant swells because of both the loss of regulation and the loss of cooling effectiveness. e RRS response includes the movement of adjuster rods and control absorbers, which adds absorbing material to the top half of the core, resulting in a top/bottom flux tilt. e combination of these effects leads to an overall increase in system damping, suppressing the oscillations. At this point, the outlet header void fractions are higher and the system pressure is lower than during the initial steady state. e average flow quality at this point is 4%. e oscillations begin to grow again as the outlet header void fraction decreases toward the initial steady-state value, as shown in Figure 23. e average flow quality is approximately 3% at the time which oscillations begin to grow again. It repeats every several minutes for as long as the simulation continues to run. e period of oscillations is dependent on the amplitude, ranging from 15 seconds at low amplitude to 12 seconds at high amplitude. ese are close to the 14 seconds predicted in literature [19]. e period is reduced at high amplitude as one voided region is fully compressed on each half-cycle, resulting in the oscillations "bouncing" off the incompressible liquid. Figure 24 shows the same case, except that the reactor power setpoint is reduced to 100%FP at 200 seconds. It reduces the amplitude of the oscillations, and they eventually stop after another 400 seconds. is figure also shows how the total reactor power begins to oscillate when the amplitude of the flow oscillations is large. Figure 25 includes a second power decrease to 90%FP at 300 seconds, which stops the flow oscillations more quickly than maintaining 100%FP. Figure 26 shows the base case until 200 seconds, at which point the header interconnects are unblocked. e result is that the oscillations are reduced to a low amplitude, as expected.  Figure 27 shows the base case for the reduced pressure scenario. e oscillations begin once the system reaches the 9300 kPa setpoint and behave similarly to the increased power case in Figure 22. Once again, the oscillations are eventually dampened by the increasing coolant voiding and the control rod movements. Figure 28 shows that reducing the reactor power to 90%FP stops the oscillations, even at the reduced system pressure. Figure 29 shows that connecting the balance headers also suppresses the oscillations. e final simulation using this model with TRACE_-Mac1.0 and PARCS_Mac1.0 is the loss of flow event initiated by a loss of Class IV power, as detailed in the methodology section. e most significant effects on this are the shutdown of the heat transport circulating pumps, along with a turbine trip. e turbine trip has the dual effect of quickly stopping steam flow out of the steam generators along with triggering a reactor stepback to 60%. e heat transport system feed pumps and the main boiler feed pumps also shut down. Table 3 summarizes the subsequent events that occur as a consequence of the transient. e low flow trip timing is very similar to the referenced study [18]. One significant difference is that the pressure in this study never reaches the threshold to actuate the HTS liquid relief valves. Figure 30 shows the effect of the transient on the core neutronics along with the circulating flow rate. Initially, the reactivity and fission power begin to increase as the void fraction near the channel outlets increases. However, the stepback quickly adds negative reactivity to being the reactor power down to 60%FP. e behavior of the stepback depends on the time extrapolation of the reactor power and the kinematics of the control absorbers. With the assumptions made in this model, the stepback is triggered briefly for a second time shortly after the first stepback ends as the power remains over 60%FP for a brief period. e impact of this secondary RRS action was minimal, however. A low flow trip is triggered at 3.4 seconds into the transient, with the shutdown rods being dropped into the core 0.3 seconds later. It reduces the reactor power to decay heat levels within 2 seconds. Figure 31 shows the effect of the transient on the heat transport pressure. e flow reduction results in coolant swelling, which rapidly increases the system pressure until the reactor trip reduces the rate of heat generated in the core, at which point the coolant shrinks and the system pressure decreases. e peak pressure in this study is lower than that in the previous studies, however, it is still comparable to the available station data. Figure 32 shows the effect on the steam generator pressures. e pressure rapidly increases as the turbine stop valves close and prevent the removal of energy from the steam generators. is initial behavior of the steam generator pressure (first 15 seconds) is similar to the results from Ref. [18], though at a somewhat higher pressure. e pressure peaks at approximately 5300 kPa as the ASDVs and CSDVs open. en, it further increases when the condenser becomes unavailable, for which the behavior in this study is closer to the station data than to the reference simulations, for both the value and timing of the peak pressure. It occurs at roughly 30 seconds into the transient, at 5440 kPa, at which point the energy input from the HTS matches the energy output through the ASDVs. e pressure then decreases at a similar rate to prior studies. However, while the pressure levels off at 5200 kPa in the station data, it continues to decrease toward the ASDV setpoint of 5085 kPa in this study. Figure 33 shows the effect of the transient on the reactor inlet header (RIH) temperature. e steady-state temperature is a function of the heat transfer efficiency of the steam generator. Typically, the steam generator model is tuned to achieve the desired RIH temperature, whether matching station data or a design value. In this study, no specific tuning was performed. For the real event, the measured temperature immediately prior to the transient was 263.8°C, which is very close to the temperature in this study and almost 2°C lower than the temperature in the previous studies. e temperature initially increases and peaks 10 seconds into the transient, consistent with the previous studies.
is increase is caused by the interruption of feedwater flow, which greatly reduces the heat transfer efficiency in the preheater section of the steam generators. Figure 34 shows the effect of the transient on the steam generator level. Overall, the trend is dominated by shrinkage in the steam generator as the rate of heat input from the HTS decreases, and the water level follows a similar trend to the station data and RELAP results. e steam generator water level also decreases because of the mass imbalance caused by the loss of the boiler feed pumps. e auxiliary feed pumps activate, however, their flow rate is lower than the rate of steam flow through the ASDVs.

RD-14M Model Flow Instability Tests.
e RD-14M flow transient simulations [20] using TRACE_Mac1.0 were performed to model the figure-of-eight flow oscillations similar to those simulated in the 900 MW CANDU model. In this case, the reference data are the results from the experiments performed at the actual facility [21]. As in the 900 MW CANDU model, flow oscillations could occur once the system state was changed in a way that induced an increase in void fraction in the outlet headers. In this case, this was done either by reducing the pump speed or by reducing the system pressure. e initial pump speed and system pressure are 75% nominal speed and 10 MPa, respectively.
With no header interconnect, when the pump speed is reduced, flow oscillations can be induced as in the experiments. One significant difference between the simulations and experiments was that a gradual pump speed reduction would not induce growing oscillations, while an abrupt change would. It suggests that in the simulation, the initial perturbation is important, where, for a given pump speed, small oscillations die out while larger oscillations grow. It is not observed in the experimental results that induce large oscillations when the pump speed is reduced in several smaller increments [21]. However, there are no experimental results available for a true gradual reduction of the pump speed. us, while a gradual pump speed reduction does not reproduce the experimental results, an abrupt reduction in the pump speed results in growing oscillations up to a large amplitude, similar to the experimental data, as shown in Figure 35. e oscillation period is roughly 15 seconds, comparable to the results for the 900 MW CANDU model, however, it is significantly different from the 19-second value from the RD-14M experiments [21]. Figure 36 shows a test, where, with no header interconnect, the system pressure was reduced in increments, temporarily isolating the surge tank after each increment.
is test does not directly match any of the experimental tests but evaluates the same results as one of the tests in a more quantitative manner. e test shows that the system is unstable between 9.3 and 9.6 MPa while the surge tank is isolated. In the experiments, this range was instead between 9.5 and 9.8 MPa [21]. While differing quantitatively, there is a qualitative agreement, showing that there is a pressure range, and hence a corresponding outlet header void fraction range makes the system unstable.
e system is stable outside of this pressure and void fraction range in either direction.
is observation also supports the theory discussed in [19].
When the "geometric similarity" header interconnect is used, the simulations and experiments are comparable, showing that the interconnect is effective when the system pressure is reduced but not when the pump speed is reduced. e results of the system pressure reduction are shown in Figure 37. A temporary oscillation occurs when the surge  tank is isolated but dies out after several oscillation periods. However, in the simulations, a pump speed reduction to 68%, as performed in the experiments, did not result in growing oscillations. Instead, a simulation was performed, which reduced the pump speed first to 60%, then to 55%, then increasing it to 64% [20], with Figure 38 showing the oscillations subsequent to this final increase. At this pump speed, the oscillation amplitude gradually decreases until the surge tank is isolated at t � 1320 s, which causes the oscillations to grow again until the surge tank is reconnected at t � 1440 s. e pump speed is increased to 66% at t � 1600 s, which causes the oscillations to decay more rapidly.
When the "dynamic similarity" header interconnect is used, the simulations and experiments show that the system is stable to the figure-of-eight oscillations for both a system pressure reduction and a pump speed reduction. Figure 39 shows that this interconnect is more strongly stabilizing for system pressure reduction than the "geometric similarity" interconnect ( Figure 37). For the pump speed reduction, in addition to the experimental procedure of reducing the pump speed gradually to 66% nominal speed, a more extreme simulation was performed, with a step reduction to 55% nominal speed. e initial oscillations resulting from this perturbation decay over several oscillation periods are as shown in Figure 40.

22
Science and Technology of Nuclear Installations        Overall, these results, combined with the results on the 900 MW CANDU model, suggest that TRACE_Mac1.0 code can simulate the CANDU figure-of-eight flow oscillation phenomenon (while the work was carried out using TRACE_Mac1.0, the ability for TRACE as distributed to model this phenomenon should be identical.). However, quantitative differences between the RD-14M simulations and experiments suggest that the current RD-14M TRACE model contains inaccuracies that affect the simulation of flow oscillations. e investigation of these inaccuracies is beyond the scope of this work.

Conclusions
is work demonstrates that the built-in coupling between PARCS_Mac1.0 and TRACE_Mac1.0, when combined with coupling additional models using ECI, can be used to model a CANDU unit with minimal changes and additions to PARCS and TRACE themselves.
is work was able to quantitatively reproduce the results obtained for the same CANDU-6 model running with external coupling while improving the computational efficiency of the model on the whole, particularly for transient analysis. In addition, successful coupling with a reactor regulating system model using ECI was demonstrated as it was able to qualitatively reproduce the behavior of CANDU RRS and show similar results to prior analyses with RELAP5 and other codes.
Certain source code modifications were required to achieve the desired coupling, most notably, additional signal variable functionality for coupled PARCS_Mac1.0-TRACE_Mac1.0 simulations, direct ECI coupling to signal variables, enabling noncoupled POWER components to function in a coupled simulation, and enabling the PARCS-TRACE coupling of control device positions in TRACE_Mac1.0.
As additional programs were developed modularly, the programs can be adapted to different applications. In particular, the ECI Python package is completely independent of the RRS module and its components, and therefore, they may be developed separately, and several other auxiliary programs were developed utilizing the ECI Python package. Both programs may be adapted to future applications, including the simulations of reactor operation, core follow analysis, other safety analysis cases, and uncertainty analysis.

Data Availability
e TRACE and PARCS model data used to support the findings of this study have not been made available due to the confidentiality of the models.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper. Figure 40: Flow oscillations with reduced pump speed with dynamic similarity interconnect. Experiment inset [21].