CANDU-6 REACTOR SIMULATOR  · 2064 MW_th / 728 MW_e · pressure-tube heavy-water reactor

1 · Neutron kinetics

Point kinetics with six delayed-neutron precursor groups plus three lumped photoneutron groups (D₂O photodisintegration by fission-product γ — a defining feature of heavy-water reactors, it dominates the source term after a trip and limits how negative the reactor period can become):

dn/dt = (ρ−βtot)/Λ·n + ΣλiCi + ΣλpjCpj + S₀
dCi/dt = βin/Λ − λiCi

Integrated with a fully implicit (backward-Euler) scheme so the solution stays stable at high time acceleration. Λ = 0.90 ms — roughly 40× a PWR, because thermalisation in D₂O takes many more scattering collisions. This long generation time is why CANDU tolerates slow-acting liquid-zone control.

Reactivity is reported in milli-k (1 mk = 10⁻³ Δk/k), Canadian convention. Prompt criticality is at ρ = β_eff.

2 · Reactivity feedback

The positive void coefficient is intrinsic to a lattice of natural-uranium fuel in separate pressure tubes: voiding the coolant hardens the in-channel spectrum and removes D₂O-coolant absorption, and the resulting increase in fast fission and resonance escape outweighs the loss of moderation (the bulk moderator is a separate, unvoided inventory). It is the reason CANDU carries two fast, independent, diverse shutdown systems. Drive the "imposed void (demo)" slider up and watch power run away against Doppler alone.

The MTC is the coefficient that swings hardest with burn-up: in a fresh natural-uranium core it is comfortably positive (+0.077 mk/°C), but as plutonium builds in and the lattice moves toward its equilibrium composition it collapses to a few thousandths of a mk per °C. The equilibrium fuel state therefore uses +0.005 mk/°C — a 40 K moderator excursion is worth only ~0.2 mk, which is why moderator temperature is an operational concern for the calandria and the reactivity-device calibration rather than for reactor control.

3 · Xenon-135

dI/dt = γIΣfφ − λII
dX/dt = γXΣfφ + λII − λXX − σaXXφ

γ_I=0.0639, γ_X=0.00237, T½(I)=6.57 h, T½(Xe)=9.14 h, σ_a^X=2.65×10⁶ b, φ̄=9×10¹³ n·cm⁻²·s⁻¹. Because σφ ≫ λ_X, burn-out dominates at power, so after a trip xenon grows from the iodine inventory and peaks ≈9 h later — several times the equilibrium load. The 21 adjuster rods hold ≈15 mk in reserve, which buys roughly 30 minutes of poison-override before the reactor is guaranteed to stay shut down for ~40 h. Run the simulation at ×600 after a trip to see the poison-out curve.

4 · Thermalhydraulics

Two-node fuel/coolant lump with 4.5 % of fission energy deposited directly in the moderator:

Mfcf·dTf/dt = 0.955·Pth − UAf(Tf−Tc)
Mccc·dTc/dt = UAf(Tf−Tc) − 2Ẇcp(Tc−Tin)

The steam generator is treated as a constant-temperature (boiling) sink, which has the exact solution Tin = Tsat + (Tout−Tsat)e−NTU, NTU = UA_SG/Ẇc_p. UA_f and UA_SG scale with flow^0.6 and flow^0.8. Decay heat uses three precursor groups summing to 6.7 % of full power. Saturation line: Tsat = 179.9·p0.2381 (°C, MPa) — within 1 K over 1–10 MPa.

Note the reactor outlet at ~311 °C sits essentially on the 9.9 MPa saturation line — real CANDU channels run with a few per cent exit quality, so the model generates boiling void of its own. Reduce PHT flow, raise power or drop PHT pressure and the void — and hence reactivity — climbs.

4b · Void onset

Void is not switched on at a hard T_out = T_sat threshold. Subcooled nucleate boiling starts at the sheath while the bulk coolant is still a few K subcooled, and the void fraction saturates once the channel is well into quality, so the model uses a C¹ (smooth-derivative) curve in the bulk subcooling Δ_sub = T_sat − T_out:

drive = k·ln(1 + e−Δsub/k)
αvoid = αmax·(1 − e−drive/s)

with k = 2.5 K (the subcooled-boiling band), s = 45 K and α_max = 0.45. The softplus reproduces max(0, −Δsub) without the kink, so dα/dT rises continuously from ≈0 when the channel is 10 K subcooled to ≈0.006 K⁻¹ at saturation — worth about 0.07 mk/K of positive feedback at the operating point. The "ROH subcooling" meter shows Δ_sub directly: watch it go negative as you throttle the pumps, and watch the void reactivity bar follow it smoothly rather than snapping on.

5 · Control & protection

RRS computes a demanded net reactivity ρdem = Kpe − Kdṅ (clamped to ±0.4 mk), subtracts every other reactivity component, and solves for the required zone-controller fill, rate-limited to 1.2 %/s. If the zones saturate (<15 % or >85 %) the adjuster bank is driven, exactly as the real regulating program does. SDS-1: 28 spring-assisted shutoff rods, −80 mk, 0.3 s delay then 1.4 s insertion, S-shaped worth curve. SDS-2: gadolinium nitrate injection, −300 mk in ~2 s. Trips modelled: neutron overpower 120 %FP, log-rate +10 %/s, high PHT pressure 10.5 MPa, low PHT flow 70 %, high moderator temperature 90 °C, manual. The turbine governor holds SG pressure; on a turbine trip the condenser steam discharge valves take up to 70 % steam flow so the reactor can stay at power.

6 · Preset lab scenarios

Each preset resets the plant and then plays timed steps against the simulation clock, driving the same sliders and buttons an operator would — so you can take over mid-script at any time (or hit stop script). Narration lands in the event log; the time acceleration is stepped automatically so the interesting parts are not spent waiting.

• Xenon override & poison-out — trip from equilibrium full power. The override milestones are not hard-coded: the script watches the live reactivity balance and fires when the best case with the adjusters fully out and the zones drained goes negative (and again when it recovers). In this model that lands just under an hour after the trip, with the peak near −115 mk at ≈10 h and a restart margin returning at ≈32 h.
• Positive void reactivity & SDS trip — RRS to manual, then void ramped in at ≈0.09 mk/s until SDS-1 answers. Worth noting which trip fires: a ramp this slow reaches the 120 %FP overpower setpoint (at ≈14 s, with the log rate still under 4 %/s) before it reaches the +10 %/s rate setpoint — a step insertion of the same worth would trip on rate first. The scenario ends with the coolant still voided and the core deeply subcritical, which is the shutdown margin the two SDSs are sized for.
• Power manoeuvre 100 → 70 %FP — the liquid zones take the step, then have to drain again as transient xenon builds; if they saturate the RRS backs the adjusters out.

7 · Deliberate simplifications

• Point kinetics — no spatial flux tilts, no xenon oscillations (a real CANDU core is ~20 neutron migration lengths across and is spatially unstable; that is why there are 14 independently controlled zones).
• One lumped channel; the drawn channel-power distribution is a cosine shape for illustration only.
• No ¹⁴⁹Sm (≈−5 mk at equilibrium), no ²³⁹Np/Pu transients, no burn-up distribution — burn-up is a single lumped reactivity trim.
• SG level, pressuriser level and feedwater dynamics are not modelled.
• Two-phase flow is a void-fraction correlation, not a thermalhydraulic network; nothing here is a substitute for RFSP/CATHENA.