In a nuclear fusion reaction two atomic nuclei combine to form a single nucleus of lower total mass, the difference in mass, $\Delta m$ being released as energy in accordance with $E = \Delta m c^2$. It is this process which powers stars (in our own sun, hydrogen nuclei are fused into helium), and nuclear fusion has been actively pursued as a potential clean and cheap energy source in reactors on Earth for over 50 years.
In deciding which nuclei should be the "fuel" in a nuclear fusion reactor, an important quantity is the "cross section", $\sigma$, a measure of the reaction probability as a function of the kinetic energy of the reactant nuclei. The nuclei are positively-charged and so repel each other (Coulomb repulsion): fusion is only possible if they are able to overcome this repulsion to approach each other at close range. This, in turn, demands that they have a large relative speed and hence kinetic energy. For a viable fusion reactor on Earth, this corresponds to a temperature of many millions of Kelvin.
The IAEA's Nuclear Data Section provides a database of evaluated nuclear data, ENDF, that can be used to illustrate the fusion cross section. The code below uses data from ENDF to plot the fusion cross section for three fusion reactions:
where the D + D reaction has two possible channels with approximately equal probability of occuring; the orange line below depicts the total cross section for both processes.
The reason that the majority of research on experimental fusion reactors has focused on the D + T process is that it has the largest cross section, which peaks at the lowest temperature. This means that the plasma in which fusion occurs does not have to be heated as much: indeed, ITER aims to operate at 150 million K. The 17.6 MeV released in the D + T reaction corresponds to an energy release many times greater than a typical chemical reaction (the combustion of octane releases only 57 eV, which is more than 300,000 times less energy).
The code below requires the following data files:
A small detail arises in the transformation of energy from the "lab-fixed" frame of the ENDF format (in which the heavier, target, reactant is assumed to be stationary) to a frame moving with the centre-of-mass of the system. This is handled in the function which reads in the cross section files if the constant COFM
is set to True
. The cross sections provided by ENDF are not all on the same energy grid, so they are also interpolated onto a common grid by this function.
import numpy as np
from matplotlib import rc
import matplotlib.pyplot as plt
from scipy.constants import e, k as kB
rc('font', **{'family': 'serif', 'serif': ['Computer Modern'], 'size': 14})
rc('text', usetex=True)
# To plot using centre-of-mass energies instead of lab-fixed energies, set True
COFM = True
# Reactant masses in atomic mass units (u).
masses = {'D': 2.014, 'T': 3.016, '3He': 3.016}
# Energy grid, 1 – 1000 keV, evenly spaced in log-space.
Egrid = np.logspace(0, 3, 100)
def read_xsec(filename):
"""Read in cross section from filename and interpolate to energy grid."""
E, xs = np.genfromtxt(filename, comments='#', skip_footer=2, unpack=True)
if COFM:
collider, target = filename.split('_')[:2]
m1, m2 = masses[target], masses[collider]
E *= m1 / (m1 + m2)
xs = np.interp(Egrid, E*1.e3, xs*1.e-28)
return xs
# D + T -> α + n
DT_xs = read_xsec('D_T_-_a_n.txt')
# D + D -> T + p
DDa_xs = read_xsec('D_D_-_T_p.txt')
# D + D -> 3He + n
DDb_xs = read_xsec('D_D_-_3He_n.txt')
# Total D + D fusion cross section is due to equal contributions from the
# above two processes.
DD_xs = DDa_xs + DDb_xs
# D + 3He -> α + p
DHe_xs = read_xsec('D_3He_-_4He_p.txt')
fig, ax = plt.subplots()
ax.loglog(Egrid, DT_xs, lw=2, label='$\mathrm{D-T}$')
ax.loglog(Egrid, DD_xs, lw=2, label='$\mathrm{D-D}$')
ax.loglog(Egrid, DHe_xs, lw=2, label='$\mathrm{D-^3He}$')
ax.grid(True, which='both', ls='-')
ax.set_xlim(1, 1000)
xticks= np.array([1, 10, 100, 1000])
ax.set_xticks(xticks)
ax.set_xticklabels([str(x) for x in xticks])
if COFM:
xlabel ='E(CofM) /keV'
else:
xlabel ='E /keV'
ax.set_xlabel(xlabel)
ax.set_ylabel('$\sigma\;/\mathrm{m^2}$')
ax.set_ylim(1.e-32, 1.e-27)
ax.legend()
plt.savefig('fusion-xsecs.png')
plt.show()
Comments
Comments are pre-moderated. Please be patient and your comment will appear soon.
Mark Paris 6 years, 1 month ago
Thanks for this useful and informative post -- much appreciated!
Link | Replychristian 6 years, 1 month ago
Glad you found it useful – thanks for the feedback!
Link | ReplyJohn Sharma 5 years, 5 months ago
Excellent code and textbook. Thanks. John.
Link | ReplyGeorg Harrer 4 years, 7 months ago
Super helpful thanks.
Link | ReplyI'd like to include the p+p -> d + positron + neutrino reaction in the plot allthough the cross section is soo much lower (i have an old plot from lecture nots but with no data source unfortunately)
The problem is that i'm not able to get the p + p cross section data from the ENDF database.
Do you maybe have any idea where or how to get it?
christian 4 years, 7 months ago
Yes, I suppose it might not be in ENDF. Have you looked in the solar physics literature, e.g. the papers by Adelberger et al. such as Rev. Mod. Phys. 70 1265 (1998)?
Link | ReplyMiguel 4 years, 3 months ago
Thanks Christian for the code!
Link | ReplyI have found some differences between the ENDF cross-sections of the D+D reaction and other ENDF reported in the literature
https://doi.org/10.1017/S026303460404011X
How have you processed them or where have you found them?
Kind regards,
Miguel
christian 4 years, 3 months ago
Hi Miguel – I did not process the data any further than described in the article (transformation from lab-fixed to CofM coordinates). The data used come from the IAEA's ENDF database (linked provided). I can't access the article you cite, I'm afraid.
Link | ReplyJarmo Kanerva 4 years, 2 months ago
Just participating a university introductory plasma physics course. Had a clue to answer "why to use D+T in a fusion reactor", and found now the quite complete basic explanation for that question! Thank you!
Link | Replychristian 4 years, 2 months ago
I'm glad it was helpful – good luck with your course!
Link | ReplyVasundhara 3 years, 9 months ago
where to find the data files for other fusion reactions like 3 Helium-3Helium fusion
Link | Replychristian 3 years, 9 months ago
As indicated in the post, you can get this sort of data from the IAEA's ENDF database: https://www-nds.iaea.org/exfor/endf.htm – in the case of 3He-3He fusion, search for He-3 and select the reaction HE-3(HE3,2P),SIG.
Link | ReplyJaves 3 years, 9 months ago
Thanks for the post and one question is where can I find the data on averaged reaction rate as <\sigma v> ?
Link | Replychristian 3 years, 9 months ago
I think you probably have to calculate the rate numerically as the integral of the (normalized) velocity distribution times the cross section.
Link | ReplyFrank 3 years, 4 months ago
hi christian, nice, very useful.
Link | Replydid you also provide the code to plot the reactivities somewhere?
thanks.
fw
christian 3 years, 4 months ago
Thanks, Frank – is this what you're looking for? https://scipython.com/blog/nuclear-fusion-cross-sections/ – I've added a link between these posts now.
Link | ReplyCheers, Christian
Sarah 3 years, 2 months ago
Hi! How would you obtain a relationship for the fusion probability?
Link | ReplyNew Comment