Reactor Criticality & Geometric Buckling in Core Geometries

1. Derivation of the Steady-State Reactor Equation

Consider a bare, homogeneous nuclear reactor operating at steady state. Thermal neutrons are lost by leakage ($-D \nabla^2\phi$) and absorption ($\Sigma_a \phi$). In a multiplying medium, thermal fission generates $k_\infty \Sigma_a \phi$ fast neutrons, which slow down to thermal energy. In one-group diffusion theory: $$-D \nabla^2 \phi(\vec{r}) + \Sigma_a \phi(\vec{r}) = k_\infty \Sigma_a \phi(\vec{r})$$ Rearranging terms: $$-D \nabla^2 \phi(\vec{r}) = (k_\infty - 1) \Sigma_a \phi(\vec{r})$$ Dividing both sides by $D$: $$\nabla^2 \phi(\vec{r}) + \frac{(k_\infty - 1) \Sigma_a}{D} \phi(\vec{r}) = 0$$ Recalling that the thermal diffusion area is $L^2 = D / \Sigma_a$: $$\mathbf{\nabla^2 \phi(\vec{r}) + B^2 \phi(\vec{r}) = 0}$$ This is the famous Helmholtz Reactor Critical Wave Equation!

2. The Separation of Variables

The eigenvalue problem $\nabla^2\phi + B^2\phi = 0$ requires that the neutron flux $\phi(\vec{r})$:

1. Must be non-negative everywhere inside the physical core: $\phi(\vec{r}) \ge 0$.
2. Must vanish at the extrapolated boundaries: $\phi(\vec{r}_{\text{extrapolated}}) = 0$.
3. Must be finite, real, and continuous throughout the interior.

The lowest eigenvalue $B^2$ corresponding to the fundamental strictly positive eigenmode defines the geometric buckling $B_g^2$ of the reactor core!

1. Material Buckling $B_m^2$

The material buckling $B_m^2$ depends strictly on the nuclear, isotopic, and material composition of the fuel-moderator mixture: $$\mathbf{B_m^2 \equiv \frac{k_\infty - 1}{M^2} = \frac{k_\infty - 1}{L^2 + \tau}}$$ It characterizes the intrinsic neutron production capacity per unit migration area of the nuclear composition, completely independent of the core shape or size!

2. Geometric Buckling $B_g^2$

The geometric buckling $B_g^2$ depends strictly on the physical size and geometric shape of the reactor core, determined by the spatial Laplacian curvature: $$\mathbf{B_g^2 \equiv - \frac{\nabla^2 \phi}{\phi}}$$ Large cores have small geometric curvature (low leakage, small $B_g^2$). Small cores have sharp curvature (high surface-to-volume ratio, intense leakage, large $B_g^2$).

3. The Criticality Balance

Recall the Six-Factor formula $k_{\text{eff}} = \frac{k_\infty}{1 + M^2 B_g^2}$. Rearranging: $$k_{\text{eff}} = 1 \iff 1 + M^2 B_g^2 = k_\infty \iff \mathbf{B_g^2 = \frac{k_\infty - 1}{M^2} = B_m^2}$$ Hence, the fundamental reactor criticality criterion is:

• Subcritical ($k_{\text{eff}} < 1$): $B_g^2 > B_m^2$ (Core is too small; leakage dominates production).
• Critical ($k_{\text{eff}} = 1$): $\mathbf{B_g^2 = B_m^2}$ (Exact balance between geometry and nuclear material).
• Supercritical ($k_{\text{eff}} > 1$): $B_g^2 < B_m^2$ (Core is larger than critical size).

1. Infinite Slab Reactor (Thickness $a$)

Consider an infinite slab bounded by planes at $x = \pm a/2$. With extrapolated thickness $\tilde{a} = a + 2d$: $$\frac{d^2\phi}{dx^2} + B^2 \phi = 0$$ General solution: $\phi(x) = A \cos(B x) + C \sin(B x)$. By spatial symmetry across the midplane $x = 0$, $C = 0$. Boundary condition: $\phi(\pm \tilde{a}/2) = 0 \implies \cos(B \tilde{a}/2) = 0 \implies B \frac{\tilde{a}}{2} = \frac{\pi}{2}$. $$\mathbf{B_g^2 = \left( \frac{\pi}{\tilde{a}} \right)^2, \qquad \phi(x) = \phi_0 \cos\left( \frac{\pi x}{\tilde{a}} \right)}$$

2. Bare Spherical Reactor (Radius $R$)

Consider a sphere of physical radius $R$ and extrapolated radius $\tilde{R} = R + d$. In spherical coordinates: $$\frac{1}{r^2} \frac{d}{dr}\left( r^2 \frac{d\phi}{dr} \right) + B^2 \phi = 0$$ Substitute $w(r) = r \phi(r)$: $$\frac{d^2 w}{dr^2} + B^2 w = 0 \implies w(r) = A \sin(B r) + C \cos(B r)$$ Since the flux at the center must be finite, $\lim_{r \to 0} \frac{w(r)}{r} < \infty \implies C = 0$: $$\phi(r) = A \frac{\sin(B r)}{r}$$ Boundary condition: $\phi(\tilde{R}) = 0 \implies \sin(B \tilde{R}) = 0 \implies B \tilde{R} = \pi$. $$\mathbf{B_g^2 = \left( \frac{\pi}{\tilde{R}} \right)^2, \qquad \phi(r) = \phi_0 \frac{\sin(\pi r / \tilde{R})}{\pi r / \tilde{R}}}$$ where $\phi_0 = A B$ is the central peak flux.

1. Infinite Cylindrical Reactor (Radius $R$)

In cylindrical coordinates with axial and azimuthal symmetry: $$\frac{1}{r} \frac{d}{dr}\left( r \frac{d\phi}{dr} \right) + B^2 \phi = 0$$ This is Bessel's differential equation of order zero. The general solution is: $$\phi(r) = A J_0(B r) + C Y_0(B r)$$ Since the Bessel function of the second kind diverges logarithmically at the origin ($Y_0(0) \to -\infty$), physical regularity requires $C = 0$: $$\phi(r) = \phi_0 J_0(B r)$$ Boundary condition: $\phi(\tilde{R}) = 0 \implies J_0(B \tilde{R}) = 0$. The first positive zero of the zeroth-order Bessel function is $\nu_{0, 1} \approx 2.40483 \approx 2.405$: $$B \tilde{R} = 2.405 \implies \mathbf{B_g^2 = \left( \frac{2.405}{\tilde{R}} \right)^2, \qquad \phi(r) = \phi_0 J_0\left( \frac{2.405 r}{\tilde{R}} \right)}$$

2. Rectangular Parallelepiped (Cuboid $a \times b \times c$)

Separating variables $\phi(x, y, z) = X(x) Y(y) Z(z)$ with origin at the center: $$\frac{X''}{X} + \frac{Y''}{Y} + \frac{Z''}{Z} + B^2 = 0 \implies -B_x^2 - B_y^2 - B_z^2 + B^2 = 0$$ Enforcing zero flux at $\pm \tilde{a}/2, \pm \tilde{b}/2, \pm \tilde{c}/2$: $$\mathbf{B_g^2 = \left(\frac{\pi}{\tilde{a}}\right)^2 + \left(\frac{\pi}{\tilde{b}}\right)^2 + \left(\frac{\pi}{\tilde{c}}\right)^2}$$ $$\mathbf{\phi(x, y, z) = \phi_0 \cos\left( \frac{\pi x}{\tilde{a}} \right) \cos\left( \frac{\pi y}{\tilde{b}} \right) \cos\left( \frac{\pi z}{\tilde{c}} \right)}$$ For a perfect cube ($\tilde{a} = \tilde{b} = \tilde{c}$): $B_g^2 = 3 (\pi/\tilde{a})^2$.

1. Solution for the Finite Cylinder (Radius $R$, Height $H$)

Combining radial Bessel dependence with axial cosine dependence: $$\nabla^2 \phi = \frac{1}{r} \frac{\partial}{\partial r}\left( r \frac{\partial \phi}{\partial r} \right) + \frac{\partial^2 \phi}{\partial z^2}$$ Separating variables $\phi(r, z) = \mathcal{R}(r) \mathcal{Z}(z)$: $$\mathbf{B_g^2 = B_r^2 + B_z^2 = \left( \frac{2.405}{\tilde{R}} \right)^2 + \left( \frac{\pi}{\tilde{H}} \right)^2}$$ The 2D spatial flux profile is: $$\mathbf{\phi(r, z) = \phi_0 J_0\left( \frac{2.405 r}{\tilde{R}} \right) \cos\left( \frac{\pi z}{\tilde{H}} \right)}$$

2. The Optimum Cylinder Dimensions ($H/D = 1.082$)

For a given critical material buckling $B_m^2$, what cylinder radius $R$ and height $H$ minimize the total physical volume $V = \pi R^2 H$? We seek to minimize $V(\tilde{R}, \tilde{H}) = \pi \tilde{R}^2 \tilde{H}$ subject to constraint $\left( \frac{2.405}{\tilde{R}} \right)^2 + \left( \frac{\pi}{\tilde{H}} \right)^2 = B^2$. Using Lagrange multipliers or substitution $\tilde{H} = \frac{\pi}{\sqrt{B^2 - (2.405/\tilde{R})^2}}$: $$\frac{dV}{d\tilde{R}} = 0 \implies \tilde{H}^2 = 2 \left( \frac{\pi}{2.405} \right)^2 \tilde{R}^2 = 2 \left( \frac{3.14159}{2.405} \right)^2 \tilde{R}^2 \approx 3.409 \, \tilde{R}^2$$ $$\frac{\tilde{H}}{\tilde{R}} = \sqrt{3.409} \approx 1.846 \implies \mathbf{\frac{\tilde{H}}{\tilde{D}} = \frac{1.846}{2} \approx \mathbf{1.082}}$$ A finite cylinder achieves minimum critical volume when its height is approximately $8.2\%$ greater than its diameter!

3. Canonical Geometry Comparison Table

Geometry	Buckling $B_g^2$	Flux Distribution $\phi(\vec{r}) / \phi_0$	Critical Volume $V_c$
Sphere (Radius $\tilde{R}$)	$(\pi/\tilde{R})^2$	$\frac{\sin(\pi r / \tilde{R})}{\pi r / \tilde{R}}$	$\mathbf{129.88 / B^3}$ (Absolute Minimum!)
Optimum Cylinder ($\tilde{H} = 1.082 \tilde{D}$)	$(2.405/\tilde{R})^2 + (\pi/\tilde{H})^2$	$J_0(2.405 r/\tilde{R}) \cos(\pi z/\tilde{H})$	$\mathbf{148.28 / B^3}$ ($+14.2\%$ vs Sphere)
Cube (Side $\tilde{a}$)	$3(\pi/\tilde{a})^2$	$\cos(\pi x/\tilde{a})\cos(\pi y/\tilde{a})\cos(\pi z/\tilde{a})$	$\mathbf{161.45 / B^3}$ ($+24.3\%$ vs Sphere)

The sphere has the smallest surface-to-volume ratio, minimizing neutron leakage and thus requiring the smallest critical mass of fissile fuel.

1. The Flux Peaking Factor $\Omega$

Because neutron fission power is proportional to local flux $P(\vec{r}) \propto \phi(\vec{r})$, the thermal margin (critical heat flux, fuel centerline melting limit) is dictated by the peak-to-average flux ratio $\Omega$: $$\mathbf{\Omega \equiv \frac{\phi_{\max}}{\bar{\phi}} = \frac{\phi_0}{\frac{1}{V} \int_V \phi(\vec{r}) \, dV}}$$ Values for bare geometries:

• Infinite Slab: $$\bar{\phi} = \frac{1}{\tilde{a}} \int_{-\tilde{a}/2}^{\tilde{a}/2} \phi_0 \cos\left( \frac{\pi x}{\tilde{a}} \right) dx = \frac{2}{\pi} \phi_0 \implies \mathbf{\Omega_{\text{slab}} = \frac{\pi}{2} \approx 1.571}$$
• Bare Sphere: $$\bar{\phi} = \frac{1}{\frac{4}{3}\pi \tilde{R}^3} \int_0^{\tilde{R}} \phi_0 \frac{\sin(\pi r/\tilde{R})}{\pi r/\tilde{R}} 4\pi r^2 dr = \frac{3}{\pi^2} \phi_0 \implies \mathbf{\Omega_{\text{sphere}} = \frac{\pi^2}{3} \approx 3.290}$$
• Finite Cylinder: $$\Omega_{\text{cyl}} = \Omega_{\text{radial}} \times \Omega_{\text{axial}} = \left[ \frac{2.405}{2 J_1(2.405)} \right] \times \left[ \frac{\pi}{2} \right] \approx 2.316 \times 1.571 \approx \mathbf{3.638}$$

A bare finite cylinder has a peak power density $3.64$ times greater than its core average! Real power reactors flatten this distribution using reflectors, fuel zoning, and burnable poison loading.

1. Physics of the Reflector

Surrounding a bare multiplying core with an unmultiplied moderating shell (the reflector, e.g. $\text{H}_2\text{O}, \text{D}_2\text{O}, \text{Be}, \text{C}$) fundamentally improves reactor economics:

1. Neutron Economy: Fast and thermal neutrons leaking from the core scatter in the reflector and return to the core (neutron albedo effect).
2. Critical Mass Reduction: Returning neutrons reduce the required critical core dimensions from $R_{\text{bare}}$ to $R_{\text{core}}$.
3. Flux Flattening: Reflected neutrons thermalize in the reflector, producing a characteristic thermal flux peak near the core-reflector interface, greatly reducing peak-to-average power peaking!

2. Two-Region Diffusion Formulation

For a reflected spherical core with core radius $R$ and reflector outer radius $R_R$:

• Core Region ($0 \le r \le R$): Multiplying medium: $$\nabla^2 \phi_c + B^2 \phi_c = 0 \implies \phi_c(r) = A \frac{\sin(B r)}{r}$$
• Reflector Region ($R \le r \le \tilde{R}_R$): Non-multiplying medium: $$\nabla^2 \phi_r - \frac{1}{L_r^2} \phi_r = 0 \implies \phi_r(r) = C \frac{\sinh\left( \frac{\tilde{R}_R - r}{L_r} \right)}{r}$$

Enforcing flux and current continuity at the core-reflector interface $r = R$: $$\phi_c(R) = \phi_r(R), \qquad D_c \left.\frac{d\phi_c}{dr}\right|_R = D_r \left.\frac{d\phi_r}{dr}\right|_R$$ Dividing the current boundary condition by the flux condition yields the logarithmic derivative matching condition: $$D_c \left[ B \cot(B R) - \frac{1}{R} \right] = D_r \left[ -\frac{1}{L_r} \coth\left(\frac{\tilde{R}_R - R}{L_r}\right) - \frac{1}{R} \right]$$

3. Reflector Savings $\delta$

The reflector savings $\delta$ is the reduction in core dimension enabled by the reflector: $$\mathbf{\delta \equiv R_{\text{bare}} - R_{\text{reflected}}}$$ For an infinite reflector ($R_R \to \infty$) with matched diffusion properties ($D_c \approx D_r$): $$\mathbf{\delta \approx L_r}$$ The reflector savings is approximately equal to the thermal diffusion length $L_r$ of the reflector material! For graphite ($L_r \approx 50\text{ cm}$), reflector savings can exceed $30\text{ to }40\text{ cm}$, shrinking critical core volume by over $50\%$!