Analysis of Shutdown System Effectiveness in the Canadian Super Critical Water Reactor Using Coupled Thermal Hydraulics and Three-Dimensional Neutron Kinetics

The pressure tube supercritical water-cooled reactor (PT-SCWR) is a generation IV reactor concept developed by Canada Nuclear Laboratories (formerly Atomic Energy of Canada Limited (AECL)) in collaboration with Natural Resources Canada and the Natural Sciences and Engineering Research Council of Canada. As an evolution of the Canada Deuterium Uranium (CANDU) reactor, the PT-SCWR features both pressure tubes and a low-temperature heavy water moderator. Unlike CANDU, however, the PT-SCWR concept makes use of a batch refueling scheme (as opposed to online refueling), vertical fuel channels, and light water coolant [1].

A key component of the PT-SCWR concept is the high efficiency re-entrant channel (HERC) [2]. As illustrated in Fig. 1, coolant enters a central flow tube at the top of each channel and travels downward before reversing and flowing upward around the fuel elements. Isolation between the supercritical light water coolant and the low temperature heavy water moderator is provided by a ceramic insulator and zirconium-alloy pressure tube. The fuel itself is a PuO₂–ThO₂ mixture (15 wt % and 12 wt % PuO₂ in the inner and outer rings, respectively) contained within zirconium modified stainless steel cladding [3].

While the reactivity effect of entirely voiding the light water coolant within the channels is negative (as required by the design specifications), the same cannot be said for void exclusively around the fuel. In fact, lattice-level neutron transport calculations have shown that the reactivity worth of void around the fuel is positive (as high as 10 mk), and only the very large negative reactivity of void within the center flow tube (>40 mk) makes the net effect of voiding both regions negative as required. Previous work demonstrated that this positive reactivity worth leads to brief power excursions in transients with flow reductions (where the coolant density around the fuel decreases before changes propagate to the center flow tube), including loss-of-coolant accidents (LOCAs) and loss-of-flow accidents [4]. A fast-acting shutdown system is potentially necessary to counteract this positive reactivity insertion and mitigate consequences of these design basis accidents.

The objective of this work was thus to model the response of the PT-SCWR core and assess a fast-acting shutdown system's capability to arrest power excursions resulting from transients in the primary heat transport system (PHTS). This would require core-level modeling of coupled thermal-hydraulics and spatial kinetics. The characteristics of a fast-acting shutdown system widely deployed in CANDU reactors (in particular, shut down system 1—rod insertion) would serve as the basis for a similar system in the PT-SCWR for the purpose of this study.

A core-level coupled model of PT-SCWR thermal-hydraulics and spatial kinetics was previously created to study the transient behavior of the system in the absence of any reactivity control and shutdown systems [4]. This model, with the addition of a fast-acting shutdown system, would also serve as the basis for this work.

The PT-SCWR core neutronics were modeled with the DRAGON and DONJON computer codes. DRAGON is an open-source code developed at École Polytechnique de Montréal that solves the multigroup neutron transport equation with burnup in two and three dimensions [5]. It was used to generate homogenized and condensed cross section and a database of thermal-hydraulic feedback for input in to a DONJON model. DONJON is also an open-source code developed at École Polytechnique de Montréal, and it is used to solve neutron diffusion and spatial kinetics equations in three dimensions [6]. The DONJON model was used to calculate the steady-state and transient power distributions in the PT-SCWR core.

The PT-SCWR reference lattice cell (as shown in Fig. 1) was modeled with DRAGON 3.06 J using the International Atomic Energy Agency's 172-group evaluated nuclear data library [7]. To account for the axial variation in material temperatures and densities, as well as the heterogeneity near the boundary of the core and the heavy water reflector, 60 separate sets of homogenized and condensed cross section as functions of burnup were generated at the reference operating conditions (corresponding to 20 axial locations for each of the infinite lattice/interior, side, and corner lattice cells). The reference temperatures and densities at each axial location and burnup step were then perturbed to generate a database of first- and second-order coefficients that describe how each homogenized and condensed value varies with the thermal-hydraulic conditions. The reference cross section and feedback database were necessary inputs to the core-level DONJON model.

Exploiting the quarter-symmetry of the reference batch refueling sequence, the DONJON 3.02 G model of the PT-SCWR consisted of 84 channels with 20 axial nodes, each node corresponding to a single lattice cell calculation as described previously. A 100 cm radial reflector and 75 cm axial reflector (both D₂O) were also included. This geometry is shown in Fig. 2. The model consisted of 3460 uniformly sized cubic nodes (1680 for fuel, 1820 for reflector) and used eight energy groups for the diffusion calculation. The same kinetics parameters (Table 1), taken from a previous standardized computer analysis for licensing evaluation/transport rigor implemented with time-dependent operation for neutronic depletion (SCALE/TRITON) calculation, were used for each core snapshot [8].

Note that this reference model contains no control devices or other neutron absorbers for reactivity hold-down during the batch cycle. No such devices were finalized within the PT-SCWR reference conceptual design at the time of this work. Nevertheless, this work assumed that some reactivity hold-down and channel power-shaping system would mitigate the local power variation that naturally occurs during a batch cycle as fuel depletes (an example of a proposed reactivity hold down and power shaping system is presented in Ref. [9]). Such a system would ensure that the thermal-hydraulic conditions along the length of the channel remain at the reference values that were the basis for the cross section database calculation (i.e., in each lattice cell depletion calculation, the thermal-hydraulic conditions were assumed to be constant over the length of the batch cycle). With this assumption, an “equilibrium” core could be determined where the transient value of burnup and power in each node was the same from cycle to cycle. A “snapshot” of the equilibrium core (corresponding to the beginning, middle, or end of the refueling cycle, respectively 0, 224, and 448 full power days) serves as the steady-state initial condition for subsequent kinetics calculations.

CATHENA (Canadian Algorithm for Thermalhydraulic Network Analysis) is a thermal-hydraulics code developed by AECL primarily for LOCA analysis of CANDU reactors [10]. It uses a one-dimensional, two-fluid representation to calculate transient flows in piping networks. The included generalized heat transfer package (GENHTP) calculates convective heat transfer to the fluid, conduction within walls, and radiative heat transfer between surfaces. Version CATHENA MOD-3.5d/Rev 3, used in this work to model the thermal-hydraulics of the PT-SCWR heat transport system, contains an expanded set of fluid properties for supercritical conditions but no specific heat transfer or friction correlations (standard correlations for subcritical water were used in lieu of dedicated models) [11]. CATHENA MOD-3.5d/Rev 3 has nevertheless been used for much of the safety analysis of the preconceptual PT-SCWR, so its continued use in this study is considered warranted [4,11].

The CATHENA idealization of the PT-SCWR is shown in Fig. 3. Pressure and temperature boundary conditions are placed at the core inlet and outlet, omitting the coolant pumps, turbines, and other elements of the primary circuit for simplicity. The 84 channels in each quarter core were modeled as four identical parallel pipes, thus giving the correct total flow for the entire core. Each channel was modeled with 20 axial nodes to ensure one-to-one matching with the DONJON model, from which the axial and radial power distributions could be imported as heat deposited in the fuel. Convective heat transfer was modeled from fuel to coolant, as well as from coolant to the center flow tube and pressure tube (with conduction through the liner and ceramic insulator). Radiation heat transfer was also modeled between the cladding surfaces, center flow tube, and liner tube in three select channels. These included the highest power channel (E6), lowest power channel (C4), and an average power channel (G8).

Flow-limiting orifices are located at the inlet of each channel (i.e., the connection between the inlet plenum and center flow tubes) that match each channel's flow to its power, thus ensuring a uniform channel outlet temperature, $Tout$⁠, of 625 °C (note that 625 °C and 25 MPa at the outlet are design goals for the PT-SCWR [1]). This flow-power matching is a key design criterion for the PT-SCWR concept, wherein each channel's power is held nearly constant by a reactivity control and power shaping system that compensates for fuel depletion during the length of the fuel cycle. Such a system would likely consist of some combination of burnable neutron absorbers and control rod movements. This work achieved flow-power matching by modeling each orifice as a valve attached to a proportional-integral controller, which automatically finds the orifice size that provides $Tout$ = 625 °C at each channel outlet. The orifice sizing calculation was repeated for each core snapshot and then frozen for subsequent kinetics calculations that use that snapshot as the initial condition. Fixed orifices would be present in the real PT-SCWR, but without a final power shaping system being part of the conceptual design, variable orifices were necessary to model the PT-SCWR core at this stage (note that similar modeling approaches have been applied previously [4,8]).

The procedure for coupling DONJON and CATHENA in transient simulations was described in Ref. [4]. Within a single time-step, DONJON first calculates the new power in each node based on the current macroscopic cross section and kinetics parameters. A script reads this power distribution and creates new CATHENA input with the corresponding heat generation in the fuel. CATHENA is then executed to generate new material densities and temperatures in each node over the same time-step. The script reads these distributions and creates new DONJON input. The feedback module in DONJON uses these distributions and the feedback database created with DRAGON to generate new macroscopic cross section in each node. The process is then repeated for the next time-step. In this way, DONJON's kinetic calculation is thus “leading” the calculation over a single time-step, with CATHENA “following” the time-step.

The size of each time-step was determined dynamically by the coupling script based on the maximum power rate of change in each node. If the maximum difference in node power from the previous step exceeded 0.5%, the current step was repeated with a smaller step size. If this criterion was satisfied for multiple consecutive steps the step size was increased. Minimum and maximum time steps of 2.5 × 10⁻⁴ s and 1.0 s, respectively, were used in this work.

A conceptual design for the shutdown system was added to the DONJON diffusion and kinetics model in the form of control blades inserted from the bottom of the core. These control blades, based on Advanced Boiling Water Reactor designs, were initially conceived as part of a reactivity hold-down and channel power shaping system [9]. The cruciform blades are inserted between channels (Fig. 4) and consist of a zircaloy structure with 304 stainless steel neutron absorbers in the blade edges. According to Ref. [9], this configuration reduces the flux/power tilt within the fuel assembly. This is desirable when the blade is inserted during normal operation, as would be the case with a channel power-shaping system. By using these blades as part of a shutdown system, however, it is explicitly assumed that the blades are normally outside of the core and only inserted following a reactor trip. This is suitable for the purposes of this study, which aims to study a shutdown system in isolation of other control systems.

The infinite lattice cell with the control blade inserted was modeled using DRAGON at each axial location for each burnup step and thermal-hydraulic condition. Additional fuel types were then added to the cross section feedback database, which correspond to the “bladed” infinite lattice cell at each axial location. The control blades could then be “inserted” within the DONJON model by dynamically changing the fuel type in each axial plane (from bottom to top) to refer to the “bladed” database directory instead of the reference directory (i.e., without the control blade). The side and corner cells near the heavy water reflector were not modeled with the control blade due to the computational complexity of the geometry, so no blades could be inserted at the edge of the core. Including only the interior locations is nevertheless equivalent to inserting 209 blades in the full core.

A series of static core calculations were performed to determine the reactivity worth of this blade design. The difference in full core $keff$ from the reference condition was calculated for each fraction of insertion (from node 1 to node 20 along the length of the channel), at three initial core states (beginning of cycle (BOC), middle of cycle (MOC), and end of cycle (EOC)), and at both the reference thermal-hydraulic conditions and “voided” conditions. The totally “voided” core in this case is with the outer (fueled) region of each channel filled completely with 20 kg m⁻³ coolant, with no change to the reference coolant density in the center flow tube. The results of these calculations are shown in Fig. 5. The total reactivity worth of the blades is between −44.42 mk and −54.39 mk depending on the core conditions. Variation between BOC and EOC can be reasonably attributed to the evolution of the core power distribution and fuel composition over the length of the cycle. For comparison, the full-core reactivity worth of void around the fuel is between +2.59 mk and +2.35 mk (BOC to EOC), and the delayed neutron fraction used in this study was $β$ = 2.82 mk (other studies have shown minimal variation in $β$ between BOC and EOC, so using a single low value is relatively conservative) [12]. The relative reactivity worth of the shutdown system is thus between 12.62$ and 19.29$, which is more than suitable for this analysis.

In this study, only neutronic trips were considered, eschewing thermal-hydraulic trips such as low pressure or flow, since historically neutronic trips are the fastest responding parameters for LOCAs with positive void reactivity. The “signal” for the shutdown system was thus the sum of node powers in the DONJON model. Seven such signals were considered: the total core power and six separate zone powers, with each zone comprising exactly one sixth of the core. In the 84 channel model, these zones correspond to three separate 28 channel groups (one group for the core interior and two for the core exterior in to equal sectors) split evenly in half along the length of the channel, each relating to 280 nodes in the model. The number of zones and their precise boundaries do not have any justification, but in this way, if the core power were to rise first in only the top or bottom of the core, or only the core interior or exterior, the power increase would still be captured by a trip signal.

where $P$ (W) is the reactor power, $t$ (s) is the time, and $i$ is the current time-step. $Pnom$ refers to the nominal (i.e., initial) reactor power. A reactor trip was initiated if either of these criteria were satisfied for the entire core (i.e., the sum of all 1680 node powers) or any individual zone. Similar trip criteria are implemented in existing reactors with positive void reactivity (e.g., CANDU), although compared to the values that may be used in practice, the trip points used in this study are conservatively high [13].

There are many sources of delay between the time of neutronic trip and the insertion of the shut off rod assemblies. For example, instrumentation, set point comparison, and initiation of the shutdown blade mechanism all have some time response function or fixed delay. Within this work, these responses were summed into a single trip delay, defined as the time interval between which the neutronic-based signal exceeds its limit within the core and when the shutdown blades begin their insertion.

In literature, trip delay times used for CANDU LOCA analysis range from 10.0 ms (the lowest likely achievable) to 100.0 ms (the largest acceptable) [13]. These served as the reference bounding values for this work, with the low delay time applied solely to the NOP trip and the high delay applied solely to the HLR trip to account for the additional processing time. The trip delays were also varied around their reference values (7.5 ms and 15.0 ms for NOP, 75.0 ms and 150.0 ms for HLR) to examine the results sensitivities to these assumptions.

where $IF$ is the rod insertion fraction and $τ$ (s⁻¹) is a time constant. Three different values of $τ$ were used in this work to study the effect of rod insertion rate on the power transient (Fig. 6). The shapes of these curves and the total duration of the insertion are consistent with CANDU's established shutdown system, albeit the devices in this case enter from the bottom of the core rather than the top [15]. While the CANDU devices enter the low-pressure moderator from the top of the core through a combination of springs and gravity, some other mechanism would be required to drive the devices in from the bottom as assumed here for the PT-SCWR. This is not considered to be a significant technical challenge given the low pressure of the moderator, and possible options are not elaborated further in this study. Note also that, implemented as they were in the DONJON model, it was not possible to have partial insertion of a blade within a node, i.e., it must either be completely inserted or withdrawn. The “steps” in the figure thus show the actual modeled insertion within the 20 axial nodes of the core.

Given the variation in rod reactivity worth, trip delay time, and rod insertion rate, it is possible to establish the “performance envelope” of the shutdown system (Fig. 7). The net result of all the above assumptions is a large variation in the rate of negative reactivity insertion following a reactor trip. Studying PT-SCWR transients within this wide range of operation will provide a good indication of a fast-acting shutdown system's efficacy in arresting power transients.

The simulated transient is similar to an inlet LOCA and is driven by the specified core inlet and outlet pressure boundary conditions shown in Fig. 8. These pressure transients are based on models created with the Reactor Excursion and Leak Analysis Program (RELAP) version RELAP5/SCDAPSIM/MOD4, since limitations in CATHENA MOD-3.5d/Rev3 prevented simulation of an actual LOCA transient (i.e., due to technical limitations with this version of the code, the fluid properties could not transition from supercritical water to below 22.1 MPa) [16]. A flow transient was thus imposed where the pressure stayed above 23 MPa for the first 5 s. The drop in the inlet pressure causes the flow through the fuel channels to slow and then reverse, resulting in low density coolant from the outlet plenum traveling into the fuel region of the channel (the slight differences in total core flow between BOC, MOC, and EOC are attributable to the different channel power profiles at the initiation of the transient). The corresponding coolant density decrease around the fuel results in a positive reactivity insertion, raising the core power. A design requirement for the shutdown systems is that a reactor trip should be initiated prior to the coolant conditions dropping below the critical pressure, ensuring that phenomena such as the critical heat flux are precluded while at power. The coupled simulation was terminated 5 s after initiation of the break, allowing sufficient time for a fast-acting shutdown system to operate, while the pressure remained above 22.1 MPa with considerable margin (23 MPa in the imposed LOCA-like transient).

Previously, it was demonstrated that this unmitigated transient would raise the power to >150% of the nominal value before negative feedback brought the core to a low power state [4]. A fast-acting shutdown system, as implemented in this work, should limit the peak and duration of this power pulse. For brevity, bounding cases of “Fast,” “Slow,” and “Intermediate” reactivity insertion, as shown in Table 2, were modeled.

The transient core power for each initial core state and each effective shutdown system speed is shown in Fig. 9. A summary of the trip times is shown in Table 3. Each case tripped on HLR, and the time at which each trip occurred was solely a function of the imposed delay time and not the core state (i.e., BOC, MOC, or EOC). Given that the void worth is approximately the same for each core state, and each used the same value of $β$⁠, it is not surprising that each tripped at the same time. Each trip also occurred on the total core power rather than a specific zone, indicating that in the imposed transient, there were no substantial top/bottom or interior/exterior tilts in the core power distribution.

The peak core power, however, was only determined by the initial core state and not the rod insertion curves. At any reasonable insertion rate, the rods are evidently capable of mitigating the pulse, which is an intuitive result since the reactivity worth of the rods (on the order of 12$ to 19$) vastly exceeds the full-core void reactivity of 2.35 mk to 2.59 mk (Fig. 5). The rods thus effectively terminate the power pulse upon their initial entry into the core region, well before complete insertion.

The total core power is observed not only to peak highest at MOC but also to drop the most rapidly. Interpreting this result requires weighting the positive feedback from the decreasing coolant density against the negative worth of the control rods and the changes in core flux distribution. Previously, it was demonstrated that coolant density reactivity feedback in the fuel region was greater further along the batch cycle (i.e., EOC > MOC > BOC) [4]. All other things being equal, this would suggest that the peak power during this transient should increase during the cycle. However, during the cycle the flux distribution tends to shift from the interior channels to those on the periphery as fuel burnup increases. Given the larger number of rods in the periphery of the core, the shutdown system's effectiveness also increases during a cycle. This is evident by the worth of the rods increasing from 12$ to 19$ as the cycle progresses. The net result of these two competing processes is that the peak power pulse occurs sometime near the MOC conditions.

The peak cladding temperatures in the core (considered for only those channels which had radiation heat transfer modeled, E6, C4, and G8) are shown in Fig. 10. In each case, the peak cladding temperature occurred in the highest power channel, E6. The rapidity of the shutdown has a noticeable but minor effect on the peak cladding temperature. The absolute value of the temperature is demonstrably more correlated with the initial core state than the shutdown system performance. In each case, the temperature continues to slowly rise since the core power is nonzero but little heat is removed by the coolant (heat from radioactive decay was not included in this model, but the fission power has not yet reached zero). The temperature increase is shown to be leveling off by the time the coupled transient is terminated, which is an expected result with radiation heat transfer occurring through the liner tubes, insulator, and pressure tube to the moderator. Previous work has demonstrated the effectiveness of the heat removal postaccident via the moderator system [16].

The power transients shown in Fig. 9 show that, as expected, the case with the slowest shutdown system is most limiting. Axial profiles in the highest power channel, E6, are shown for this limiting transient at BOC, MOC, and EOC in Figs. 11, 12, and 13, respectively. The axial profiles depicted average fuel temperature, coolant density in the fuel region, and coolant density in the center flow tube represent the largest reactivity feedback.

The increase in fuel temperature is shown to be insubstantial and, most notably, the channel power is shown to have decreased to a very low level before significant voiding can occur in the center flow tube. This demonstrates the efficacy of the shutdown system. In the earlier work (transients executed without shutdown rods), only the negative reactivity feedback of center channel void would lower the core power, and hence the net energy deposited (NED) in the fuel would be much larger in cases where the shutdown system was not credited.

The shutdown system implemented in the coupled PT-SCWR core model, consisting of neutron absorbing blades inserted from the bottom of the core, has been demonstrated to be capable of mitigating power excursions initiated by flow transients resembling a large LOCA at the core inlet. This is evidenced by the lower peak core power (131% at EOC, versus >150% in the unmitigated case [4]) and the shortened duration of the excursion (i.e., the core power is lowered before significant negative feedback occurs from center flow tube void). The NED in the fuel is therefore reduced by approximately a factor of two as compared to the unmitigated case.

Since the attributes of this shutdown system closely resemble the fast-acting shutdown system implemented in CANDU reactors (in terms of reactivity worth and insertion rates), it is reasonable to conclude that power excursions in the PT-SCWR can be as effectively arrested as in CANDU large LOCA. Of course, the requirements for a fast-acting shutdown system in the PT-SCWR should be established by an analysis of the fuel and its risk of failure in such a power excursion. Given that the power transient is inherently self-terminating (as a result of coolant void in the center flow tube), and the peak power of the excursion is much lower than a CANDU large LOCA, it is possible that such a fast-acting system may be unnecessary. Nevertheless, this work demonstrates that a fast-acting system is a feasible option to limit peak power, NED, and fuel temperatures during transients in the PT-SCWR.

• Funding to the Canada Gen-IV National Program was provided by Natural Resources Canada through the office of Energy Research and Development, Canadian Nuclear Laboratories (formerly AECL), and the Natural Sciences and Engineering Research Council of Canada (Project NNAPJ 422784-11; Funder ID: 10.13039/501100007178).

Symbols

$IF$ =

shutdown rod insertion fraction

$keff$ =

effective neutron multiplication factor, mk

$LR$ =

log-rate of change in power, s⁻¹

$P$ =

core neutron power, W

$t$ =

time, s

$Tout$ =

channel outlet temperature, °C

Greek Symbols

$β$ =

delayed neutron fraction

$λ$ =

delayed neutron source decay constant, s⁻¹

$τ$ =

shutdown rod insertion time constant, s⁻¹

Subscripts

$i$ =

time step

$j$ =

delayed neutron group

$nom$ =

nominal

Acronyms and Abbreviations

AECL =

Atomic Energy of Canada Limited

BOC =

beginning of cycle

CANDU =

Canada Deuterium Uranium

CATHENA =

Canadian Algorithm for Thermalhydraulic Network Analysis

EOC =

end of cycle

GENHTP =

generalized heat transfer package

HERC =

high efficiency re-entrant channel

HLR =

high log rate

LOCA =

loss of coolant accident

MOC =

middle of cycle

NED =

net energy deposited

NOP =

neutron over power

PHTS =

primary heat transport system

PT-SCWR =

pressure tube supercritical water-cooled reactor

RELAP =

Reactor Leak Analysis Program

SCALE =

standardized computer analysis for licensing evaluation

TRITON =

transport rigor implemented with time-dependent operation for neutronic depletion

wt % =

weight percent