Lattice physics

Hana's first task is to rerun the lattice calculation for the 37-element natural uranium bundle used in the refurbished Darlington units. "The core simulator never sees a fuel pin," Teodora tells her. "It sees homogenized bundles, each with a handful of cross sections. Everything about pins, sheaths, pressure tubes and the moderator between them has to be folded into those numbers by the lattice code. If the lattice step is wrong, no amount of care in the core model can fix it."

The lattice and the two-step method

A power reactor core is built from repeated units. A light-water reactor core is an array of square fuel assemblies, each holding a lattice of fuel pins, guide tubes and sometimes absorber pins, all bathed in water. A CANDU core is an array of horizontal channels on a square pitch in a tank of heavy-water moderator; each channel holds a string of fuel bundles in a pressure tube, separated from the moderator by a gas-filled calandria tube. In both cases, the core is far too detailed to model pin by pin in three dimensions over a whole cycle.

Reactor physics therefore works in two steps. In the first, the lattice step, a transport calculation is made for one assembly or bundle cell, usually in two dimensions with reflective boundaries, as if it were part of an infinite array of identical units. This calculation uses dozens or hundreds of energy groups and represents every pin, sheath and coolant region explicitly. Its output is a set of few-group, homogenized cross sections for the whole unit, together with factors that let the detailed pin powers be reconstructed later. In the second step, the core step, a diffusion or simplified transport code solves for the whole core using those homogenized constants node by node.

The pin cell and the Wigner-Seitz approximation

The simplest lattice unit is a single fuel pin and its share of moderator. In a square lattice, each pin sits at the centre of a square cell whose side is the pitch. The traditional simplification, due to Wigner and Seitz, replaces the square cell by a cylinder of the same area, with a reflective or white boundary condition on its outer surface. The fuel occupies the inner region, followed by a gap, the cladding and the moderator. Because the cell is symmetric, the flux depends only on the radius.

Within the cell, the thermal flux is not uniform. Thermal neutrons are created by slowing down in the moderator, diffuse toward the fuel and are absorbed there. The flux is therefore highest in the moderator and lowest at the centre of the fuel. Two quantities describe this flux shape:

The thermal utilization enters the four-factor formula from Reactor Physics I directly. Written with flux-weighted averages, it is:

where ζ = φ̄_M/φ̄_F is the disadvantage factor. A larger disadvantage factor means that the moderator sees relatively more flux and therefore absorbs a larger share of neutrons. Lattice design always balances two effects: adding moderator improves slowing down and reduces resonance absorption, but also increases parasitic absorption in the moderator. The best ratio of moderator to fuel depends strongly on the moderator. Light water absorbs enough that light-water lattices are tight, with a pitch only about 1.3 times the pin diameter. Heavy water absorbs so little that CANDU lattices can be widely spaced, with a pitch of about 28.6 cm between channels.

Homogenization: flux-volume weighting

To give the core simulator a single cross section for a cell or assembly, the lattice code averages each cross section over space, weighted by flux and volume, so that the homogenized cell absorbs and produces neutrons at the same total rates as the detailed one:

A plain volume average would give the fuel too much weight, because the fuel actually sees less flux than the moderator. The program below solves a one-speed thermal diffusion problem across a Wigner-Seitz cell with a uniform slowing-down source in the moderator, then computes the disadvantage factor, the thermal utilization and the flux-weighted absorption cross section. It runs one light-water-like pin cell and one heavy-water-like cell in which the fuel cluster has been smeared into a single cylinder.

# Lattice cell: one-speed thermal diffusion in a Wigner-Seitz cylinder.
# Neutrons are slowed down (uniform source) in the moderator and absorbed mostly
# in the fuel. Finite differences in r give the flux shape; flux-volume weighting
# gives homogenized constants. Values are illustrative.
import math
def cell(r_fuel, r_cell, D_f, sa_f, D_m, sa_m, n=200):
    h = r_cell / n; r = [(i + 0.5) * h for i in range(n)]
    D = [D_f if x < r_fuel else D_m for x in r]
    sa = [sa_f if x < r_fuel else sa_m for x in r]
    q = [0.0 if x < r_fuel else 1.0 for x in r]          # slowing-down source
    a, b, c = [0.0] * n, [0.0] * n, [0.0] * n
    for i in range(n):
        rw, re_ = i * h, (i + 1) * h
        Dw = 2 * D[i] * D[i - 1] / (D[i] + D[i - 1]) if i else 0.0
        De = 2 * D[i] * D[i + 1] / (D[i] + D[i + 1]) if i < n - 1 else 0.0   # reflective outer edge
        vol = r[i] * h
        a[i] = -Dw * rw / h; c[i] = -De * re_ / h
        b[i] = -(a[i] + c[i]) + sa[i] * vol; q[i] *= vol
    cp, dp = [0.0] * n, [0.0] * n
    for i in range(n):
        m = b[i] - (a[i] * cp[i - 1] if i else 0.0)
        cp[i] = c[i] / m; dp[i] = (q[i] - (a[i] * dp[i - 1] if i else 0.0)) / m
    phi = [0.0] * n
    for i in reversed(range(n)):
        phi[i] = dp[i] - (cp[i] * phi[i + 1] if i < n - 1 else 0.0)
    V = [r[i] * h for i in range(n)]
    fF = sum(phi[i] * V[i] for i in range(n) if r[i] < r_fuel) / sum(V[i] for i in range(n) if r[i] < r_fuel)
    fM = sum(phi[i] * V[i] for i in range(n) if r[i] >= r_fuel) / sum(V[i] for i in range(n) if r[i] >= r_fuel)
    vf = r_fuel**2 / r_cell**2
    f = sa_f * fF * vf / (sa_f * fF * vf + sa_m * fM * (1 - vf))
    sa_hom = (sa_f * fF * vf + sa_m * fM * (1 - vf)) / (fF * vf + fM * (1 - vf))
    return fM / fF, f, sa_hom
# light-water-like pin cell (cm): small pitch, strong moderator absorption
dis, f, sah = cell(0.41, 0.71, 0.36, 0.30, 0.16, 0.020)
print(f"light-water pin cell: disadvantage factor {dis:.3f}, thermal utilization f = {f:.3f}, homogenized sigma_a = {sah:.4f} /cm")
# heavy-water-like cell: fuel cluster smeared into a 5 cm radius, large lattice
dis2, f2, sah2 = cell(5.0, 16.0, 0.70, 0.060, 0.85, 0.00008)
print(f"heavy-water cluster cell: disadvantage factor {dis2:.3f}, thermal utilization f = {f2:.3f}, homogenized sigma_a = {sah2:.5f} /cm")
# what a simple volume-weighted average (no flux weighting) would give
vf = 0.41**2 / 0.71**2
print(f"volume-only average for the pin cell = {0.30*vf + 0.020*(1-vf):.4f} /cm (overestimates fuel absorption)")

light-water pin cell: disadvantage factor 1.054, thermal utilization f = 0.877, homogenized sigma_a = 0.1101 /cm
heavy-water cluster cell: disadvantage factor 1.790, thermal utilization f = 0.978, homogenized sigma_a = 0.00350 /cm
volume-only average for the pin cell = 0.1134 /cm (overestimates fuel absorption)

The light-water cell has a small disadvantage factor because its moderator region is thin. The heavy-water cell has a much larger one, because neutrons slow down far from the fuel and the flux builds up in the large moderator region. Even so, the heavy-water cell has the higher thermal utilization, because heavy water absorbs so little. The volume-only average for the light-water cell overestimates the absorption cross section by about 3 percent, which would bias k by several hundred pcm in a core calculation.

From pin cell to assembly

Real lattice codes go well beyond the Wigner-Seitz pin cell. A modern light-water lattice calculation models the full assembly in two dimensions: each fuel pin with its own composition, the guide tubes and instrument tube filled with water, any burnable absorber pins, the grid spacers smeared into the coolant and the water gap between assemblies. The transport equation is solved by the method of characteristics, collision probabilities or discrete ordinates in many groups. The resulting flux gives homogenized constants for the whole assembly, and also the pin-by-pin power within it.

For a CANDU bundle, the lattice cell contains the 37 fuel elements of the bundle, the heavy-water coolant inside the pressure tube, the pressure tube, the gas gap, the calandria tube and the surrounding square of heavy-water moderator. Because the coolant and moderator are separate, a lattice code can study what happens if the coolant density changes while the moderator does not, which is the origin of the coolant void reactivity effect studied in Reactor Physics I. The Canadian industry standard lattice code for CANDU has been WIMS-AECL, with DRAGON also widely used, and supercell calculations are added to model devices such as adjusters and shutoff rods between channels.

Branches, histories and the cross-section table

A lattice calculation is not done once. The fuel burns, and its composition changes; the operating state changes, and the spectrum changes with it. A lattice code therefore performs a depletion calculation through the life of the fuel at nominal conditions, called the history, and at selected burnup points it branches off to calculate constants at other conditions: higher or lower fuel temperature, moderator temperature or density, boron, void or control insertion. The results fill a multidimensional table, and the core simulator interpolates in that table for each node at each step.

Effects that depend on history rather than the current state are especially important. Fuel burned with the moderator at high void fraction in a boiling water reactor, for example, builds up more plutonium than fuel burned at low void, because the harder spectrum favours capture in uranium-238. Two nodes with the same burnup but different void histories therefore have different reactivities. Codes handle this by tabulating constants against both the instantaneous and the historical conditions, or by tracking important nuclides directly in the core calculation.

Pin power reconstruction and discontinuity factors

Homogenization throws away detail that matters for safety: the power of the hottest pin. Core simulators recover it by pin power reconstruction. The lattice code stores, for each assembly type, the ratio of each pin's power to the assembly average in an infinite lattice. The core simulator computes a smooth intra-assembly flux shape from its own nodal solution and multiplies it by these form factors to estimate each pin's power in the actual core environment.

A second correction keeps the homogenized model consistent. When a detailed heterogeneous assembly is replaced by a homogeneous block, the flux at the block surface no longer matches the true surface flux. Assembly discontinuity factors, the ratio of the heterogeneous surface flux to the homogeneous one, are stored by the lattice code and used by the core simulator to allow a controlled discontinuity in the homogeneous flux at assembly faces. With these factors, a coarse nodal model can reproduce the reaction rates of the detailed calculation much more closely.