Chapter 8

Nuclear Fission

Fission Physics

Nuclear fission is the splitting of a heavy nucleus into two (or rarely three) lighter fragments, releasing approximately 200 MeV of energy per event. The process was discovered by Hahn and Strassmann in 1938, with the theoretical explanation provided by Meitner and Frisch using the liquid drop model.

$$^{235}_{92}\text{U} + n \to ^{236}_{92}\text{U}^* \to \text{fragments} + 2\text{-}3\,n + \sim200\text{ MeV}$$

Energy Budget per Fission

Fission fragment kinetic energy:

~170 MeV

Prompt neutron kinetic energy:

~5 MeV

Prompt gamma rays:

~7 MeV

Beta particles from products:

~8 MeV

Gamma rays from products:

~7 MeV

Neutrinos (lost):

~12 MeV

Total:

~200 MeV

Fission Barrier

The liquid drop model predicts a fission barrier arising from the competition between surface energy (resists deformation) and Coulomb energy (favors deformation). The barrier height depends on the fissility parameter:

$$x = \frac{E_C^{(0)}}{2E_S^{(0)}} = \frac{Z^2/A}{(Z^2/A)_{\text{crit}}} \approx \frac{Z^2/A}{49}$$

For $^{236}$U* (compound nucleus): $Z^2/A \approx 36$, $x \approx 0.73$, barrier height ~6 MeV. The neutron separation energy of $^{236}$U is 6.5 MeV, so thermal neutrons can induce fission of $^{235}$U.

For $^{239}$U* (from $^{238}$U + n): the neutron separation energy is only 4.8 MeV, less than the barrier height (~6 MeV). Hence $^{238}$U requires fast neutrons ($E_n > 1$ MeV) for fission.

Chain Reactions and Criticality

The neutron multiplication factor $k$ determines whether a chain reaction is self-sustaining:

$$k = \frac{\text{neutrons in generation } n+1}{\text{neutrons in generation } n}$$

- $k < 1$: Subcritical -- chain reaction dies out
- $k = 1$: Critical -- steady-state chain reaction (reactor operation)
- $k > 1$: Supercritical -- exponentially growing chain reaction

The four-factor formula for an infinite medium is:

$$k_\infty = \eta\,f\,p\,\epsilon$$

where $\eta$ = neutrons per absorption in fuel, $f$ = thermal utilization factor,$p$ = resonance escape probability, $\epsilon$ = fast fission factor.

Python Simulation: Fission Product Yields

The characteristic double-humped fission product mass distribution and neutron multiplicity as a function of incident neutron energy.

Fission Product Yield Distribution

Python

Double-humped mass distribution for U-235 and Pu-239 thermal fission

script.py97 lines

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt

# ──────────────────────────────────────────────
# Fission Product Yield Distribution
# Double-humped distribution for U-235 thermal fission
# ──────────────────────────────────────────────

def fission_yield(A, A_parent=236):
    """Model fission product mass yield (percent)
    Double Gaussian approximation for thermal U-235 fission"""
    A_half = A_parent / 2.0
    # Light fragment peak
    A_L = 95.0
    sigma_L = 6.0
    # Heavy fragment peak
    A_H = A_parent - A_L
    sigma_H = 6.0

yield_L = 6.0 * np.exp(-(A - A_L)**2 / (2 * sigma_L**2))
    yield_H = 6.0 * np.exp(-(A - A_H)**2 / (2 * sigma_H**2))
    # Symmetric fission (small contribution)
    yield_sym = 0.01 * np.exp(-(A - A_half)**2 / (2 * 10**2))

return yield_L + yield_H + yield_sym

A_range = np.arange(70, 170)
yields_U235 = np.array([fission_yield(A, 236) for A in A_range])

# Pu-239 fission (slightly different peaks)
def fission_yield_Pu(A, A_parent=240):
    A_L = 100.0
    sigma_L = 6.5
    A_H = A_parent - A_L
    sigma_H = 6.5
    yield_L = 6.5 * np.exp(-(A - A_L)**2 / (2 * sigma_L**2))
    yield_H = 6.5 * np.exp(-(A - A_H)**2 / (2 * sigma_H**2))
    yield_sym = 0.015 * np.exp(-(A - A_parent/2)**2 / (2 * 10**2))
    return yield_L + yield_H + yield_sym

yields_Pu239 = np.array([fission_yield_Pu(A, 240) for A in A_range])

fig, axes = plt.subplots(1, 2, figsize=(14, 6))
fig.patch.set_facecolor('#0a0a1a')

for ax in axes:
    ax.set_facecolor('#0a0a1a')
    for spine in ax.spines.values():
        spine.set_edgecolor('#334155')

# Left: U-235 yield
axes[0].semilogy(A_range, yields_U235, color='#ef4444', linewidth=2.5, label='U-235 (thermal)')
axes[0].semilogy(A_range, yields_Pu239, color='#3b82f6', linewidth=2.5, label='Pu-239 (thermal)')
axes[0].set_xlabel('Fragment Mass Number A', color='#94a3b8', fontsize=12)
axes[0].set_ylabel('Fission Yield [%]', color='#94a3b8', fontsize=12)
axes[0].set_title('Fission Product Mass Distribution', color='white', fontsize=14)
axes[0].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[0].tick_params(colors='#94a3b8')
axes[0].set_ylim(1e-3, 10)
axes[0].grid(True, alpha=0.15, color='#334155')

# Annotate key products
for A, name in [(90, 'Sr-90'), (95, 'Mo-95'), (131, 'I-131'), (137, 'Cs-137'), (140, 'Ba-140')]:
    y = fission_yield(A)
    if y > 0.01:
        axes[0].annotate(name, (A, y), textcoords="offset points", xytext=(5, 10),
                         color='#fbbf24', fontsize=8)

# Right: Neutron multiplicity
E_incident = np.linspace(0, 20, 100)  # MeV
nu_bar_U235 = 2.43 + 0.065 * (E_incident - 0.025)  # approximate
nu_bar_Pu239 = 2.88 + 0.065 * (E_incident - 0.025)
nu_bar_U238 = 2.48 + 0.065 * (E_incident - 1.0)

axes[1].plot(E_incident, nu_bar_U235, color='#ef4444', linewidth=2.5, label='U-235')
axes[1].plot(E_incident, nu_bar_Pu239, color='#3b82f6', linewidth=2.5, label='Pu-239')
axes[1].plot(E_incident, nu_bar_U238, color='#22c55e', linewidth=2.5, label='U-238')

axes[1].set_xlabel('Incident Neutron Energy [MeV]', color='#94a3b8', fontsize=12)
axes[1].set_ylabel('Mean Neutrons per Fission', color='#94a3b8', fontsize=12)
axes[1].set_title('Average Neutron Multiplicity', color='white', fontsize=14)
axes[1].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[1].tick_params(colors='#94a3b8')
axes[1].grid(True, alpha=0.15, color='#334155')

plt.tight_layout()
plt.savefig('fission_yield.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())
print("U-235 thermal fission: nu_bar = 2.43 neutrons/fission")
print("Energy per fission: ~200 MeV total")
print("  Fission fragments KE: ~170 MeV")
print("  Prompt neutrons: ~5 MeV")
print("  Prompt gammas: ~7 MeV")
print("  Beta/gamma from products: ~18 MeV")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation

Criticality calculation: computes k-effective using the four-factor formula.

Reactor Criticality Calculation

Fortran

Four-factor formula k-effective calculation for various enrichment levels

criticality_calc.f9068 lines

program criticality_calc
  implicit none
  integer, parameter :: dp = selected_real_kind(15)

real(dp) :: eta, f, p, epsilon_ff, k_inf, k_eff
  real(dp) :: L_th, L_f, B2, enrichment
  real(dp) :: sigma_f_235, sigma_a_235, sigma_a_238, sigma_a_mod
  real(dp) :: N_235, N_238, N_mod, nu_bar
  integer :: i

write(*,'(A)') '=== Nuclear Reactor Criticality Calculator ==='
  write(*,'(A)') '(Four-factor formula with leakage correction)'
  write(*,'(A)') ''

! Scan enrichment levels
  write(*,'(A)') 'Enrichment(%)  eta    f      p      eps    k_inf   k_eff'
  write(*,'(A)') '-------------  ---    -      -      ---    -----   -----'

do i = 1, 8
    enrichment = 0.5_dp + 0.5_dp * i  ! 1% to 4.5%

! Cross sections at thermal energy (barns)
    sigma_f_235 = 583.0_dp    ! U-235 fission
    sigma_a_235 = 680.0_dp    ! U-235 absorption (fission + capture)
    sigma_a_238 = 2.68_dp     ! U-238 absorption
    sigma_a_mod = 0.664_dp    ! H2O moderator absorption per molecule

nu_bar = 2.43_dp  ! neutrons per fission

! Number densities (relative)
    N_235 = enrichment / 100.0_dp
    N_238 = 1.0_dp - N_235
    N_mod = 40.0_dp  ! moderator-to-fuel ratio

! eta: neutrons produced per absorption in fuel
    eta = nu_bar * sigma_f_235 * N_235 / (sigma_a_235 * N_235 + sigma_a_238 * N_238)

! f: thermal utilization (fraction absorbed in fuel)
    f = (sigma_a_235 * N_235 + sigma_a_238 * N_238) / &
        (sigma_a_235 * N_235 + sigma_a_238 * N_238 + sigma_a_mod * N_mod)

! p: resonance escape probability (approximate)
    p = exp(-N_238 * 25.0_dp / (N_mod * 1.0_dp + N_238 * 0.1_dp))

! epsilon: fast fission factor
    epsilon_ff = 1.0_dp + 0.03_dp * N_238

k_inf = eta * f * p * epsilon_ff

! Leakage correction: k_eff = k_inf / (1 + L^2 * B^2)
    ! Thermal diffusion length ~ 2.85 cm for H2O
    ! Fast diffusion length ~ 5.0 cm
    L_th = 2.85_dp  ! cm
    L_f = 5.0_dp     ! cm
    B2 = (3.14159_dp / 100.0_dp)**2  ! geometric buckling for R~100 cm sphere

k_eff = k_inf / (1.0_dp + (L_th**2 + L_f**2) * B2)

write(*,'(F10.1, 5F8.4, F9.4)') enrichment, eta, f, p, epsilon_ff, k_inf, k_eff
  end do

write(*,'(A)') ''
  write(*,'(A)') 'For k_eff = 1.0 (critical), need ~3-5% enrichment'
  write(*,'(A)') 'in a light-water reactor with H2O moderator.'
  write(*,'(A)') ''
  write(*,'(A)') 'Natural uranium (0.7% U-235) cannot sustain a chain'
  write(*,'(A)') 'reaction in H2O, but can in heavy water (D2O) or graphite.'
end program criticality_calc

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Cross Sections Fusion →