One-group diffusion in bare reactors

With the diffusion equation in hand, Teodora sets Hana the classic exercise: take a uniform core with no reflector, a bare reactor, and find the size at which it becomes critical and the shape the flux takes. The answer is idealized, she warns, because real cores are reflected and far from uniform. But it introduces the two ideas that run through all of core design, buckling and the balance between material and geometry.

The one-group reactor equation

In a multiplying medium the source in the diffusion equation is fission. Writing the fission source divided by \(k\), so that the equation can always be balanced, gives the one-group criticality equation:

Rearranged, it has the form of a Helmholtz equation:

For a given shape, the equation with zero-flux conditions at the extrapolated boundary has a non-negative solution everywhere only for one value of \(B^2\), the smallest eigenvalue of the Laplacian for that shape. That value depends only on the geometry, so it is called the geometric buckling \(B_g^2\). The name comes from the mathematics of buckled plates: \(B^2\) measures how sharply the flux shape curves, and therefore how much leakage the shape implies.

Material buckling and criticality

The composition defines a second quantity. Setting \(k = 1\) in the definition of \(B^2\) gives the material buckling,

The reactor is critical when the two bucklings are equal: the curvature the material can support matches the curvature the geometry demands. If \(B_m^2 > B_g^2\), the core is larger than critical and supercritical; if \(B_m^2 < B_g^2\), it is subcritical. In general,

which identifies the one-group non-leakage probability as \(1/(1 + L^2 B_g^2)\). With fast neutron migration included through the Fermi age \(\tau\), the corresponding modified one-group result uses the migration area \(M^2 = L^2 + \tau\) in place of \(L^2\). This formula is the most useful single result in elementary reactor physics, because it splits \(k_{\text{eff}}\) into a material part and a geometric part.

Geometric buckling for simple shapes

Solving the Helmholtz equation in standard shapes gives the fundamental flux modes and their bucklings. All dimensions are extrapolated dimensions, including the distance \(d\) beyond the physical surface.

For a fixed volume, the shape with the smallest buckling leaks least. A sphere is optimal, but power reactors are cylinders because fuel comes in long vertical rods or horizontal channels. For a cylinder of fixed volume, minimum buckling occurs at a height of about 1.08 times the diameter. Power reactors depart from that optimum for engineering reasons: light water cores are about as tall as they are wide, while CANDU cores are short, wide cylinders lying on their side, about 6 m long and 7 m in diameter.

Peak-to-average ratios

The flux shape gives a peaking factor at once. In a bare finite cylinder the maximum flux at the center is about 3.64 times the volume average: 2.32 from the radial Bessel shape and 1.57 from the axial cosine. A uniform core at that peaking would be limited by its hottest fuel at the center while the edges run cold. That is the motivation for reflectors, for differential fuel loading and for the flattening devices described in lesson 3, which together bring the overall peaking factor of a real core down to between about 1.5 and 2.5.

# One-group bare cores: geometric buckling and critical size for a given material buckling.
import math
D, Sa, nuSf = 0.95, 0.0165, 0.0170    # homogenized core constants, cm and 1/cm (illustrative)
L2 = D / Sa
k_inf = nuSf / Sa
Bm2 = (k_inf - 1) / L2               # material buckling
print(f"k_inf = {k_inf:.4f}, L^2 = {L2:.1f} cm2, material buckling = {Bm2:.3e} cm^-2")
d = 2.13 * D                         # extrapolation distance, about 0.71 transport mfp
# critical sizes (extrapolated dimensions, then subtract d)
a_slab = math.pi / math.sqrt(Bm2)
R_sph = math.pi / math.sqrt(Bm2)
# cylinder with height equal to diameter: B2 = (2.405/R)^2 + (pi/H)^2, H = 2R
R_cyl = math.sqrt(2.405 ** 2 + (math.pi / 2) ** 2) / math.sqrt(Bm2)
print(f"slab thickness: {a_slab - 2 * d:6.1f} cm")
print(f"sphere radius:  {R_sph - d:6.1f} cm")
print(f"cylinder (H = D): radius {R_cyl - d:6.1f} cm, height {2 * R_cyl - 2 * d:6.1f} cm")
# k_eff of a cylinder 3.4 m across and 3.7 m tall made of this material
R, H = 170 + d, 370 + 2 * d
Bg2 = (2.405 / R) ** 2 + (math.pi / H) ** 2
keff = k_inf / (1 + L2 * Bg2)
print(f"\n3.4 m x 3.7 m cylinder: geometric buckling {Bg2:.3e}, k_eff = {keff:.4f}")
print(f"non-leakage probability = {1 / (1 + L2 * Bg2):.4f}")

k_inf = 1.0303, L^2 = 57.6 cm2, material buckling = 5.263e-04 cm^-2
slab thickness:  132.9 cm
sphere radius:   134.9 cm
cylinder (H = D): radius  123.2 cm, height  246.4 cm

3.4 m x 3.7 m cylinder: geometric buckling 2.660e-04, k_eff = 1.0148
non-leakage probability = 0.9849

The program finishes by taking a cylinder 3.4 m in diameter and 3.7 m tall, roughly the size of a large pressurized water reactor core, and computing \(k_{\text{eff}} = 1.015\) for this material, with a non-leakage probability above 0.98. In large power cores leakage is a small correction. In a small research core, such as the SLOWPOKE-2 or TRIGA cores used for teaching, leakage is large and a reflector is essential, as lesson 6 shows.

Material and geometry together

The formula \(k_{\text{eff}} = k_\infty/(1 + L^2 B_g^2)\) can be read in two directions. A fuel designer reads it from left to right: given the lattice, how big must the core be? A core designer reads it from right to left: given the size the plant layout allows, how much multiplication must the fuel supply? Both readings show that leakage matters most when the core dimensions are not many times larger than the migration length.

Hana tries the second reading on two cores. For a large light water core, \(M^2\) is about 60 cm² and \(B_g^2\) about \(2.7 \times 10^{-4}\ \text{cm}^{-2}\), so the leakage term \(M^2 B_g^2\) is about 0.016: less than 2 percent of neutrons leak. For a CANDU core, \(M^2\) is several hundred square centimeters, but the core is about twice as wide and the buckling correspondingly smaller, so leakage is a few percent. For a small research core of the SLOWPOKE-2 or TRIGA type, with dimensions of tens of centimeters, leakage would be tens of percent without a reflector. The same equation describes all three, and the difference is the ratio of migration length to size.

Reading a core map

Teodora gives Hana a radial power map from a core-follow calculation, normalized so that the average assembly power is 1.00. The center reads about 1.3, the edge assemblies about 0.5, and the hottest assembly 1.42. Hana compares these with the bare-core numbers she has just derived. A uniform bare cylinder would have a radial peak of 2.32 relative to the average, so the real core is far flatter, because fuel management has placed fresh fuel and burnable absorbers to even it out and because the reflector returns neutrons to the edge.

The comparison is useful in reverse as well. If a calculated map ever came out close to the bare-core shape, it would suggest an error in the fuel loading input or a missing reflector model. Reactor physicists keep the simple analytic results in mind precisely so they can recognize when a large calculation has gone wrong. The bucklings and peaking factors of this lesson are the yardsticks.

The axial direction tells a similar story. A bare cosine axial shape has a peak of 1.57, but real cores show a flatter shape with features: a dip at spacer grids in a light water assembly, a shift toward the bottom early in a cycle when the coolant entering at the bottom is coolest and densest, and a shift toward the top later as fuel at the bottom burns faster. Each departure from the cosine has a physical cause that diffusion theory explains once the materials are no longer uniform.

Higher modes and why they decay

The Helmholtz equation in a slab has a whole family of solutions, \(\cos(n\pi x/a)\) for odd \(n\) and sines for even \(n\), with bucklings \((n\pi/a)^2\). Only the fundamental mode is positive everywhere. If the core is perturbed, for example by moving a control rod, the flux temporarily contains higher modes, but because their bucklings are larger they leak more and have lower effective multiplication. After a few generations they die away and the fundamental shape returns. This is why a reactor has a stable flux shape, and why computer codes find it by repeated iteration: each pass through the fission source suppresses the higher modes relative to the fundamental.

Power iteration

The criticality problem is solved numerically by the power method. Guess a flux shape, compute the fission source \(\nu\Sigma_f\phi\), solve the diffusion equation with that source divided by the current estimate of \(k\), and update

Each iteration is one solve of the kind built in lesson 4. The ratio of the second eigenvalue to the first, the dominance ratio, controls how fast the iteration converges. Large, loosely coupled cores such as CANDU have dominance ratios close to one and need many iterations or acceleration methods.

That coupling also matters physically. In a large light water core, xenon can set up slow axial power oscillations, a higher-mode effect, which operators at TVA's pressurized water reactor fleet control with rod adjustments while watching the axial offset of the power shape.