Lecture 12 of the MIT course 22.01: Introduction to Nuclear Engineering and Ionizing Radiation covers the manufacture of \(\mathrm{^{60}Co}\) by neutron irradiation of \(\mathrm{^{59}Co}\) in a research reactor. The thermal neutron absorption cross section of \(\mathrm{^{59}Co}\) is \(\sigma_0 \approx 20 \; \mathrm{barn}\), so for a neutron flux \(\Phi\), the rate of change of the number of \(\mathrm{^{59}Co}\) nuclei, \(N_0\), may be written
The \(\mathrm{^{60}Co}\) nuclei produced in this process undergo natural beta-decay with a half-life \(t_{1/2} = 1925.4\;\mathrm{days}\) to the stable isotope, \(\mathrm{^{60}Ni}\), and, potentially, neutron absorption with a cross section \(\sigma_1\). Denoting the number of \(\mathrm{^{60}Ni}\) as \(N_1\):
where \(\lambda_1 = \ln 2 / t_{1/2}\). This differential equation can be solved with the introduction of the integrating factor \(\exp[(\sigma_1\Phi + \lambda_1)t]\) since:
(There appears to be a sign error in the videoed lecture for this derivation).
In principle, the \(\mathrm{^{60}Ni}\) could capture a neutron as well, but the cross section for this process is rather small.
This Jupyter Notebook calculates and plots the amounts of these isotopes as a function of time in a reactor with isotropic neutron flux, \(\Phi = 10^{14}\;\mathrm{n\,cm^{-2}\,s^{-1}}\), and determines the maximum amount of \(\mathrm{^{60}Co}\) created from an initial amount of 100 g of \(\mathrm{^{59}Co}\). \(\sigma_1\) is small and taken to be 0.
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt
from scipy.constants import N_A
# Neutron absorption cross section for 59Co, cm2.
sigma0 = 20 * 1.e-24
# Neutron absorption cross section for 60Co.
sigma1 = 0
# Half-life (days) and decay constant for 60Co, s-1.
thalf1 = 1925.4
lam1 = np.log(2) / (thalf1 * 24 * 60 * 60)
# Neutron flux, cm2.s-1.
Phi = 1.e14
# Mass of 59Co sample, g.
m59Co = 100
# Atomic mass of 59Co, Da.
M59Co = 58.9
# Number of 59Co atoms
N00 = m59Co / M59Co * N_A
def N0(t):
"""Amount of 59Co as a function of time."""
return N00 * np.exp(-sigma0 * Phi * t)
def N1(t):
"""Amount of 60Co as a function of time."""
d = Phi * (sigma1 - sigma0) + lam1
return sigma0 * Phi * N00 / d * (np.exp(-sigma0*Phi*t) - np.exp(-(sigma1*Phi + lam1)*t))
def N2(t):
"""Amount of 60Ni as a function of time."""
return N00 - N0(t) - N1(t)
def dN1dt(t):
"""Rate of change of 60Co as a function of time."""
return sigma0 * Phi * N0(t) - N1(t) * (sigma1 * Phi + lam1)
# Plot the isotope amounts up to tmax seconds.
tmax = 3.e9
t = np.linspace(0, tmax, 1000)
plt.plot(t, N0(t), label='$\mathrm{^{59}Co}$')
plt.plot(t, N1(t), label='$\mathrm{^{60}Co}$')
plt.plot(t, N2(t), label='$\mathrm{^{60}Ni}$')
plt.xlabel('$t$ /s')
plt.ylabel('$N$ /atoms')
The problem set associated with the course also suggests finding the optimum amount of time to irradiate the \(\mathrm{^{59}Co}\) to maximise the profit from the \(\mathrm{^{60}Co}\) generated, assuming a reactor time cost of USD 1000 per day and a \(\mathrm{^{60}Co}\) price of USD 100 per \(\mathrm{\mu Ci}\).
# Cost of reactor, USD.s-1.
CR = 1.e3 / 24 / 60 / 60
# Price of 60Co, Bq-1.
CI = 100 * 1.e6 / 3.7e10
def I(t):
"""The predicted net income from the 60Co after time t in the reactor (USD)."""
# Calculate the activity of 60Co as A1 = lam1.N1.
return -CR * t + CI * lam1 * N1(t)
def dIdt(t):
"""dI/dt in USD/sec."""
return -CR + CI * lam1 * dN1dt(t)
The maximum profit can be found as the minimum of \(-I\).
from scipy.optimize import minimize
minimize(lambda t: -I(t), 0.5e9, jac = lambda t:-dIdt(t))
fun: -2806888567702.389
hess_inv: array([[42765.86838808]])
jac: array([-7.46765967e-06])
message: 'Optimization terminated successfully.'
nfev: 21
nit: 6
njev: 21
status: 0
success: True
x: array([3.38753637e+08])
This turns out to be the same time as would produce the maximum amount of \(\mathrm{^{60}Co}\) possible, presumably because USD 1000 per day is an unrealistically cheap cost for the reactor time during production and omits other considerations (such as isotopic separation and opportunity costs).
topt = 3.388e8 / 60 / 60 / 24 / 365
print(f'{topt:.1f} yrs')
10.7 yrs