Chapter 12

Nuclear Reactors

Fission Reactor Physics

A nuclear fission reactor sustains a controlled chain reaction. The key physics involves neutron moderation, absorption, and multiplication. The reactor must maintain criticality ($k_{\text{eff}} = 1$) during steady-state operation.

Thermal Reactors

Use a moderator (water, graphite, heavy water) to slow fast fission neutrons (~2 MeV) to thermal energies (~0.025 eV), where the U-235 fission cross section is large (~583 barns). Most commercial reactors are thermal.

Fast Reactors

No moderator -- fission is sustained by fast neutrons. Can breed new fissile material (Pu-239 from U-238) and burn long-lived actinide waste. Examples: sodium-cooled fast reactors (SFRs).

Neutron Moderation

Neutrons lose energy through elastic scattering with moderator nuclei. The average logarithmic energy decrement per collision is:

$$\xi = 1 + \frac{(A-1)^2}{2A}\ln\!\left(\frac{A-1}{A+1}\right)$$

The number of collisions to thermalize from $E_0 = 2$ MeV to $E_{th} = 0.025$ eV:

$$n = \frac{\ln(E_0/E_{th})}{\xi} = \frac{18.2}{\xi}$$

Moderator	A	$\xi$	n (collisions)
Hydrogen (H)	1	1.000	18
Deuterium (D)	2	0.725	25
Carbon (C)	12	0.158	115
Uranium (U)	238	0.0084	2172

Reactor Kinetics

The time-dependent behavior of a reactor is governed by the point kinetics equations:

$$\frac{dn}{dt} = \frac{\rho - \beta}{\Lambda}\,n + \sum_{i=1}^{6} \lambda_i C_i$$

$$\frac{dC_i}{dt} = \frac{\beta_i}{\Lambda}\,n - \lambda_i C_i$$

where $n$ is the neutron population, $C_i$ are delayed neutron precursor concentrations, $\rho = (k-1)/k$ is the reactivity, $\beta \approx 0.0065$ is the delayed neutron fraction, and $\Lambda \sim 10^{-4}$ s is the prompt neutron generation time.

Delayed neutrons are crucial for reactor control. Without them, the reactor period would be $\Lambda/\rho \sim 0.03$ s, far too fast for mechanical control systems. Delayed neutrons effectively lengthen the generation time to ~0.1 s, making the reactor controllable.

Python Simulation: Reactor Kinetics

Numerical solution of the point kinetics equations with six delayed neutron groups, showing reactor response to various reactivity insertions.

Reactor Neutron Population Dynamics

Python

Point kinetics equation solver showing subcritical, critical, and supercritical behavior

script.py102 lines

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

# ──────────────────────────────────────────────
# Reactor Neutron Population Dynamics
# Point kinetics with delayed neutrons
# ──────────────────────────────────────────────

# Point kinetics equations:
# dn/dt = (rho - beta)/Lambda * n + sum(lambda_i * C_i) + S
# dC_i/dt = beta_i/Lambda * n - lambda_i * C_i

# Delayed neutron parameters for U-235 thermal fission
beta_total = 0.0065  # total delayed neutron fraction
Lambda = 1.0e-4      # prompt neutron generation time (s) for thermal reactor
beta_groups = np.array([0.000215, 0.001424, 0.001274, 0.002568, 0.000748, 0.000273])
lambda_groups = np.array([0.0124, 0.0305, 0.111, 0.301, 1.14, 3.01])  # s^-1

def solve_point_kinetics(rho, t_max, dt=0.001, n0=1.0):
    """Solve point kinetics equations with 6 delayed neutron groups"""
    N_steps = int(t_max / dt)
    t = np.zeros(N_steps)
    n = np.zeros(N_steps)
    C = np.zeros((6, N_steps))

n[0] = n0
    # Initial equilibrium precursor concentrations
    for i in range(6):
        C[i, 0] = beta_groups[i] * n0 / (lambda_groups[i] * Lambda)

for k in range(N_steps - 1):
        t[k+1] = t[k] + dt

# Apply reactivity (can be time-dependent)
        rho_k = rho(t[k]) if callable(rho) else rho

dn = ((rho_k - beta_total) / Lambda * n[k] + np.sum(lambda_groups * C[:, k])) * dt
        n[k+1] = n[k] + dn

for i in range(6):
            dC = (beta_groups[i] / Lambda * n[k] - lambda_groups[i] * C[i, k]) * dt
            C[i, k+1] = C[i, k] + dC

return t, n

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

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

# Case 1: Subcritical (rho < 0)
t1, n1 = solve_point_kinetics(-0.001, 10.0, dt=0.01)
axes[0,0].plot(t1, n1, color='#3b82f6', linewidth=2)
axes[0,0].set_title('Subcritical (rho = -0.001)', color='white', fontsize=12)
axes[0,0].set_xlabel('Time [s]', color='#94a3b8')
axes[0,0].set_ylabel('Relative Power n/n0', color='#94a3b8')
axes[0,0].tick_params(colors='#94a3b8')

# Case 2: Supercritical below prompt critical (0 < rho < beta)
t2, n2 = solve_point_kinetics(0.003, 30.0, dt=0.01)
axes[0,1].plot(t2, n2, color='#22c55e', linewidth=2)
axes[0,1].set_title('Supercritical (rho = 0.003 < beta)', color='white', fontsize=12)
axes[0,1].set_xlabel('Time [s]', color='#94a3b8')
axes[0,1].set_ylabel('Relative Power n/n0', color='#94a3b8')
axes[0,1].tick_params(colors='#94a3b8')

# Case 3: Prompt critical (rho = beta)
t3, n3 = solve_point_kinetics(0.0065, 0.5, dt=0.0001)
axes[1,0].plot(t3, n3, color='#f59e0b', linewidth=2)
axes[1,0].set_title('Prompt Critical (rho = beta)', color='white', fontsize=12)
axes[1,0].set_xlabel('Time [s]', color='#94a3b8')
axes[1,0].set_ylabel('Relative Power n/n0', color='#94a3b8')
axes[1,0].tick_params(colors='#94a3b8')

# Case 4: Step reactivity insertion comparison
rho_values = [0.001, 0.003, 0.005]
colors_rho = ['#3b82f6', '#22c55e', '#ef4444']
for rho_val, color in zip(rho_values, colors_rho):
    t4, n4 = solve_point_kinetics(rho_val, 20.0, dt=0.01)
    axes[1,1].semilogy(t4, n4, color=color, linewidth=2, label=f'rho = {rho_val}')

axes[1,1].set_title('Reactivity Comparison', color='white', fontsize=12)
axes[1,1].set_xlabel('Time [s]', color='#94a3b8')
axes[1,1].set_ylabel('Relative Power n/n0', color='#94a3b8')
axes[1,1].legend(fontsize=10, facecolor='#1a1a2e', edgecolor='#334155', labelcolor='white')
axes[1,1].tick_params(colors='#94a3b8')

plt.tight_layout()
plt.savefig('reactor_kinetics.png', dpi=130, bbox_inches='tight',
            facecolor=fig.get_facecolor())
print(f"Delayed neutron fraction beta = {beta_total}")
print(f"Prompt neutron lifetime Lambda = {Lambda} s")
print(f"Reactor period for rho=0.003: T ~ {Lambda/(0.003-beta_total+beta_total*0.003/0.003):.2f} s")
print(f"Without delayed neutrons, period would be: {Lambda/0.003:.4f} s")
print(f"Delayed neutrons slow reactor response by ~100x")

Click Run to execute the Python code

Code will be executed with Python 3 on the server

Fortran Implementation

Point kinetics equation solver with six delayed neutron groups.

Point Kinetics Equation Solver

Fortran

Solves reactor point kinetics with delayed neutrons for step reactivity insertion

point_kinetics.f9077 lines

program point_kinetics
  implicit none
  integer, parameter :: dp = selected_real_kind(15)
  integer, parameter :: n_groups = 6
  real(dp), parameter :: beta_total = 0.0065_dp
  real(dp), parameter :: Lambda = 1.0d-4  ! prompt generation time (s)

real(dp) :: beta_i(n_groups), lambda_i(n_groups)
  real(dp) :: n_pop, C_i(n_groups), rho, dt, t, t_max
  real(dp) :: dn, dC
  integer :: i, k, n_steps

! Delayed neutron group parameters (U-235)
  beta_i = (/ 0.000215_dp, 0.001424_dp, 0.001274_dp, &
              0.002568_dp, 0.000748_dp, 0.000273_dp /)
  lambda_i = (/ 0.0124_dp, 0.0305_dp, 0.111_dp, &
                0.301_dp, 1.14_dp, 3.01_dp /)

write(*,'(A)') '=== Point Kinetics Equation Solver ==='
  write(*,'(A,F8.4)') 'Beta (total) = ', beta_total
  write(*,'(A,ES10.3,A)') 'Lambda = ', Lambda, ' s'
  write(*,'(A)') ''

! Case: step reactivity insertion rho = 0.003
  rho = 0.003_dp
  write(*,'(A,F8.4)') 'Reactivity rho = ', rho
  write(*,'(A,F8.4)') 'rho/beta = ', rho/beta_total
  write(*,'(A)') ''

n_pop = 1.0_dp
  do i = 1, n_groups
    C_i(i) = beta_i(i) * n_pop / (lambda_i(i) * Lambda)
  end do

dt = 0.001_dp
  t_max = 20.0_dp
  n_steps = nint(t_max / dt)

write(*,'(A)') '  Time(s)    n/n0       Period(s)'
  write(*,'(A)') '  -------    ----       ---------'

do k = 0, n_steps
    t = k * dt

! Print at selected times
    if (mod(k, nint(1.0_dp/dt)) == 0) then
      write(*,'(F9.1, ES14.4)') t, n_pop
    end if

! Euler step
    dn = ((rho - beta_total) / Lambda * n_pop + &
          dot_product(lambda_i, C_i)) * dt
    n_pop = n_pop + dn

do i = 1, n_groups
      dC = (beta_i(i) / Lambda * (n_pop - dn/dt*dt) - lambda_i(i) * C_i(i)) * dt
      C_i(i) = C_i(i) + dC
    end do
  end do

write(*,'(A)') ''
  write(*,'(A)') 'Final neutron population (relative):'
  write(*,'(A,ES12.4)') '  n/n0 = ', n_pop
  write(*,'(A)') ''

! Inhour equation: reactor period calculation
  write(*,'(A)') '--- Reactor Period from Inhour Equation ---'
  write(*,'(A)') 'rho = beta_eff * sum(beta_i*Lambda_i / (1 + Lambda_i*omega))'
  write(*,'(A)') ''
  write(*,'(A)') 'For rho < beta (delayed supercritical):'
  write(*,'(A)') '  Period is dominated by longest-lived group'
  write(*,'(A,F8.2,A)') '  T ~ ', 1.0_dp/lambda_i(1), ' s (group 1 half-life ~80s)'
  write(*,'(A)') ''
  write(*,'(A)') 'For rho > beta (prompt supercritical):'
  write(*,'(A,ES10.3,A)') '  T ~ Lambda/(rho-beta) = ', Lambda/(rho - beta_total), ' s'
  write(*,'(A)') '  (Extremely fast - must be avoided!)'
end program point_kinetics

Click Run to execute the Fortran code

Code will be compiled with gfortran and executed on the server

← Collective Models Nucleosynthesis →