Abstract
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 the interactions between thermal-hydraulic and physic phenomena with the CANDU reactor regulating system (RRS). To date, no or limited development has taken place in TRACE-PARCS in this area. In this work, the system thermal-hydraulic code 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 900 MW 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 900 MW 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 2x to 3x less execution time. The 900 MW CANDU simulations successfully demonstrate the RRS functionality for the loss of flow events, and the coupled solutions demonstrate adequate performance for figure-of-eight flow instability modeling.
1. Introduction
In nuclear power plants, including CANDU stations, many phenomena arise because of coupled interactions between the reactor physic (nuclear chain reaction) and system thermal-hydraulic phenomena. These range from simple reactivity feedbacks to 3D spatial power changes to coupled flow instabilities. For this reason, many core physic and system thermal-hydraulic codes include coupling capabilities. This 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. The 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. This 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. The 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]. The 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. This study focuses on the implementation of RRS response via the ECI interface.
This 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). This 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. This 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. The changes are relatively minor, primarily to facilitate code coupling for this work’s applications, and they do not add or modify any physical models.
2. Materials and Methods
2.1. PARCS-TRACE Coupling Background
The PARCS core physic code is included within the distribution of TRACE. The two codes are compiled into a single executable.
The 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. The MAPTAB file is a user-created 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. The updated core power data is passed back to TRACE for the next time step. The accuracy of this coupling is equivalent to externally controlled coupling for a given data transfer frequency. This 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 thermal-hydraulic 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] in the MAPTAB file, followed by an arbitrary number of signal variables to rod bank pairs. On each time step, the values of signal variables are mapped to the position values of the corresponding rod banks. The PARCS upper and lower limit values are not utilized. Thus, the appropriate conversion must be incorporated on the TRACE/ECI side if rod banks have different amounts of travel, which is the case for CANDU reactivity devices.
Therefore, there are three ways to manipulate or control banks during a PARCS transient, which are as follows:(1)Specified bank movement table that updates bank positions as a function of time. This function runswithin PARCS without interaction from TRACE (other than TRACE determining the time step size).(2)Specified SCRAM bank movement table that drops in all banks or a subset of banks. The SCRAM may be triggered within PARCS by relative power and/or be triggered by TRACE when the PARCS input specifies a TRACE trip signal.(3)Fully coupled to TRACE signal variables using the CRSIG card.
This research focuses on method 3, which is the only method that can dynamically react to the rest of the system in a generalized fashion.
2.2. Exterior Communication Interface (ECI) Overview
The 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.
The two main components of the ECI are as follows [12]:(1)The 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)The 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.
The 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. Thus, 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.
The 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]. The 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)The ECI library distributed with TRACE(2)A completed version of the SpecExTrans module if its variables are to be visible to other processes(3)The 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. They may be found in the ECI manual [12]. The program structure is summarized in Figure 2.
2.3. The Python ECI Package and RRS Module
This 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. This package contains items 1 and 2 as outlined above and allows for the development of items 3–5 as a Python program.
The actual package consists of the following components:(1)The 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. These subroutines call the ECI library’s internal subroutines. This 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 VariableData, which contains and manages ECI-linked variables. This 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. Then, it executes a Python 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. This 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. The 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.
The 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.
The 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. The 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. These 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. 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 JSON-formatted file at the end of the simulation and reading the state at the start of a restarted run.
The 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]. Therefore, 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. The low flow trip uses the total flow rather than individual channel flows since the TRACE model aggregates the channels in groups of 60.(v)The platinum detectors are treated as perfect detectors that only respond to the thermal neutron flux with no delayed components.(vi)The ion chambers are not modeled, and the reactor power and log-neutron-rate signals from PARCS are used in its place. The 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 tab-delimited data file for analysis.(3)A profiler that serves no purpose, except to measure the time spent waiting at each synchronization point. This data can be used to optimize the model’s runtime by revealing which steps take up the most CPU time.
2.4. 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)TABLE cards: used for the manual mapping of PARCS mesh cells to TRACE fluid and HTSTR components. For CANDU analysis, it is necessary to use the TABLE cards.
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. This 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)The 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)The 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)The 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 D2O properties was modified to force a time step backup instead of terminating the simulation.
2.5. The CANDU-6 IAEA Code Coupling Benchmark
The 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. The purpose of the benchmarks is to provide test problems that may be implemented in different physic and thermal-hydraulic codes. The 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. The 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. The 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. The simulation is run for 25 seconds. There is no credit for SDS. The key result is the set of channels that exceeds its respective critical channel power [17].(3)Coolant pump rundown: in this loss of flow (LOF) transient, heat transport pump #2 begins a rundown with a specified pump speed profile. The simulation is run for 25 seconds. There is no credit for SDS. The 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 m2 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. The simulation is run for 5 seconds. The key result is the maximum bundle enthalpy [17].
The reactor regulating system (RRS) is not credited for any of these transients, and thus, it does not need to be simulated in this model. The adjuster rod withdrawal can be carried out using the MOVE_BANK card in PARCS.
The benchmark is not described in full detail here, however, a summary of key features follows. The 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. The 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. The 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. There is no axial reflection.
For the core neutronic model, the nominal thermal power is 2000 MW. The model is highly simplified, and most structural materials are not modeled. The burnup distribution provided in the benchmark includes a significant flux tilt. The 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. The adjuster rods and liquid zone controller positions are not changed during steady state or transient analysis, except if specified by the specific transient. The 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. The transients were also modeled with the benchmark cross-sections that utilize a reduced branch structure for comparison purposes.
The 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 leap-frogging 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. This 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. The 900 MW CANDU Model and Simulations
The 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. Thus, the comparison of the model results to prior work is more qualitative in nature.
Figure 7 shows the nodalization of the neutronic calculation and its mapping to the thermal-hydraulic model in Figure 8, showing how the eight representative fuel channels map to the 480 actual fuel channels. Axially, both models have fuel channels that are divided into 12 nodes of one bundle length each, with no axial reflector regions in the neutronic model. On the primary side, shown in Figure 8, the fluid boundary conditions are the feed and bleed flows along with the pressurizer steam bleed. The steam generator models are shown in Figure 9. The primary components of each steam generator, as shown in Figure 10, are the upper and lower boiler regions, the preheater, the downcomer, and the steam drum. The details of steam separation are not implemented in this model. Instead, the steam drum is modeled as a single large node, modeling the steam generator liquid inventory. Only dry steam is permitted to flow into the steam dome node, while only liquid is permitted to flow into the downcomer.
The reactor physic model includes liquid zone controllers, adjuster rods, mechanical control absorbers, and shutoff rods. All devices are in their normal positions during a steady state calculation, except for the liquid zone controllers, which are adjusted to control reactivity and power distribution according to the liquid zone control algorithm, which is part of the RRS model. During transients, the reactivity devices are controlled by the RRS and respond based on the RRS program rules.
For steady-state calculations, where PARCS is running an eigenvalue calculation, the goal is to converge the zone levels to their equilibrium positions as quickly as possible, not to simulate the actual dynamics of the zone control system. Since the power level is fixed, keff replaces reactor power as the variable for bulk control to converge. To achieve the goals for steady-state analysis, there are several differences in the operation of the RRS model in a steady-state analysis.(1)The 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. Therefore, 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. Therefore, reactivity devices will respond to converge the core reactivity to zero. Under transients, the RRS control goes back to ensure convergence on power.(3)The setback, stepback, and reactor trip modules are disabled to avoid spurious triggers while the steady state is being converged.
The 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)Two-phase flow instability tests that induce figure-of-eight flow oscillation in an off-normal state(3)A loss of flow station transient, where the initiating event is a loss of Class IV power
Prior to testing the coupled system against data, some simple tests were performed on the functionality of the coupled system. These include a zone control failure (fill valve fails closed for one zone) and an adjuster rod pull (rods 1 and 9). There 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.
The 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]. The figure-of-eight flow oscillation is an instability characterized by low-frequency density waves (T = 14 s) propagating through the system [19]. Under certain conditions, normal perturbations can result in diverging oscillation, eventually leading to reactor setback, stepback, or trip. The instability is very sensitive to reactor configuration, with the most recent CANDU plants designed to preclude this instability using a loop interconnect. To generate the instability, the 900 MW CANDU model was modified as follows:(i)Removal of the interconnects within each loop. The loop interconnect is a balance line that promotes stability in the two-loops and was adopted in later designs like CANDU-900. The pressurizer connection between the two loops is maintained for all cases.(ii)Setback, stepback, and reactor trip are disabled.(iii)The increasing of the outlet header quality above the nominal levels to roughly 3–4% of the flow quality. It can be done by increasing reactor power (to 108% full power in this study) or reducing system pressure. Both cases are demonstrated.
In the TRACE model, the interconnects are replaced with valve components. These are open during the initial steady state convergence but closed at the start of the transient. The interconnects remain closed, except for the subcases, where the interconnects are reopened to test the effect of the interconnects at dampening flow oscillations.
The first case begins with a steady-state calculation at 108% full power (FP) and then disconnects the interconnect line. Once large oscillations develop, two subcases are tested to see if the oscillations can be stopped: unblocking the interconnects and lowering the reactor power.
The second case begins with nominal conditions (at 100%FP) and decreases the system pressure setpoint to 9300 kPa along with blocking the interconnects.
The 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.
The 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. The auxiliary feedwater pump is started two seconds later based on available information from the station. The 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. The 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. The 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.
The following model simplifications and assumptions are present in this model:(i)The feedwater pumps are not explicitly modeled. Hence, the feedwater flow rate is modeled as a boundary condition instead. The trip of the main feedwater pump was simulated separately, and the results were applied to the feedwater flow rate in the main model. The auxiliary feedwater pump is modeled as a flow rate of 10 kg/s per steam generator, starting from 12 seconds after they are started up (i.e., 14 seconds after the start of the transient).(ii)Only the main circuit of the HTS is modeled, with the feed and bleed modeled as flow rate boundary conditions. Feed flow is disabled at the start of the transient. Bleed flow uses values taken from station data.(iii)The stepback forecasts the reactor power by 0.25 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)The pressurizer level setpoint is fixed at 6.5 m to match the station data and reference simulations.
2.7. The RD-14M Model
The 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). The facility is used to model various phenomena related to the CANDU heat transport system and provides experimental data on these phenomena. The 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. The primary and secondary sides of the model are shown in Figures 11 and 12, respectively. The 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. The surge tank functions as a pressure boundary condition using the TRACE pressurizer component. This 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. These remain fixed throughout the simulations performed in this work.
These 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. These 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. The “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.
3. Results and Discussion
3.1. CANDU-6 IAEA Code Coupling Benchmark
The CANDU-6 benchmark model described in the methodology section was simulated using the coupled PARCS_Mac1.0-TRACE_Mac1.0 for the following cases:(1)Steady state(2)Inlet header break—loss of coolant accident (LOCA)(3)Single coolant pump rundown—loss of flow (LOF)(4)Adjuster rod withdrawal—loss of regulation (LOR)
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. The model and cases were run for the full SCALE-generated branch structure and the reduced “benchmark specifications” branch structure from [3].
The quantitative results, shown in Table 2, agree well with the referenced externally coupled results [3], with the differences of at most a few percent. The 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. The 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.
(a)
(b)
(c)
(d)
Figure 14 shows the end state of the broken flow pass for the LOCA transient. The 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. The clad temperatures in CHAN24 exceed 1000 K.
Figure 16 shows the end state of the adjuster rod withdrawal transient. The reactor power is elevated throughout the entire core but with localized peaking where the adjuster rods have been withdrawn.
(a)
(b)
(c)
The original externally coupled model and the internally coupled model were run to compare the performance of the two models. The greatest performance benefit is found when 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. This benefit arises from avoiding the need for frequent restarts as is needed for the externally coupled case.
There 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.
3.2. 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.
The 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. The 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. The 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. The 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%.
The adjuster driven transient begins with adjuster rods 1 and 9 being withdrawn at the maximum rate possible based on the adjuster drive motor design. The 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. The 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. The 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. The 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.
The flow instability tests using TRACE_Mac1.0 and PARCS_Mac1.0 were carried out by running a steady-state 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. Then, 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. Then, 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 lower-pressure flow pass. The result is oscillations in both the end-to-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. The average system pressure increases as the coolant swells because of both the loss of regulation and the loss of cooling effectiveness. The 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. The 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. The average flow quality at this point is 4%. The oscillations begin to grow again as the outlet header void fraction decreases toward the initial steady-state value, as shown in Figure 23. The 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.
The period of oscillations is dependent on the amplitude, ranging from 15 seconds at low amplitude to 12 seconds at high amplitude. These are close to the 14 seconds predicted in literature [19]. The 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. This 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. The result is that the oscillations are reduced to a low amplitude, as expected.
Figure 27 shows the base case for the reduced pressure scenario. The 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.
The 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. The most significant effects on this are the shutdown of the heat transport circulating pumps, along with a turbine trip. The turbine trip has the dual effect of quickly stopping steam flow out of the steam generators along with triggering a reactor stepback to 60%. The 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. The 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. The 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. The 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. The 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. The 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. The pressure rapidly increases as the turbine stop valves close and prevent the removal of energy from the steam generators. This initial behavior of the steam generator pressure (first 15 seconds) is similar to the results from Ref. [18], though at a somewhat higher pressure. The pressure peaks at approximately 5300 kPa as the ASDVs and CSDVs open. Then, 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. The 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. The 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. The temperature initially increases and peaks 10 seconds into the transient, consistent with the previous studies. This 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. The steam generator water level also decreases because of the mass imbalance caused by the loss of the boiler feed pumps. The auxiliary feed pumps activate, however, their flow rate is lower than the rate of steam flow through the ASDVs.
3.3. RD-14M Model Flow Instability Tests
The 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. The 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. Thus, 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. The 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. This 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. The 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. The system is stable outside of this pressure and void fraction range in either direction. This 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. The 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. The 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. The initial oscillations resulting from this perturbation decay over several oscillation periods are as shown in Figure 40.
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. The investigation of these inaccuracies is beyond the scope of this work.
4. Conclusions
This 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. This 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
The 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
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This work was supported by the University Network of Excellence in Nuclear Engineering (UNENE) and the Natural Sciences and Engineering Research Council of Canada (NSERC). The authors would like to thank Feng Zhou, Kai Groves, and Michael Tucker for their prior work in developing the TRACE and PARCS CANDU models which were used in this research. They would also like to thank Anatol Mysen and David Hummel for their prior work in developing the TRACE RD-14M model used in this research.