Multigroup diffusion

After her rotation in the materials and aging group, Hana Okoye returns to the core physics group at Meridian Nuclear Engineering. Dr. Teodora Vasquez, who taught her one- and two-group theory a year earlier, has a new assignment for her: join the team that designs reload cores for TVA's pressurized water reactors and supports the core design of the BWRX-300 at Clinch River. "Last year you learned why a reactor is critical," Teodora says. "This year you learn how we predict, cycle by cycle, where the power is in a real core, and how we arrange fresh and burned fuel to keep it where we want it. That starts with more energy groups than two."

Why more than two groups

In Reactor Physics I, the neutron population was described by one or two energy groups: a fast group for neutrons born in fission and slowing down, and a thermal group for neutrons in equilibrium with the moderator. Two groups capture the main physics of a thermal reactor surprisingly well. They show why a reactor needs a moderator, how leakage depends on size, why the thermal flux peaks in a reflector and how a soluble absorber changes the multiplication factor.

Two groups are not enough, however, when the answers must be accurate to a fraction of a percent, or when the spectrum changes strongly from place to place. Neutron cross sections vary by many orders of magnitude between 10 MeV and a hundredth of an electronvolt. The capture cross section of uranium-238 has huge resonances between about 6 eV and a few keV. Plutonium-239 and plutonium-240 add resonances of their own near 0.3 eV and 1 eV as fuel burns. Boron and gadolinium absorb most strongly at thermal energies, while structural materials matter mostly at higher energies. Collapsing all of this into two numbers per material hides effects that move the multiplication factor by hundreds of pcm and shift power between assemblies by several percent.

The answer is to divide the energy range into G groups, write a diffusion equation for each group and couple the groups through scattering and fission. Lattice codes use dozens to hundreds of groups for the detailed cell and assembly calculation. Core simulators then usually work with two groups, but with constants that were generated from a many-group calculation of the same assembly. Multigroup theory is the bridge between the two.

The multigroup diffusion equations

Number the groups from g = 1 at the highest energies to g = G at the lowest. For each group, a balance holds in every small volume of the core: neutrons leak out by diffusion, are removed by absorption or by scattering to other groups, and are supplied by fission and by scattering in from other groups. In steady state, with an eigenvalue k that adjusts the fission source so that a solution exists, the equation for group g is:

Here D_g is the diffusion coefficient in group g, Σ_{r,g} is the removal cross section (absorption plus scattering out of the group), χ_g is the fraction of fission neutrons born in group g, νΣ_{f,g'} is the fission neutron production cross section in group g' and Σ_{s,g'→g} is the cross section for scattering from group g' into group g. In a thermal reactor, fission neutrons are born almost entirely in the top one or two groups, so most χ_g are zero.

Two features of the scattering terms simplify the problem. First, neutrons above about 1 eV almost never gain energy in a collision, because the moderator atoms are moving much more slowly than the neutrons. Scattering in those groups goes only downward, and the matrix of scattering cross sections is lower triangular. Second, in the thermal range, neutrons can gain energy from the moving moderator atoms. Upscattering couples the thermal groups to each other in both directions, which is why thermal group structures need more care and why codes iterate over the thermal groups.

Choosing a group structure

The group boundaries are not arbitrary. A good structure places boundaries where the cross sections or the flux change character, so that the averaging inside each group is gentle. For thermal reactors, common choices include a boundary near 0.625 eV, roughly where the cadmium cutoff lies and where neutrons above it are treated as epithermal; a boundary near 5.5 keV or 9.1 keV between the resolved resonance region and the fast range; and a boundary near 0.82 MeV or 1.35 MeV near the top of the slowing-down region. Lattice codes use much finer structures, such as 69, 172 or several hundred groups, especially in the resonance region where self-shielding matters.

The flux used to weight the averages inside each group is the key to good group constants. If the group constants are defined so that each group conserves the reaction rates of the fine-energy calculation, the multigroup solution reproduces those reaction rates exactly for the case the constants were made from. For a removal cross section, for example:

The difficulty is that the weighting flux is the answer the calculation is trying to find. In practice a lattice code computes a fine-group flux for a representative piece of the reactor, such as a pin cell or a whole assembly with reflective boundaries, and uses that flux to collapse the cross sections to a few groups. The few-group constants are then valid for that composition, temperature and burnup, and are tabulated against those conditions for the core simulator. Lesson 2 describes the lattice step in detail, and lesson 11 shows how the two steps fit together in production codes.

Solving the equations: the power iteration

The multigroup equations form a generalized eigenvalue problem. Written in operator form, with M the migration and loss operator and F the fission production operator, the problem is M φ = (1/k) F φ. The fundamental mode is the solution with the largest k, and its flux is positive everywhere. The standard solution method, used in essentially every diffusion code, is the power iteration, also called the source iteration:

1. Guess a flux shape and a value of k. 2. Compute the fission source from the current flux. 3. Solve the group equations one after another, from the fastest group to the slowest, using the fission source and the scattering sources from groups already solved. 4. Compute a new fission source and update k as the ratio of the new to the old total fission source, multiplied by the old k. 5. Repeat until k and the fission source shape stop changing.

Each group equation in step 3 is an ordinary fixed-source diffusion problem, which in one dimension becomes a tridiagonal linear system that can be solved directly. In two and three dimensions, codes use iterative inner solvers, but the outer structure is the same. The rate of convergence depends on the dominance ratio, the ratio of the second largest eigenvalue to the largest. Large, loosely coupled cores such as CANDU reactors and big light-water reactors have dominance ratios close to one, so the plain power iteration converges slowly. Production codes accelerate it with methods such as Chebyshev extrapolation, Wielandt shifts or coarse-mesh rebalancing.

A multigroup slab code

The program below solves the multigroup diffusion equations for a bare slab, 300 cm thick, by finite differences and the power iteration. It runs a two-group case with the constants from the example and a four-group case with constants chosen to resemble a light-water lattice. The four-group case shows how the spectrum is shared among the groups at the center of the slab.

# Multigroup diffusion in a bare slab: finite differences and power iteration.
# Group constants are illustrative homogenized values for a light-water lattice
# (cm, 1/cm). Group 1 is fast, the last group is thermal.
def solve(G, D, sa, nsf, chi, ss, width=300.0, n=150, tol=1e-8):
    h = width / n
    def tri(g, src):
        # Solve -D phi'' + sr phi = src with zero flux at both extrapolated edges.
        sr = sa[g] + sum(ss[g][j] for j in range(G) if j != g)
        a = [-D[g] / h**2] * n; b = [2 * D[g] / h**2 + sr] * n; c = a[:]
        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] = (src[i] - (a[i] * dp[i - 1] if i else 0.0)) / m
        x = [0.0] * n
        for i in reversed(range(n)):
            x[i] = dp[i] - (cp[i] * x[i + 1] if i < n - 1 else 0.0)
        return x
    phi = [[1.0] * n for _ in range(G)]; k = 1.0
    for it in range(500):
        F = [sum(nsf[g] * phi[g][i] for g in range(G)) for i in range(n)]
        new = []
        for g in range(G):
            src = [chi[g] * F[i] / k + sum(ss[j][g] * (new[j][i] if j < g else phi[j][i])
                   for j in range(G) if j != g) for i in range(n)]
            new.append(tri(g, src))
        F2 = [sum(nsf[g] * new[g][i] for g in range(G)) for i in range(n)]
        k_new = k * sum(F2) / sum(F)
        if abs(k_new - k) < tol: k = k_new; phi = new; break
        k = k_new; phi = new
    return k, phi, it + 1

# Two-group constants
D2 = [1.40, 0.38]; sa2 = [0.0100, 0.0850]; nsf2 = [0.0070, 0.1350]
chi2 = [1.0, 0.0]; ss2 = [[0, 0.0180], [0, 0]]
k2, phi2, it2 = solve(2, D2, sa2, nsf2, chi2, ss2)
mid = len(phi2[0]) // 2
print(f"two-group k = {k2:.5f} after {it2} iterations; thermal/fast flux at center = {phi2[1][mid]/phi2[0][mid]:.3f}")

# Four-group constants: fast, epithermal, resonance, thermal (illustrative)
D4 = [2.10, 1.10, 0.80, 0.38]
sa4 = [0.0030, 0.0040, 0.0150, 0.0850]
nsf4 = [0.0045, 0.0010, 0.0060, 0.1350]
chi4 = [0.98, 0.02, 0.0, 0.0]
ss4 = [[0, 0.060, 0.0, 0.0], [0, 0, 0.045, 0.0], [0, 0, 0, 0.040], [0, 0, 0, 0]]
k4, phi4, it4 = solve(4, D4, sa4, nsf4, chi4, ss4)
spec = [phi4[g][mid] for g in range(4)]; s = sum(spec)
print(f"four-group k = {k4:.5f} after {it4} iterations")
print("center spectrum shares:", ", ".join(f"g{g+1} {v/s:.2f}" for g, v in enumerate(spec)))

two-group k = 1.26368 after 228 iterations; thermal/fast flux at center = 0.212
four-group k = 1.18710 after 179 iterations
center spectrum shares: g1 0.27, g2 0.33, g3 0.27, g4 0.13

The two-group slab result is lower than the infinite-medium value because neutrons leak from the faces of the slab. The ratio of thermal to fast flux at the center, about 0.2, is typical of a tight light-water lattice, where the thermal flux is depressed by strong absorption in the fuel and in the water itself. In a heavy-water CANDU lattice, the same ratio is several times larger, because heavy water absorbs very weakly and the lattice is widely spaced. The four-group case gives a different k because its constants were chosen independently; in a real methodology the four-group and two-group constants would both come from the same fine-group calculation, and the two-group result would be checked against the four-group one.

How many groups does a core need?

For the core calculation itself, two groups have been the industry standard for light-water reactors for decades, with the group boundary near 0.625 eV. The reason is economy: a full three-dimensional core model may have tens of thousands of nodes, the calculation is repeated for many burnup steps and many operating states, and two groups are enough when the assembly-level constants are generated carefully. Some modern codes use four, eight or more groups in the core, especially for cores with strong spectral interfaces, such as fresh fuel next to highly burned fuel, mixed fuel types or reflector boundaries, where two-group constants generated in an infinite lattice are least accurate.

For boiling water reactors such as the BWRX-300, the void fraction changes strongly along each assembly, so the spectrum hardens toward the top of the core. Two-group constants must therefore be tabulated not only against burnup but also against void fraction and void history. The PWR cores at TVA's plants see a more uniform moderator density, but boron concentration and moderator temperature still enter the tables.