Abstract
In CANDU reactors, shim operation is used when the online refuelling capability becomes temporarily unavailable. Adjuster rods, normally in-core to provide flux flattening, are withdrawn in sequence to provide additional reactivity as the fuel is depleted. In a CANDU 900 reactor, up to three of the eight adjuster banks may be withdrawn, with the power derated accordingly. In this study, the shim operation was modelled using a combination of TRACE_Mac1.1, PARCS_Mac1.1, and scripts modelling the reactor regulating system, all running as a single coupled simulation. A driver script simulated the operation as a sequence of steady-state, depletion, and transient models. The results were compared to operational data from a nuclear power plant, evaluating the key figures of merit. The simulation was extended beyond the operational data by reducing the power to 59% FP and withdrawing the remaining adjusters, to observe the behaviour of the simulated reactor for a deeper reactivity-driven transient. Sensitivity cases, including adjuster rod depletion and nuclear data uncertainty, were also evaluated. This study was able to successfully reproduce the general results of the shim operation. Some discrepancies were observed between the simulation and dataset for the duration of the shim, particularly for the one bank out phase of the shim. Several potential causes for the early phase behaviour were identified. When the simulation was extended, the model predicted that a power reduction below 60% FP would lead to xenon poison out when the adjusters were depleted, with the timing sensitive to the adjuster depletion. Nodalisation of the PARCS model also had a significant impact, due to the effect on adjuster nodalisation and its area-of-effect with respect to the actual adjuster locations. Nuclear data uncertainty had a lesser but still noticeable effect. Other parameters, such as the distribution of fuel burnups in the core, only had a small effect on the shim operation.
1. Introduction
The CANDU reactor has unique characteristics in its operation, primarily due to the use of natural uranium fuel and heavy water moderators. Most notably, CANDU reactors are fuelled online, with the capability to refuel individual fuel channels multiple times per day while at power. This “continuous refuelling” allows for excess reactivity to be kept at a minimum, improving the overall neutron economy to achieve a higher exit burnup out of the natural uranium fuel. CANDU reactors are also loosely spatially coupled, thus requiring spatial reactivity control to maintain a nominal flux distribution with no tilts. When modelling long-term reactor operation, it is important to model the core evolution, including depletion, saturating fission products, and the actions of the reactor regulating system (RRS) [1].
During normal reactor operation, the power of the reactor is controlled primarily using 14 liquid zone compartments. These may be filled and drained with light water, which acts as an absorber in a heavy water-moderated reactor. The RRS liquid zone control algorithm acts to control both the bulk reactor power as well as its spatial distribution and can counteract the reactivity imbalances that result from refuelling along with other minor reactivity perturbations. Adjuster rods, normally inserted, are used to flatten the flux profile in the core and may be withdrawn to provide additional reactivity when necessary. Mechanical control absorbers, normally withdrawn, are used to provide negative reactivity, and may also be partially or fully dropped to the core to rapidly reduce power without fully shutting down the reactor. Moderator poison may also be used to introduce negative reactivity.
Online refuelling is achieved using two fuelling machines, one on either end of the reactor. When refuelling a channel, one fuelling machine will push new fuel bundles into a channel, while the other machine will receive the used bundles. One possible operational occurrence is a fuelling machine outage, where a fuelling machine is unavailable due to maintenance and thus the reactor cannot be refuelled. In this scenario, the reactor may be operated in a “shim operation” mode, using the adjuster rods for reactivity shim. During the initial stage of shim operation, the reactivity loss due to depletion is compensated by draining the liquid zone compartments. When additional reactivity is required beyond what the liquid zone control can provide, banks of adjuster rods are sequentially withdrawn. This is continued until the reactor is either shut down or refuelling resumes.
A number of previous studies have been performed on adjuster rods, shim operation, and/or other adjuster-driven transients under various conditions. Some examples include(i)Evaluating reactivity devices using DUPIC fuel compared to natural uranium fuel [2]. The study modelled CANDU-6 reactivity devices and shim operation using coupled RFSP and NUCIRC, finding that the maximum shim duration was significantly shorter for DUPIC fuel.(ii)Evaluating reactivity devices using DUPIC fuel for xenon override in a CANDU-6 [3], using RFSP, finding comparable xenon override capability to the natural uranium cycle.(iii)Optimisation of adjuster rods in thorium-based fuel cycles to achieve comparable xenon override capability and shim duration to natural uranium fuel cycles [4], using DRAGON and DONJON.
Unlike these studies, which used CANDU-specific tools and evaluated theoretical shim operation performance for advanced fuel cycles when compared to natural uranium, this study focused on implementing the methodology using the PARCS and TRACE codes, which are typically used for the analysis of light water reactors rather than CANDU reactors, as a basis, and modifying them to support the necessary coupling capabilities. This study involved simulating the shim operation of a 900 MW class CANDU reactor using a standard natural uranium fuel cycle and comparing with operational data, to demonstrate that these tools could be adapted for modelling long CANDU transients and considered for analysis of either actual or theoretical reactor transients. In addition, the use of coupled transient simulations to supplement quasi-steady-state calculations is a novel contribution, as it allows for the reactor dynamics to be properly captured at key points in an otherwise quasi-static reactor evolution.
The simulation and analysis of reactor operation may be performed using one of two methodologies based on the evolution of the reactor state. If the reactor fuel composition and transients evolve very slowly, then the state of the reactor at any given time can be modelled as a quasi-steady state. Under such a simulation, the burnup and saturating fission product state of the core, which evolves slowly, is held constant, while more rapidly evolving variables such as the neutron flux and thermal-hydraulic conditions are solved for their equilibrium values. The flux is then held constant and used to advance the burnup and saturating fission product state by a given period of time. If the core state evolves rapidly, the reactor is simulated as a transient, with the neutron flux and thermal-hydraulic conditions evolving with time, and the effects of fuel composition changes are secondary.
2. Materials and Methods
2.1. Background
The CANDU reactor is a pressurised heavy water reactor (PHWR) design, differing significantly from the more popular light water reactor (LWR) designs. Some key traits of the CANDU design that differ from the PWR design include [5](i)Heavy water is used as the coolant and moderator(ii)Natural (unenriched) uranium as fuel(iii)Horizontal fuel channel arrangement, with fuel channels serving as the pressure boundary. Coolant and moderator are physically separated, with the fuel channels surrounded by the moderator in an unpressurised calandria.(iv)Fuel channels have a significant degree of separation, to allow for moderation by heavy water, access by the fuelling machines, and the insertion of vertically oriented reactivity devices(v)Online refuelling, through the use of two fuelling machines which can latch onto either end of a fuel channel to insert and discharge fuel bundles while the reactor is operating at power(vi)Reactivity control primarily uses 14 liquid zone compartments, whose light water levels can be varied independently to control reactor power and its spatial distribution
Online refuelling is important due to the use of natural uranium fuel. This fuel has far less excess reactivity when compared to enriched uranium fuel, leading to lower exit burnups and core residence times for CANDU fuel when compared to LWR fuel. A batch refuelling scheme for a CANDU reactor would require a refuelling outage every few months [5], which would greatly diminish the station’s capacity factor. As well, the continuous refuelling allows for the excess reactivity on a core-wide basis to be kept at a consistent minimal value, improving the overall neutron economy and thus the exit burnup.
Online refuelling is achieved using a pair of fuelling machines, with one on each end of the reactor. Each fuelling machine is capable of latching onto one end of a fuel channel, removing the channel closure and shield plug, then either inserting new fuel bundles or receiving used fuel bundles, before replacing the shield plug and channel closure [5]. Typically, not all fuel bundles in a channel are replaced, but rather only 4 or 8 bundles (out of 12 or 13 bundles) are replaced, with the remaining bundles closer to the inlet end of the channel remaining in the core but being shifted towards the outlet end of the channel. The horizontal orientation of the fuel channels allows for bidirectional refuelling, where adjacent channels are refuelled in opposite directions, along with coolant flow being in opposite directions, to achieve symmetry in the neutron flux distribution.
This introduces additional complexities in simulating reactor operation which is not present in the simulation of LWR operation. For an LWR, the fuel's physical location remains the same over the course of a fuel cycle, with the fuel being depleted over the course of the cycle and the core kept critical through the use of moderator poison, reactivity devices, and/or negative power feedbacks, depending on the type of reactor. In a CANDU reactor, the fuel distribution changes on a daily basis as individual channels are refuelled. Thus, high-fidelity simulations require a full accounting for both effects of burnup and fuel relocation due to fuelling prior to the event of interest, as well as the transient of interest. Different methods for performing the core follow are as follows:(i)Using actual station data, the operating history of the reactor is simulated, from the start of operations or from a given snapshot in time to a later snapshot in time using the fuelling and operating histories of that core. This may be performed using any diffusion code, with the industry typically using RFSP, but studies have also been done using DONJON [6], to model the core state leading up to a transient.(ii)Starting with an initial core of fresh fuel, reactor operation is simulated, including the simulation of fuel channel selection for refuelling, until the core reaches an equilibrium. This typically takes 400 to 500 full-power days [5]. Such fuelling simulations have been previously performed using PARCS [7].
When simulating reactor operation, it is thus often important to have frequent snapshots which capture the salient core features at each discrete point in time. For a given snapshot, the slow-changing quantities, such as burnup, may be held constant, and the equilibrium state of the reactor is determined. Then, a depletion step may be performed, where the slow-changing quantities are updated under the assumption that the equilibrium state of the reactor changes little between successive snapshots. The burnups in one or more channels would then be updated to simulate refuelling. Depending on the model and particularly on the time between snapshots, the saturating fission product densities may be treated as a slow-changing quantity to calculate transiently, or as a fast-changing quantity to calculate as an equilibrium. For Xe-135, the density reaches equilibrium on the order of a day for a reactor at full power, based on an approximately 6.7-hour half-life of I-135 [5]. Thus, for daily or less frequent snapshots, an equilibrium xenon model is reasonable except for recently refuelled channels.
In addition to determining the flux distribution in a core snapshot, it is also important to model the reactivity devices, since reactivity devices have an important effect on the flux distribution. The following are used to control reactivity in the CANDU design:(1)14 liquid zone compartments (LZCs) distributed throughout the core. These are filled with light water, which acts as a neutron poison due to its high absorption cross section compared to heavy water. The unfilled portion of each compartment is filled with helium gas. The liquid levels are varied automatically by the reactor regulating system (RRS) running on a digital control computer (DCC). Both bulk control and spatial control are used. Bulk control adjusts the liquid levels equally in all 14 compartments to add positive or negative reactivity and control the reactor power. Spatial control adjusts the liquid levels independently to balance the spatial distribution of the flux [1].(2)Adjuster absorbers, which are normally inserted in the inner part of the core to provide flux flattening. Flattening of the flux distribution is desirable as it allows for the total reactor power to be increased without exceeding power limits for individual fuel channels or fuel bundles. Adjuster absorbers are grouped into eight banks. The RRS will withdraw one or more banks of adjusters if additional positive reactivity is required, based on either a low average liquid zone level (indicating that the LZCs cannot provide enough reactivity) or a large negative power error (indicating that positive reactivity is needed more rapidly) [1]. Withdrawing some of the adjusters requires the reactor power to be derated due to a loss of the flux flattening effect.(3)Mechanical control absorbers (MCAs), which are normally positioned out of the core. The RRS will insert one or more banks of MCAs if additional negative reactivity is required, based on either a high average liquid zone level (indicating that the LZCs cannot provide enough reactivity) or a large positive power error (indicating that negative reactivity is needed more rapidly) [1].(4)Moderator poison, either boron or gadolinium. Boron is used on initial start-up to depress the excess core reactivity from a full load of fresh fuel during the first weeks of operation. Gadolinium is used to compensate for transient excess positive reactivity beyond the worth of the LZCs and MCAs, particularly to manage the xenon transient resulting from an outage and subsequent restart [5].
For managing the normal operation of the core, the LZCs are most important, with the other reactivity device states typically not being actuated/changed. The LZCs are used to manage small reactivity perturbations, particularly those that occur from depletion and refuelling. When a channel is refuelled, not only does the new fuel introduce positive reactivity, but it also locally increases the flux near the refuelled channel. For a large reactor such as the CANDU-9 reactor, the core radius is far greater than the neutron migration length, and the reactor is susceptible to certain flux harmonics, particularly the first azimuthal (side-to-side or top-bottom) and first axial (end-to-end) modes [8]. Therefore, spatial control, via the LZCs, is required to control these flux modes. Conceptually, the reactor is divided into 14 spatial control zones, with one LZC corresponding to each zone. Spatial control is used to maintain a design ratio of measured powers between the 14 zones, suppressing any tilts resulting from refuelling, xenon perturbations, or the movements of other reactivity devices. When calculating a snapshot, the 14 LZC water levels must be determined, both to achieve criticality as well as to eliminate any flux tilts.
It is also of interest to simulate operating conditions other than normal operation at full power. One such example is the simulation of a shim operation. CANDU online refuelling relies on a pair of fuelling machines. If the fuelling machines are not available, then the reactor will need to be operated without refuelling. In shim operation, the negative reactivity due to depletion is counteracted by the LZCs and the adjusters. Reactor power must be derated prior to withdrawing adjuster banks to counteract the resulting increase in channel and bundle power peaking factors. These power reductions also lead to xenon perturbations that must be managed by the LZC system. Once the fuelling machines are available, the reactor can gradually be refuelled until it reaches its normal level of reactivity. If the fuelling machines are unavailable for an extended period of time, the reactor may need to be shut down. In the case of planned maintenance, extra refuelling can be performed prior to shim operation, using LZCs or moderator poison to compensate for positive reactivity. This will extend the possible duration of the shim operation.
Figure 1 provides the liquid zone designation for the CANDU 900 reactor used in this work.
2.2. Methodology for Modelling and Simulation
This work uses the TRACE [9] and PARCS [10] computer codes, licensed through the NRC’s Code Application and Maintenance Program (CAMP), to simulate system thermal-hydraulics and core physics, respectively. More specifically, this work uses modified versions of the codes, hereafter named TRACE_Mac1.1 and PARCS_Mac1.1. It also uses the Exterior Communications Interface (ECI) [11] to couple an RRS model to PARCS through TRACE. The RRS model was developed in the Python programming language and uses a Python ECI library built around the original Fortran-based ECI library. This coupled model has been previously evaluated against shorter transient events such as the figure of eight flow oscillations (using a modified system model without header interconnects) and a loss of flow transient initiated by a loss of Class IV power [12]. This section covers the additional work performed to model long-duration transients such as adjuster shim operation.
For modelling of CANDU adjuster shim operation, there are two simulation regimes to consider. During the majority of the operation, the core state evolves slowly, and the only significant effects to consider are the effects of depletion and saturating fission product evolution. The second regime corresponds to a time-dependent situation where device configuration (e.g., adjuster positions and liquid zone levels) change over a relatively short period of time due to a change in RRS control actions.
During the slowly evolving quasi-steady-state periods the reactor can be modelled as a steady-state simulation. During this simulation, the depletion, saturating fission product densities, and reactor power are held at constant values, while the spatial flux distribution, thermal-hydraulic state, and control system state are computed. To model the core evolution during these steady-state steps, a PARCS depletion calculation is performed with PARCS_Mac1.1, with the reactivity device and thermal-hydraulic feedbacks held constant. The depletion calculation is used to calculate an updated depletion and saturating fission product density state for the next steady-state step. This methodology is very similar to a core follow, as it generates periodic core snapshots modelling the evolution of the core over a period of time. Core follow analysis has been previously performed using both PARCS [7] as well as other core physics tools, particularly RFSP.
For the steady-state model, certain simulation parameters have been adapted to support eigenvalue calculations, accelerate model convergence, and improve performance.(i)The power error calculation module (CEP) may override the power error with either a keff error or a zone level error for steady-state analysis, to converge either the system keff or average zone level to a target value. The zone levels may also be frozen to a set of known values.(ii)The PARCS eigenvalue calculation is performed only once every 20 TRACE time steps. This greatly accelerates the simulation as the PARCS eigenvalue calculation is significantly more computationally intensive than the TRACE time step. Sensitivity studies demonstrate that the transients predicted here are not sensitive to the exchange rate in this regard, as long as the three models (neutronics, thermal-hydraulics, RRS) are given sufficient time to converge and that the neutronics-RRS feedback loop is stable.(iii)The RRS time step size is decoupled from the TRACE simulation and set to a fixed constant. This allows for the RRS to converge more quickly when the TRACE time step size is small. The maximum step size is limited by stability due to the feedback between the RRS and PARCS models and is currently set to 4 seconds per PARCS step (0.2 seconds per TRACE time step).(iv)Convergence of thermal calibration of the ion chambers and in-core flux detectors (ICFDs) is greatly accelerated. Calibration is applied every two seconds instead of every 16 seconds, and the time constants of the filters in the thermal calibration routine are decreased to 1.875 seconds (instead of 30 seconds or 180 seconds).(v)Calculation of the delayed components in the ICFD response is performed using the flux from the steady-state calculation of the previous depletion step and the flux from the current calculation, under the assumption that the flux changed linearly over the course of the previous depletion step. The ICFD detector response in this work is based on the response characteristics calculated in reference [13], with the four normalised cases averaged. The resulting prompt fraction was 1.048 (4.8% overprompt), and the response also includes four different delayed components with time constants ranging from minutes to days.
This approach is very similar to RFSP’s asymptotic zone control functionality [2], except that the liquid zone control module functions identically in the steady-state and transient modes, only manipulating its inputs (bulk power error and time step size) to implement the steady-state mode. This guarantees that the asymptotic zone levels are consistent with transient modelling, as long as the power error converges to zero.
The second regime corresponds to time-dependent actuation of the device positions caused by RRS or operator actions. As the reactor is depleted, the reactor regulation system is expected to decrease the average liquid zone level in order to maintain criticality. Once the liquid zone level reaches a certain threshold, a bank of adjuster rods is to be withdrawn. The dynamic nature of the adjuster rod withdrawal process is simulated using a coupled transient model with transient depletion and saturating fission product evolution. In advance of withdrawing the adjuster bank, the reactor power will be derated by the operators to an appropriate level for withdrawing the adjuster bank consistent with station operating rules. This is done at a higher liquid zone level threshold, to allow the xenon-135 transient to trigger the adjuster withdrawal and provide enough reactivity margin so that only one adjuster bank is withdrawn. The full transient procedure is thus as follows:(1)As the zone level nears a threshold where action to remove the adjuster may be required, the simulations will decrease reactor power to the setpoint for the next adjuster bank to be withdrawn to mimic operator actions.(2)The simulations will then continue in time to allow the average LZC level to decrease to the adjuster withdraw threshold as the xenon transient progresses.(3)The simulations will then withdraw an adjuster bank at the time and rate dictated by the design parameters. The average LZC level will increase to maintain criticality.(4)Average LZC level will continue to decrease to a minimum value corresponding to peak xenon, then subsequently increase as xenon density returns to a new equilibrium value.(5)As xenon approaches equilibrium, the average LZC level will begin to decrease again once depletion becomes the dominant effect.
For this work, the coupled transient is simulated from the start of the power derate until shortly after each adjuster bank is withdrawn. The most rapid changes in flux distribution and reactivity device positions occur during adjuster bank movement, and it is desired to model the behaviour of the flux distribution as the bank moves out of the core. A significant transient flux tilt may degrade safety margins, and if sufficient flux tilt is detected, an RRS setback may be triggered, further reducing power below the intended setpoint. In this model, the setback threshold for flux tilt is set to 20%, based on zone power measurements excluding the centre zones 4/11. Once a setback is initiated, it continues until the flux tilt drops below 18%. Once a bank is fully out of the core, the spatial flux distribution is relatively stable, and the depletion and saturating fission product evolution may be modelled using a steady-state depletion model once again.
Therefore, the overall procedure for the shim simulation is as follows:(1)Perform steady-state simulation and depletion steps while monitoring the average zone level (AZL).(2)As AZL approaches the threshold to derate power, determine the depletion step size to reach the AZL threshold. Iterate until the calculated AZL is at the threshold.(3)Perform coupled transient simulation of the power derate and subsequent adjuster bank pull.(4)Perform steady-state simulation and depletion steps until the end of the xenon transient.(5)Repeat steps 1–3 for the next adjuster bank.
For step 1, a 6-hour depletion interval has been selected. Once step 3 has been completed, step 4 is performed until 12 hours pass from the end of the previous adjuster bank withdrawal before returning to step 1. The step size starts at 15 minutes and is doubled until it reaches the original 6-hour interval again. The 12-hour delay prevents the next adjuster bank withdrawal from being triggered by the AZL threshold due to the xenon transient, ensuring that each adjuster bank movement is triggered due to depletion and not xenon evolution.
Figure 2 outlines a high-level methodology that can be used to implement this simulation. The simulation alternates between performing a coupled steady-state calculation and a standalone depletion calculation, until the average zone level approaches the threshold to perform the adjuster bank withdrawal process. At this stage, a coupled transient run is performed to simulate the power manoeuvre and subsequent adjuster bank withdrawal.
A more detailed version of the methodology is provided in Figure 3. While similar to the high-level methodology presented in Figure 2, it includes more detailed logic for determining which action to perform after a steady-state simulation. First, it checks that there were no abnormalities in the simulation, such as a setback, stepback, or reactor trip. If such an abnormality is detected, then the simulation will either revert to before the previous depletion step and repeat it with a smaller depletion interval or will switch to coupled transient simulation. This allows for the simulation to pinpoint the occurrence of such an event and simulate it as a transient simulation. Secondly, the PARCS steady-state simulation is limited to a fixed reactor power, but the indicated power in RRS does not exactly match the power in the PARCS simulation. Therefore, at the end of a steady-state step, the indicated power in the simulated RRS is compared to its setpoint, and if their difference exceeds a certain tolerance, the PARCS power is adjusted, and the steady-state step is repeated.
The simulation must then determine whether it has reached a part of the simulation that should be performed as a fully coupled transient. This occurs when either of the following conditions is met:(1)The average zone level is at or below the threshold to perform the next manoeuvre, which consists of a power derate (if needed) followed by an adjuster pull.(2)The average zone level is above the threshold, but the estimated time until the average zone level reaches the threshold is less than the minimum depletion step size (864 seconds, or 0.01 days). This estimate is calculated by extrapolating the average zone level from the previous two steady-state calculations. If the estimate is negative, due to the average zone level increasing over the previous depletion step, then it is ignored.
If neither condition is met, then the next step performed is a depletion step. A nominal depletion step size of 6 hours was selected in this work, but the size of a given depletion step is determined as follows:(1)If the estimated time from step 2 of the previous list is positive and no greater than the nominal depletion step size plus 15 minutes, then the depletion step size is set to this estimated value. Otherwise, the nominal step size is selected.(2)In addition, for a scheduled action that is not to be simulated using a coupled transient (e.g., refuelling), the depletion intervals are sequentially shortened as the action is approached, up until the action time is reached.(3)The actual depletion step size is the lesser of this selected step size or twice the step size of the previous depletion step.
This creates the effect of performing fixed depletion intervals until the next transient portion of the simulation approaches, then performing a depletion interval which is estimated to bring the simulation to the starting time of the transient. After a transient is completed, the subsequent depletion step size is set to 1 hour, then allowed to double with each depletion step until it reaches the nominal step size. For actions not requiring transient simulation and that are based on a schedule, finer steps are taken around the time of the action, both before and after. For this transient, several refuelling events are scheduled. To simulate a refuelling, a script is used to read the binary depletion file produced after a depletion step and modify the burnups in the selected channel (pushing 4 or 8 new bundles into the channel and pushing the existing bundles towards the end of the channel, discharging channel out the other end) before rewriting the file.
Since exact controller parameters for the 900 MW CANDU design are not available in the literature, realistic approximate or representative values were selected. The following values were used:(i)KP = 31.0 (LZC bulk controller gain)(ii)KT = 5.0 (LZC spatial controller gain for flux)(iii)KH = 0.6 (LZC spatial controller gain for level, when flux control is phased out)(iv)KL = 0 (LZC spatial controller gain for level, independent of flux control being phased in)(v)Adjuster movement speed = 10.5 cm/s(vi)LZC drain rate = 5.0 cm/s(vii)LZC max net fill rate = 5.0 cm/s(viii)LZC valve lift/bias for constant level = 60%(ix)Flux tilt setback threshold = 20% FP (difference from highest to lowest zone power, excludes centre zones 4 and 11)
The LZC drain/fill rate was selected based on taking approximately one minute to fill/drain a “long” zone compartment between 0% and 100% [1], the longest of which in this model is 285.8 cm, resulting in a rate of 4.76 cm/s, which was rounded to 5 cm/s.
The most sensitive parameters are the gains and detailed sensitivity results are performed to determine the effect of the gains on the transient. Outcomes related to reactor setback are also highly sensitive to the flux tilt setback threshold.
2.3. Station Description and Modelling
The station design being simulated is a 900 MW class CANDU reactor, based on those used in previous coupled simulations of TRACE_Mac1.0 and PARCS_Mac1.0 [12]. When comparing preliminary results to station data, it was found necessary to update the core physics model for better consistency with the specific CANDU 900 under analysis. This included increasing the number of fuel bundles per channel from 12 to 13 both in the thermal-hydraulic and reactor physics models, as well as applying appropriate extrapolation distances as these are known to have a significant impact on full-core modelling outcomes [14] and was found to have a significant impact in this work.
In addition, adjuster rod ageing was accounted for by prorating adjuster incremental cross-sections, based on prior industry investigations, to account for their decrease in worth as a function of irradiation. An adjuster rod age of 182,500 effective full-power hours (EFPH) was selected for this work. The effect of ageing was to reduce the effective thermal absorption cross-section of the adjuster by approximately 10% to 20%, depending on the adjuster type. Other cross-sections were perturbed by varying amounts.
The RRS model was also updated, particularly to implement an approximate model for fully instrumented channels (FINCHes). The FINCH powers are approximated as the PARCS computed channel powers with a 30-second decay filter. The FINCH powers are compared against a nominal power distribution and are used to calibrate the zone power readings from the ICFDs as part of the RRS control features modelled in these simulations.
In the TRACE model, the 480 channels are reduced down to 40 representative channels, based on zone designation, then subdividing each zone into two regions, being the upper and lower half for the centre zone, and inner and outer regions based on average channel power for other zones. All channels in a given region that belong to the same flow loop and pass are grouped into a single channel in the TRACE model. This results in a total of 40 channels in the TRACE model to represent the flow through all 480 channels. The coolant temperature and density, as well as fuel temperature, are shared for all channels in a given grouping.
The 900 MW CANDU design uses flow-power matching to achieve a similar exit enthalpy in channels with varying average channel power, trimming the flow in lower-power channels [5]. To approximate this flow-power matching, a minor loss representing a flow orifice was added to the lower-power channel groups where the inlet feeder connects to the inlet end fitting such that the outlet enthalpy of all channels was approximately equivalent.
For most CANDU reactors, adjuster bank withdrawal typically requires power reductions to ensure thermal margins are maintained during these periods where radial core power peaking is increased. To simulate this, the simulations perform the following power changes with each adjuster bank:(i)One bank out: 93.5% FP(ii)Two or three banks out: 90.0% FP(iii)More than three banks out: 59.0% FP
2.4. Event Timeline
The transient being simulated is a period of shim operation where only intermittent refuelling was available. A summary of events is as follows:(i)Day 1: initial core state. The reactor power setpoint is at 100% FP. Reactivity banking was performed in anticipation of the shim operation. A small amount of moderator poison (modelled as equivalent boron) is used to compensate for the excess reactivity.(ii)Day 5–6: seven channels are refuelled (32 bundles total).(iii)Day 13: reactor power setpoint is decreased from 100% FP to 93.5% FP in accordance with shim operation procedures. Adjuster bank A is subsequently withdrawn.(iv)Day 18: reactor power setpoint is decreased from 93.5% FP to 90.0% FP in accordance with shim operation procedures. Adjuster bank B is subsequently withdrawn.(v)Day 20: there is a brief period where the fuelling machine is available. Four channels are refuelled (16 bundles) to extend the shim operation.(vi)Day 25: adjuster bank C is withdrawn.(vii)Day 29: end state for shim operation.
Figure 4 presents the timeline for reactor power and liquid zone level changes. The depletion trend, three adjuster bank manoeuvres, and the refuelling run are visible in the liquid zone level trend.
To simulate this sequence of events, specific conditions must be identified for triggering the initiation of each event. The simulated sequence is as follows:(i)The simulation starts at the initial known burnup distribution, with 0.35 ppm of boron used to represent the equivalent moderator poison that was deemed to be present. The moderator poison is gradually reduced to zero, following the known trend of equivalent moderator poison from station snapshots and linearly interpolating between them.(ii)The seven-channel refuellings on days 5 and 6 are performed consistent with station actions.(iii)At an average LZC level of 35%, the reactor power setpoint is set to 93.5% FP. The subsequent xenon transient causes the average LZC to decrease rapidly. At an average LZC level of 25%, adjuster bank A is withdrawn. The xenon transient is allowed to continue until reaching a new equilibrium.(iv)At an average LZC level of 30%, the reactor power setpoint is set to 90.0% FP. This causes a similar transient as before, and at an average LZC level of 25%, adjuster bank B is withdrawn.(v)While the “shim” operation typically does not have fresh fuel inserted, in this operating transient there were 4 fuelling operations. For the simulation, the first channel refuel is performed 1.9 days after the withdrawal of adjuster bank B is completed. This timing and the interval between refuelling each of the four channels are taken from the operational history of the station.(vi)At an average LZC level of 25%, adjuster bank C is withdrawn, increasing the zone level.(vii)The average liquid zone level continues to decrease from depletion until it reaches a final zone level of approximately 25%, at which point the comparison with the station data ends.(viii)The reactor power setpoint is set to 59.0% FP, as outlined in the previous subsection, to extend the simulation and evaluate the effect of such a power reduction. This causes a xenon transient that causes the average LZC level to rapidly decrease. The subsequent adjuster bank pulls to compensate are performed each time the average liquid zone level gets too low.(ix)With all adjusters out, if the average zone level drops to 10%, the reactor is shut down, as it would go subcritical anyway. The simulation is ended when the actual reactor power has dropped below 2.0% FP or when the simulation time has exceeded 6 hours from the power reduction to 59% FP.
2.5. Figures of Merit
The following key figures were identified for evaluating the performance of the simulations compared to the station data:(i)Event timings (adjuster bank A, B, C withdrawals and subsequently reaching 25% average zone level)(ii)Trends in average liquid zone levels(iii)Trends in individual liquid zone levels(iv)Channel power distribution
In addition, it is possible to calculate reactivity device worths in the simulated cases for comparative study.
For individual liquid zone level trends and channel power distribution, the average zone level of the two central zones (4 and 11) was selected as an indicator of the radial reactivity balance. More specifically, it was checked 5 minutes prior to the end of the station data. This point was selected as the average zone level would be consistent (around 25%) and the level for the two central zones was generally within the range of 30% to 70% for both the station data and the simulations.
The following phenomena have a significant effect on these figures of merit as they affect the core reactivity:(i)Depletion–the rate of reactivity change due to depletion will govern the overall length of the transient.(ii)Power coefficient of reactivity–a positive coefficient will result in a reactivity decrease when the power is derated for an adjuster withdrawal, shortening the time to the next subsequent adjuster movement. A negative power coefficient has the opposite effect.(iii)Void coefficient of reactivity–the CANDU reactor has a positive void coefficient of reactivity. Its effect on the shim transient is similar to the power coefficient, except only at high power, where voiding would be expected near the channel outlets for high-power channels.(iv)Absorption by reactivity devices–each reactivity device affects the core both through its absorption of neutrons as well as the resulting redistribution of flux in the core.(v)Saturating fission products–when the power is decreased, there is an initial negative reactivity change as the concentration of Xenon-135 increases. After a couple of days, a new equilibrium is approached with a slight net positive change in reactivity. The xenon-135 negative reactivity transient is also the driving force of the extended portion of the simulation when the power is reduced to 59% FP.
The duration of each phase of the shim operation (one bank out, two banks out, and three banks out) is determined by the reactivity worth of the adjusters. For one and two banks out, there is also an effect from the power coefficient, void coefficient, and saturating fission products, due to the power reduction immediately prior to the removal of the first two adjuster banks. The rate of reactivity change due to depletion also influences the duration of the shim operation.
For the two-bank-out portion of the shim operation, four channels are refuelled, replacing 16 depleted bundles with 16 fresh bundles. The resulting reactivity increases the delay in the withdrawal of bank C and subsequent events of the shim operation.
The outcome of the 59% FP portion of the simulation is sensitive to the reactivity worth of the remaining five adjuster banks as well as the xenon-135 transient, along with some contribution from the power coefficient of reactivity. The power coefficient affects the initial part of the transient–the first seconds to minutes after the power manoeuvre as the TH conditions settle to a new equilibrium. Beyond that, the core reactivity is dictated predominantly by the xenon-135 transient, and adjuster rod banks are withdrawn to compensate for this transient.
The average liquid zone level may be used as an approximate analogue for the reactivity balance of the core. The liquid zones contribute roughly 6 mk of reactivity. Other changes to core conditions that affect reactivity will introduce a compensating change to liquid zone levels from the RRS. However, two considerations must be kept in mind when considering the liquid zone levels as a reactivity measure.(1)The reactivity to zone level relationship depends on which zones are being filled and their levels. For example, the top of zones 3 and 10 are near the edge of the core and extend into the reflector, thus the marginal reactivity worth of these zones is low when their fill level is high.(2)The model may overestimate or underestimate the reactivity worth of the liquid zones, thus underestimating or overestimating the change in zone level resulting from a given reactivity change.
Thus, the analysis of any discrepancies in the average zone level should consider the entire simulation to determine whether a discrepancy is due to the non-LZC reactivity change or due to the LZCs themselves. For example, if the LZC level change is consistently underestimated throughout the transient, it suggests that the worth of the LZCs themselves is being overestimated.
The individual liquid zone levels may be used to evaluate the spatial distribution of reactivity, due to the distribution of fuel burnups, adjuster rods, and xenon-135. The LZCs are used to spatially control the flux distribution to compensate for perturbations due to refuelling and xenon-135 and will also attempt to compensate for adjuster movements. Any differences in individual zone levels from the data suggest a reactivity discrepancy for part of the core.
2.6. Sensitivity and Uncertainty Analysis
For this study, the sensitivities of the shim operation to different uncertainties and parameters must be identified. These may be compared to any differences identified between the model and station data to determine any potential sources of discrepancies. The sensitivities may also be evaluated by comparing the different models to one another. Some of the key parameters being evaluated are as follows.
2.6.1. Burnup Distribution
To evaluate the effect of the distribution of fuel burnups, a simulated core follow calculation was performed on the reactor. Using a random snapshot as a starting point, thousands of days of core operation and refuelling are simulated prior to the shim transient. This generates a set of burnup distributions that can be used as initial conditions for simulating the shim operation. The results from different initial core snapshots may be used to identify whether the fuel burnup distribution can potentially affect the evolution of the transient. If different snapshots behave similarly during shim operation, it suggests that the actual station burnup distribution was unlikely to have a significant effect.
2.6.2. Nuclear Data
The neutronic calculation is highly reliant on the nuclear data used in the lattice physics code (SCALE/TRITON [15]) to calculate the homogenised cross-sections for the diffusion code (PARCS [10]). SCALE includes a stochastic uncertainty analysis tool, called Sampler, which runs a SCALE sequence using one or more perturbed sets of nuclear data. These data sets are preperturbed using the covariance matrix that accompanies the nuclear data and are distributed with SCALE.
For this work, lattice calculations were performed using 60 perturbed nuclear datasets from the SCALE 252-group library based on ENDF/B-VII.1. The depletion and branch calculations are then performed using each dataset, for both the infinite lattice cell and the lattice edge cell, along with performing reflector calculations. The homogenised fuel and reflector cross-sections are then combined with the device incremental from the supercell calculations, which are not perturbed by Sampler.
This creates a set of perturbed PARCS PMAXS cross-section datasets which may be used either to simulate the station transient from its actual burnup state or to perform simulated core follows to generate snapshots to simulate the station transient from different burnup states. A new core follow analysis is performed for each dataset to construct core snapshots consistent with the perturbed fuel properties, as described in [7].
It should be noted that the Shift module of SCALE [16] was not released and hence the nuclear data uncertainty impact on the device incrementals could not be directly calculated. Instead, incremental cross-sections were perturbed directly as described in the paper.
2.6.3. Incremental Modelling and Uncertainty
While it is not possible to perform a proper uncertainty analysis for the reactivity devices in the current methodology, it is possible to perform a sensitivity study. The device incrementals can be perturbed in order to increase or decrease their worth. There are several uncertainties in the modelling of reactivity devices.(i)Uncertainties in the true device geometry due to manufacturing tolerances (assumed to be small).(ii)Uncertainties in the lattice supercell calculation, either from the microscopic cross-section data or from modelling approximations.(iii)Modelling uncertainties introduced by the application of the device incremental cross-sections to the PARCS core model, e.g., from the expansion of the device’s effective volume to fit to the PARCS mesh.(iv)Adjuster rod irradiation.
2.6.4. Incremental Cross-Section Placement in the Lattice
The lattice supercell calculation uses a repeating lattice of 2 lattice pitches across by 1 lattice pitch vertically by 1 bundle length axially (2 × 1 × 1 region) and homogenises a 1 × 1 × 1 region centred on the reactivity device. Incremental cross-sections are determined by comparing the macroscopic cross-sections with and without the device. It is assumed that this supercell geometry is representative of all possible device locations in the core, and the same incremental cross-sections are applied for a given device for all PARCS nodes, even reflector nodes, except that fission incrementals are zeroed for the reflector nodes.
While the homogenisation region is a 1 × 1 × 1 region centred on the reactivity device, since the device is vertically oriented (Y-direction in PARCS), running perpendicular to and between the fuel channels (Z-direction in PARCS), it is not possible to apply the incremental cross-sections directly to a 1 × 1 fuel channel region. Instead, the incremental cross-sections are distributed to both adjacent fuel channels, with a weight of 0.5 assigned to each channel. Furthermore, as most reactivity devices are not perfectly centred on a single axial node, the incremental cross-sections must be further distributed between two axial planes, with the relative weights determined by the actual device location. The distribution of device incrementals over a larger volume has a tendency to weaken the spatial self-shielding effect, increasing their predicted reactivity worth.
The following three models are evaluated by modifying the axial (Z-direction, parallel to the horizontal fuel channels) nodalisation:(i)12 axial plane model, with 12 fuel bundles per channel, and extrapolation distance applied axially. The burnups from the station, which has 13 fuel bundles per channel, are averaged between adjacent bundles. In this model, the liquid zone compartments are centred on an axial plane, while the centre row of adjusters is distributed between the two centre axial planes.(ii)13 axial plane model, with 13 fuel bundles per channel, with the lengths of the two end planes being equal to a ½ bundle plus an extrapolation length (while the actual length of the core is exactly 12 bundle lengths, the flux does not go exactly to zero at the ends of the core, with the extrapolated length being approximately 12 cm longer than the actual length (see page 22 of chapter 2 of [5]). For the 13 plane model, this extrapolation is applied by lengthening the nodes at either end). In this model, the centre row of adjusters is centred on an axial plane, while the liquid zone compartments are distributed between two axial planes.(iii)26 axial plane model, with 13 fuel bundles per channel, with the length of the two end planes consisting solely of the extrapolation length used in the 13 plane model. Both the liquid zone compartments and the adjusters are distributed between two axial planes, but as the length of each axial plane is half of its value in the other cases, this is equivalent to one axial plane in the other cases.
The reference case uses the 13 axial plane model. It should be noted that, while the modelling of the centre row of adjusters is significantly affected by the axial nodalisation, the effect on the other two rows is less as they do not fall neatly along a bundle or half-bundle position axially (they are offset 80 cm from the centre row on either side or 1.62 bundle lengths). For these adjusters, the 12 axial plane model has these adjusters more heavily weighted towards one axial plane. Thus, the effect of the nodalisation of these adjusters is expected to be similar to the liquid zone compartments, except with a smaller magnitude.
2.6.5. Adjuster Irradiation/Ageing
As the adjuster rods are normally in-core, they are continuously exposed to the full in-core flux and are thus irradiated over their lifetime. The neutron-absorbing materials in the rods will gradually be transmuted to less-strongly absorbing isotopes, over time decreasing the macroscopic absorption cross-section of the rods. The effect can be significant over the lifetime of the adjuster rods. For the CANDU 900 data, the adjuster rods are significantly aged and therefore ageing factors must be applied to achieve accurate results for keff and flux distribution. The incremental cross-sections for unirradiated adjusters, calculated using Serpent, are perturbed by an ageing factor dependent on the adjuster rod age, adjuster rod type, and cross-section type. These aging factors were derived from prior industry work using DRAGON. Perturbation factors are calculated separately for each reaction type and energy group.
For this work, an adjuster rod age of 182,500 effective full-power hours (EFPH) was selected. Adjuster rod ages of 195,000 EFPH and 154,008 EFPH were also selected as sensitivity cases, along with a “new adjusters” case.
3. Results
3.1. Reference Case
For the reference case, using a nodalisation of 13 axial planes (one bundle per plane, with the end bundles shortened to a ½ bundle plus extrapolation distance), using the initial CANDU 900 depletion snapshot and zone levels, unperturbed nuclear data, unperturbed device incremental cross-sections, and aged adjuster rods with irradiation of 182,500 EFPH, the coupled TRACE/PARCS model, running the modified codes TRACE_Mac1.1 and PARCS_Mac1.1, it reproduces the overall behaviour of the transient at the CANDU 900 station.
Figure 5 shows a comparison of PARCS_Mac1.1 channel powers against reference SORO channel powers obtained from independent simulations by station staff, given identical bundle depletions and equivalent boron moderator poison. The channel power discrepancies for the vast majority of fuel channels do not exceed 2%, except for edge channels (channels with fewer than four neighbours) and some channels in the vicinity of the liquid zone compartments. In the PARCS model, the edge channels use a different set of cross-sections calculated by homogenising an edge channel lattice model.
Table 1 summarises the initial zone levels calculated by the coupled RRS model compared to SORO predictions made by the operators. Most of the predicted zone levels in this simulation are very close to the SORO predictions, except for zones 3, 4, and 11. For zone 3, the discrepancy corresponds to a rather small reactivity difference as the zone is near the top of the core, and the change in fill from 70% to 84% corresponds to near the top of the core, where the flux is lower and thus the zone level change contributes less reactivity. The other two zones, 4 and 11, are in the centre of the core, and the lower zone level may indicate either excess absorption near the centre of the core or excess reflection from the periphery of the core or a combination of the two. Table 2, which summarises the liquid zone level results obtained throughout this work, shows that the centre zone level is particularly sensitive to nuclear data uncertainty. Though the values in Table 2 are calculated at the end of the shim operation rather than at the start, they show that the standard deviation is roughly 15%, thus the zone level discrepancy is well within two standard deviations when considering nuclear data uncertainty, and thus not significant.
The effect of the thermal-hydraulic channel grouping (40 channels) was evaluated by modelling 480 single channels using the conditions from the coupled model, then feeding the single-channel model conditions back into PARCS. The effect was found to be insignificant, with keff changing by −0.03 mk and the vast majority of channel powers changing by less than 0.2%.
Figure 6 presents the liquid zone level trends of the simulated reference case against the station data, with a rolling window used to reduce the noise for individual zone levels. The two subplots are synchronised based on a known offset between the initial snapshot used to start the simulation and the station data. The trends are overall similar between the station data and simulation, but some differences were noted.(i)The average liquid zone levels for the central zones (zones 4 and 11) are up to 15% lower in the simulation compared to the station data, for similar points in the transient. This indicates that there is slightly less reactivity in the centre of the core than the periphery in the simulation, compared to the station data, that the liquid zone controllers are compensating for.(ii)The maximum liquid zone levels after each adjuster pull, as well as after the refuelling run, are slightly less in the simulation than in the station data (by less than 5%).(iii)The withdrawal of adjuster bank A occurs approximately one day earlier in the simulation. The withdrawal of adjuster bank B occurs approximately 3 days earlier in the simulation, giving a discrepancy of two days for the interval between adjuster banks A and B. The remainder of the simulation shows little time discrepancy from the station data, besides the 3-day shift arising from the previous discrepancies.(iv)The zone level tilt between the two zones in each zone pair is considerably larger in the simulated core than in the real core, particularly for the 3/10 and 4/11 zone pairs.
(a)
(b)
Figure 7 shows the continuation of the simulation, where the power is reduced to 59% FP, showing the simulated reactor and zone powers along with liquid zone compartment levels. After roughly 2.5 hours, the average zone level reaches 10% with all adjusters withdrawn, and the reactor is shut down. The zone levels behave as expected, with the zone levels in the central zones (4/11) increasing, in particular, to counteract the flux peaking from the removed adjusters. The zone powers also follow expected trends, where fully drained zones decrease in power while the remaining zones increase in power. As well, towards the end of the simulation, the total power decreases by approximately 2% FP. This is due to the RRS response from the liquid zone control algorithm, from the following two effects:(1)With twelve liquid zones fully drained, the bulk controller’s effective gain is reduced, as a given negative valve lift will drain only two zones instead of all 14. Thus, a larger power error is required for a given reactivity insertion rate (dictated by the xenon transient being offset to maintain criticality).(2)The spatial controller is applying positive valve lift to the 12 drained zones (to attempt to bring their fill levels up to the average fill level) as well as applying positive valve lift to the two central zones (to attempt to bring the zoning power down to average zone power). The bulk controller must counteract this positive lift by applying a negative lift, which requires a negative power error.
(a)
(b)
This effect would not have been captured in a quasi-static calculation that assumes reactor power is held at the setpoint while controlling for keff and spatial flux, thus demonstrating the advantage of using a dynamic (transient) approach at key points in a long-term simulation.
The maximum power tilts immediately prior to shutdown, ignoring zones 4 and 11, is approximately 13%. The flux tilt increases over the course of the transient as the zone levels become fully drained or filled and cannot spatially control the flux, combined with the positive feedback effect from xenon.
3.2. Effect of Adjuster Rod Irradiation
Table 3 gives the effect of applying the adjuster depletion on the reactivity worth of the adjusters. The reactivity worth of each bank is calculated from the initial fuel depletion snapshot by calculating the reactivity change between the bank-in and bank-out states with all liquid zones filled to 50% and with all prior banks withdrawn. Overall, the reactivity worth of the adjusters is reduced by 10.92% when aged to 182,500 EFPH, with a greater effect on adjuster types 1, 2, and 4.
Table 4 summarises the differences in the timeline of the transient between the station data and the different simulations throughout this work. Overall, the length of the operation with zero and two banks out is slightly underpredicted, the length of the operation with three banks out is overpredicted, and the length of the operation with one bank out is greatly underpredicted. The greater the adjuster age, the shorter each interval becomes, with the three-bank-out interval being the most affected. With unirradiated adjuster rods, the entire operation is 1.22 days longer than for the reference case.
The adjuster depletion also influences the radial flux distribution in the core. Table 2 shows the magnitude of this effect, with the difference in zone levels in zones 4 and 11 between the simulations with irradiated versus unirradiated adjusters at 28.0% towards the end of the transient, with the unirradiated adjusters showing a greater underprediction against the station data when compared to the reference case.
Finally, Figure 8 shows the results of changing the effective adjuster rod age to the end of the simulated transient. Results are shown for fresh adjusters, as well as ages of 154,008 EFPH, 182,500 EFPH, and 195,000 EFPH. The irradiation of the adjuster rods has a significant effect on the outcome of the transient. With fresh adjusters, the reactor survives xenon poison out with a minimum zone level of 35%, achieved with all adjusters out of the core. With an adjuster rod irradiation of 154,008 EFPH, the reactor initially survives xenon poison out with a minimum zone level of 12%. However, the spatial xenon transient continues until reaching a 20% power tilt, which triggers a setback in the simulation. This reduces the power setpoint to 55.8% FP, causing the xenon transient to poison out the reactor. With adjuster rod irradiations of 182,500 EFPH (reference case) and 195,000 EFPH, the reactor poisons out without a setback, with the poison out occurring 10 minutes sooner for the 195,000 EFPH case.
The 154,008 EFPH case demonstrates that the reactor becomes unrecoverable once the average zone level is too low, as the lack of spatial control leads to a flux tilt setback. Thus, the cases that are shut down without a setback upon reaching the 10% average zone level would be expected to inevitably shut down, either from fully draining the liquid zones and becoming subcritical or from a flux tilt setback leading to the reactor becoming subcritical.
3.3. Effect of Liquid Zone Reactivity Worth
For this sensitivity case, a perturbation was applied to the liquid zone compartment incremental cross-sections, decreasing their values by 10%. This resulted in decreasing their reactivity worth by approximately 8.6%, from 5.76 mk to 5.26 mk. This change only has a small impact on the poison-out trajectory compared to the reference case, with the adjuster movement timings essentially unchanged, and the poison out occurring roughly 10 minutes earlier. The rest of the shim transient is slightly shorter. This effect is in line with expectations, as the reactivity gained from draining the liquid zones from 50% (initial state) to 25% (the threshold for adjuster movements) and then to 10% (the threshold for poison out) is slightly decreased. These results are quantified in Tables 2 and 4.
3.4. Effect of Axial Nodalisation
For this sensitivity case, the axial nodalisation is altered. The effect on the final phase of the transient is shown in Figure 9. When a 26-plane nodalisation is used instead of the 13-plane nodalisation, the simulation predicts that the minimum zone level reaches 10%, then poisons out when the flux tilt triggers a setback. When a 12-plane nodalisation is used, the simulation predicts that the reactor survives xenon poison out with a minimum zone level of approximately 18%. There is also a small effect on the duration of the shim operation; the effect is shown in Table 4.
Table 5 shows how the nodalisation affects the predicted reactivity worth of each set of reactivity devices. Note that the 13-plane model predicts a higher liquid zone controller worth than the 12-plane model. Conversely, the 12-plane model predicts higher adjuster worths than the 13-plane model, with the effect coming from banks D, F, G, and H, which are the centre row banks. This is consistent with expectations, where a model that distributes the device incrementals over a larger volume will predict a larger reactivity worth due to the weaker spatial self-shielding. This also explains the results of the previous figures. For Figure 9, the 12-plane model survives poison out as the adjuster reactivity worths are significantly greater. The results of the 26-plane model fall intermediately between the 12-plane and 13-plane models.
It should also be noted that the effect of the nodalisation on the liquid zone compartments is of a similar magnitude to the liquid zone sensitivity study. Therefore, the effect of the nodalisation of the liquid zone compartments, independent of the adjusters, will be similar to the results presented for the effect of the liquid zone reactivity worth.
3.5. Effect of Nuclear Data Uncertainty
The comparison of the shim transient timings is included in Table 4. The station data fell within roughly two standard deviations of the average simulation results, except for the one-bank-out interval. However, the durations of the intervals are highly correlated (>96%), while the differences in station data for different intervals have opposite signs, thus nuclear data uncertainty alone cannot explain all of the differences between simulations and station data for shim durations.
Figure 10 plots the effect of nuclear data perturbation on the extended portion of the transient, with each sample, plotted as a coloured line, with time zero set as the time at which the reactor power setpoint is set to 59.0% FP. Each of the 60 samples follows a similar trend, though the timing varies. The results of the samples are as follows:(i)44 samples (73.3%) reach the 10% average zone level threshold, triggering the simulation to shut down the reactor. These can be considered to poison out from bulk xenon only, even without the effect of spatial xenon.(ii)15 samples (25.0%) do not reach the 10% average zone level threshold at first, but a reactor setback on flux tilt is subsequently triggered. These can be considered to not poison out directly from bulk xenon, but poison out due to spatial xenon increasing the flux tilt to an unacceptable level.(iii)1 sample (3.3%) does not reach the 10% average zone level threshold and does not experience a setback.
Due to the negative valve lift from the bulk controller along with spatial control being phased out below the 10% zone level, it is possible for the zoning power in drained zones to exceed the zoning power in spatially controlled zones while the zones remain drained. This was observed for zones 5 and 12. This continues until either a flux tilt setback occurs, or a crossover point is reached where the liquid zone compartment starts to refill, providing positive feedback as spatial control is phased in, at which point the liquid zone levels rapidly change to a new equilibrium.
The sensitivity to liquid zone controller gains was evaluated, with a focus on the marginal cases (i.e., those closest to the threshold of whether a setback occurs) and it was found that realistic perturbations to the controller gains did not significantly alter the outcome of the cases evaluated in this work.
This indicates that, when only nuclear data uncertainty is considered, the simulation predicts a poison out with an approximately 73% confidence, or 0.62σ, if the setback is not considered. However, if setbacks within the simulation duration are also considered, then the confidence rises to 98%, or 2.1σ. In addition, the nuclear data uncertainty does not include uncertainty in the incremental cross-sections of the reactivity devices, as perturbed incremental cross-sections were not available in the current methodology. However, the sensitivity of the outcome to the adjuster and liquid zone incremental cross-sections was determined in the prior subsections.
For comparison, 24 cases were run with unirradiated adjusters modelled and showed that none of these cases poison out with unirradiated adjusters and that the nuclear data uncertainty has far less of an impact than the adjuster rod irradiation effect.
3.6. Effect of Burnup Distribution
Figure 11 plots the effect of simulating the transient with different initial fuel burnup distributions. Unlike the simulations using the original CANDU 900 snapshot, these cases do not include any refuelling runs. Therefore, a case was also run using the original CANDU 900 snapshot, but without any refuelling runs, to be the reference against which other burnup distributions are compared. All 14 snapshots lead to a poison out, with the reference case falling roughly in the middle of the 14 snapshots. Some reactor setbacks were observed earlier in the transient than in the above sections, as the generated burnup distributions were less optimal than the CANDU 900 snapshot.
When the transient was simulated with both nuclear data uncertainty and different initial fuel distributions, with each nuclear data sample having its own corresponding core snapshot, the distribution of shim operation trajectories was found to be similar to those observed for the CANDU 900 samples in Figure 10, except that earlier setbacks were observed as in Figure 11. All 20 of the perturbed cases led to a poison out. These results are quantified in Tables 2 and 4.
For one of the nuclear data samples with its corresponding snapshot, the simulated reactor tripped on regional overpower while withdrawing adjuster bank A. This was due to the snapshot having a single channel at abnormally high power (>7.1 MW) which increased CPPF to an abnormal level (>1.15), reducing the trip margin (note: the reactor would never be operated in such a state in reality). For this sample, the snapshot was replaced with a different snapshot from the same simulated core follow.
A minority of the perturbed cases also experienced setbacks during the adjuster movements. In this case, the simulation was set up to restore the reactor power as soon as the setback condition cleared.
4. Discussion
This study demonstrates that, by coupling TRACE_Mac1.1 and PARCS_Mac1.1 along with an RRS model, shim operation can be simulated, and sensitive parameters can be evaluated. These sensitivities are most interesting to examine for the extended portion of the simulation, where the reactor power is reduced to 59% FP. The importance of these sensitivities in this phase of the simulation, as identified from this study, are presented by listing the sensitivities from most to least significant as follows, quantified using the average zone level at 105 minutes.(1)Adjuster rod irradiation (−19.5%, reference versus unirradiated)(2)Axial nodalisation (+6.9%, 12 planes versus 13 planes, due to effect on adjusters)(3)Nuclear data uncertainty (σ = 2.6%)(4)Burnup distribution (σ = 1.0%)(5)Liquid zone compartment reactivity worth (−0.6% for 10% incremental change)
The adjuster rod irradiation was the greatest factor in defining the end of the transient, especially when comparing irradiated and unirradiated adjusters. When the transient is simulated with unirradiated adjusters, they had sufficient reactivity depth to maintain criticality and prevent the core from being poisoned out by xenon, with a reasonable margin. When the transient is simulated with irradiated adjusters, the core is poisoned out. The reactivity worth of the five banks of adjusters was reduced by 0.96 mk, or 11.6%, for the irradiated case compared to the unirradiated case. Varying the irradiation around the selected reference value had a proportionally smaller effect on the transient.
The axial nodalisation had a significant effect on the end of the transient, as the centre row of adjusters, which dominate the adjuster contribution to the end of the transient, were significantly affected by the doubling of their effective volume. The calculated worth of the five banks of adjusters was 0.38 mk, or 5.3%, higher for the 12-plane model compared to the 13-plane model—nearly half of the adjuster irradiation effect! Since the adjuster worth, particularly for the centre row, dominates the outcome of the end of the transient, the 13-plane model should be more accurate in this regard, as the adjuster effective volume is closer to the original homogenisation volume. As well, it should be considered that the model used in this study also doubles the effective volume of the adjusters in the horizontal (x) direction, making this an area for future investigation.
The nuclear data uncertainty was also found to have a significant effect. When the average zone level was taken 105 minutes after the reactor power setpoint is set to 59% FP, Table 2 showed that the effect of 28,500 EFPH of adjuster depletion was 2.8% zone level, while one standard deviation of the nuclear data uncertainty was 2.6% zone level, thus one standard deviation of nuclear data uncertainty is equivalent to approximately
Thus, one standard deviation of nuclear data uncertainty, not including device incremental uncertainty, is roughly 15% of the total adjuster depletion effect. The influence of the fuel distribution in the core is equal to approximately 40% of the influence of the nuclear data uncertainty, based on a standard deviation of 1.0% zone level, compared to 2.6% zone level for the nuclear data uncertainty.
The other significant figure of merit was the timing of events from the initial snapshot up to the end of three-bank-out operations at a 25% average zone level. The reference case underpredicted the total length of the shim operation by 2.75 days (10%), primarily due to an underprediction of the one-bank-out duration by 1.73 days (32%), as well as an underprediction of the zero-bank-out duration by 1.14 days (9%). The importance of the sensitivities, from most important to least important, is as follows:(1)Nuclear data uncertainty (σ = 1.31 days)(2)Adjuster rod irradiation (−1.21 days, reference versus unirradiated)(3)Axial nodalisation (−0.82 days, 12 planes versus 13 planes)(4)Burnup distribution (σ = approximately half of nuclear data uncertainty)(5)Liquid zone compartment reactivity worth (−0.29 days for 10% incremental change)
The discrepancy thus cannot be entirely explained by the uncertainties studied in this work, especially for the one-bank-out duration. The two-bank-out and three-bank-out durations are within two standard deviations (from nuclear data uncertainty) of the station data, but in opposite directions. The discrepancies must be explained by phenomena that were not evaluated in this work, with potential phenomena as follows:(i)Nodalisation of the core physics model, outside of the evaluated axial nodalisation meshes.(ii)Differences in depletions of different adjusters of the same type (up to a 5% variation in flux, and thus depletion, was identified to exist, but not applied to the model). This effect is expected to be small (equivalent to a 5% variation in adjuster depletion, but for individual banks).(iii)Errors affecting the rate that fuel depletion changes the core reactivity, either from approximations in the lattice physics and core physics model, or parameters in the lattice physics model (e.g., fuel density).(iv)Errors affecting the power coefficient of reactivity, void coefficient of reactivity, or xenon reactivity, such as the modelling of fuel temperature or coolant voiding. A negative change to the power coefficient would add positive reactivity when the reactor power is decreased, lengthening the shim operation (results move closer to the station data). A negative change to the void coefficient would have a similar effect as the power coefficient, but only at high power.(v)Changes to thermal-hydraulic setpoints over the course of the transient that are neither captured in the model nor provided with the station data. No known changes to setpoints were identified for this transient.(vi)The presence of moderator poisons and changes in their concentration during the shim operation. While the equivalent moderator poison was estimated periodically over the course of the real-world transient using SORO, no actual measurements were taken. The simulation in this work used the equivalent moderator poison values, which went to zero prior to the first adjuster bank being withdrawn. If these estimates were inaccurate, particularly for the initial state as well as for the time at which adjuster bank A is withdrawn, then this would create a discrepancy in the shim duration. More specifically, if moderator poison was present when bank A was withdrawn, and subsequently removed, it would extend the duration of the one-bank-out shim operation.
At the time of this writing, the most likely explanations for the discrepancies in shim duration are a combination of nuclear data uncertainty and errors in the estimated moderator poison concentrations.
The discrepancy in radial power (as indicated by central zone compartment levels) was also measured. The importance of the sensitivities, from most important to least important, is as follows:(1)Adjuster rod irradiation (+28%)(2)Nuclear data uncertainty (σ = 15%)(3)Axial nodalisation (−5.0%, 12 planes versus 13 planes)(4)Liquid zone compartment reactivity worth (+5.1% for 10% incremental)
The LZC reactivity worth is ranked lower as the effect is predominantly due to an increase in the amount of water required for the same amount of spatial control. Thus, zone levels are not a good indication of radial power balance when the LZC worth changes significantly. Since the LZC worth also changes for the axial nodalisation sensitivity case, the LZC worth effect partially offsets the adjuster worth effect. Note that the difference between the reference simulation and station data was −18.6%. Therefore, while other errors may be present, the discrepancy in the station data may be explained by the nuclear data uncertainty.
5. Conclusions
This work was able to demonstrate that a shim operation could be simulated using a coupled model of TRACE_Mac1.1, PARCS_Mac1.1, and an ECI-coupled RRS model, and model depletion and adjuster movements comparably to real-world data. The coupled model is able to simulate the effect of depletion, xenon-135, thermal-hydraulics, and reactor regulating system response on a CANDU reactor.
When the transient was extended to evaluate a possible continuation, with the power reduced to 59% FP, the reference case predicted that xenon poison out would occur 2.5 hours after the power reduction. The model showed that adjuster depletion had a major impact on the outcome of the simulation, with no poison out occurring with unaged adjusters. When including cases that poison out due to flux tilt setback, 59/60 nuclear data uncertainty cases led to poison out. The effect that the burnup distribution had on the simulated transient is less significant than nuclear data uncertainty, with only roughly 40% as great of an effect on the final stage of the transient.
Overall, out of all parameters evaluated, the modelling of the adjusters had the greatest impact on the results and figures of merit for the transient. Not only did the adjuster depletion have a significant impact on the results, but the nodalisation of the adjusters also played a significant role. This work demonstrated that nodalising the core such that an appropriate equivalent volume is achieved for the adjusters is important and that a larger equivalent volume increases the estimated adjuster's worth. While the modelling of the LZCs is similarly affected, the LZCs had far less impact on the results of the transient.
The nuclear data uncertainty was shown to be able to significantly perturb the results, both in terms of the duration of the shim operation, the radial flux distribution, and the outcome at the end of the extended simulation. Therefore, the effect of nuclear data uncertainties cannot be ruled out as a contributing factor to the observations made in this work or for other similar studies.
6. Future Work
While this work was able to generally model the real-world shim operation and identify several important sensitivities in the modelling of a shim transient, some unaccounted-for discrepancies were still identified, and some sensitivities were identified as candidates for future analysis. The results of this work present several areas for future investigation.
This work identified nodalisation of the core physics model as a significant contributing factor, due to the effect on the modelling of the adjuster rods. However, this work only tested a few nodalisations in the axial direction, and these nodalisations do not completely solve the issue of adjuster incremental cross-sections being diluted over a larger volume. In particular, in the horizontal (x-axis) direction, the incremental cross-sections are distributed between two adjacent channels.
A horizontal subdivision of the model was not pursued in this study due to the extra effort required in PARCS when compared to axial renodalisation. Currently, there is no method in PARCS to specify independent meshes for fuel assemblies/channels and for the macroscopic cross-sections, the latter of which is needed for reactivity device incremental cross-sections. This contrasts with some CANDU-specific core physics codes, such as RFSP. Therefore, PARCS would treat these subdivided channels as independent fuel assemblies/channels, requiring extra effort to map powers and fluxes to TRACE and particularly to map powers and fluxes to the ECI through signal variables as is done in the current version of PARCS_Mac1.1 and TRACE_Mac1.1. The execution time of the simulation would also be greatly increased.
However, due to the identified importance of nodalisation, it is recommended that future studies further quantify the effect of the nodalisation on the adjuster rod reactivity worth. The current model was able to quantify the effect of changing adjuster worths on the simulation of a shim transient, thus a study that estimates the change in adjuster reactivity worth due to other modelling parameters (such as nodalisation in the x-axis) may be used to estimate the effect these parameters would have on the shim transient. Future work may also investigate other methods for correcting nodalisation effects, such as extending the rod cusping correction methodology already present in PARCS or applying similar methodologies to calculate correction factors. Monte Carlo models of a full core or mini-core may be created to calculate reactivity device worths that are unaffected by nodalisation, to compare against PARCS models.
Secondly, there were some large discrepancies in the duration of different phases of the shim operation, particularly for the one-bank-out phase. One significant “unknown” that was identified was the state of moderator poison in the core. While an “equivalent moderator poison” was calculated at the time of the operation using SORO, no true moderator poison measurements were made. Any error in the moderator poison estimation, either at the initial state or upon withdrawing an adjuster bank, would impact the timing of each phase of the shim operation, as the coupled simulation assumed these estimates to be accurate. Another possible area for future investigation is the thermal-hydraulic feedback.
Several recommendations can thus be made for studies aiming to reproduce these results or further expand on this work. The first recommendation is to obtain data from other shim operations besides the one used in this work, ideally using a case where moderator poison can be definitively ruled out as a contributing factor. These may include shim operations from either 900 MW class or 600 MW class CANDU reactors. The second recommendation would be to obtain thermal-hydraulic data, as the dataset provided for this work did not include any thermal-hydraulic data against which to compare the simulated thermal-hydraulic model. Another related recommendation is to obtain real quantities for RRS parameters, such as liquid zone control gains, rather than estimating these parameters, to more accurately model RRS feedback, eliminating a source of uncertainty.
As an example, consider the hypothesis that the large underestimation of the duration of the one-bank-out phase of the operation was due to a discrepancy between the estimated moderator poison concentration (through equivalent boron used in SORO) and the actual moderator poison concentration. If similar discrepancies are not found when simulating different shim operations, it suggests that the hypothesis may still hold, and alternative hypotheses, such as a systematic modelling error, may be ruled out, or at least considered less likely. If the discrepancy does recur, then it suggests that either the hypothesis may be ruled out (and alternative hypotheses that may be common to both events considered), that a common discrepancy in estimating the moderator poison occurred in the events (the likelihood of which can be estimated by comparing the timelines of both events, particularly the use of moderator poison and the estimates of equivalent moderator poison), or that the discrepancy was co-incidental, with two different causes for the two events.
Another area to be investigated is the use of nuclear data and the estimation of the uncertainty of nuclear data. Since nuclear data uncertainty was found to have a significant impact, the use of newer nuclear datasets, such as ENDF/B-VIII, may impact both the best estimate as well as uncertainty results. Additionally, the methodology in this study was limited as it was unable to apply the stochastic nuclear data perturbations (from Sampler) to the reactivity device incremental cross-section calculations (using Serpent, which is outside of the SCALE framework supported by Sampler). An improved uncertainty analysis would integrate the reactivity device calculations into the nuclear data uncertainty analysis, such as by migrating the supercell calculation to CSAS-Shift in the upcoming SCALE 6.3 release [16]. Finally, the methodology described in reference [7] could be applied to the uncertainty analysis, performing a core follow from a much earlier state to the point of the shim event using the reactor’s operational history, to achieve core depletion profiles consistent with their respective nuclear data perturbations.
The methodology used in this work may also be applied to the analysis of design changes, such as the adoption of advanced fuel cycles, which had been previously studied for DUPIC fuel [2, 3] as well as thorium-based fuel [4], with PARCS_Mac serving as a viable alternative to RFSP and DONJON while allowing coupling to TRACE_Mac for thermal-hydraulic feedback as well as taking advantage of the uncertainty propagation capabilities present in SCALE using Sampler.
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. The detailed data used to support this work is also unavailable as it was provided confidentially, with approval to present the data included in this article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors would like to thank Michael Tucker for his prior work in developing the TRACE and PARCS CANDU models which were used in this research. In addition, the authors would like to thank him for his support in running the lattice physics calculations using SCALE/TRITON, Sampler, and Serpent 2, to generate the macroscopic cross-sections for PARCS, as well as in running simulated core follows to generate burnup snapshots, based on his previously developed methodology [7]. Both the macroscopic cross-sections and burnup snapshots were used in this work in addition to his own work. 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).