Nuclear fusion cross sections and reactivities

(0 comments)

A quick update on this blog article with some extra cross sections and an additional plot of Maxwell-averaged reactivities. The processes covered are:

\begin{align*} \mathrm{D} + \mathrm{T} & \rightarrow \alpha (3.5\;\mathrm{MeV}) + n (14.1\;\mathrm{MeV}) & \\ \mathrm{D} + \mathrm{D} & \rightarrow \mathrm{T} (1.01\;\mathrm{MeV}) + p(3.02\;\mathrm{MeV})& 50\%\\ & \rightarrow \mathrm{^3He} (0.82\;\mathrm{MeV}) + n(2.45\;\mathrm{MeV}) & 50\%\\ \mathrm{D} + \mathrm{^3He} & \rightarrow \alpha (3.6\;\mathrm{MeV}) + p(14.7\;\mathrm{MeV}) & \\ \mathrm{^3He} + \mathrm{^3He} & \rightarrow 2p + \alpha + 14.7\;\mathrm{MeV} & \\ \mathrm{T} + \mathrm{^3He} & \rightarrow n + p + \alpha + 12.1\;\mathrm{MeV} & 59\% \\ & \rightarrow D (9.52\;\mathrm{MeV}) + \alpha(4.8\;\mathrm{MeV}) & 41\% \\ \mathrm{T} + \mathrm{T} & \rightarrow 2n + \alpha + 11.3\;\mathrm{MeV} & \\ p + \mathrm{^{11}B} & \rightarrow 3\alpha + 8.7\;\mathrm{MeV} & \end{align*}

Fusion cross sections

Fusion reactivities

The cross section gives a probability for the reaction to occur as a function of the collision energy. In practice, the colliding particles move with a distribution of relative velocities, $f(v)$ (and hence collision energies); a useful quantity is the averaged reactivity: $$ \langle \sigma v \rangle = \int_0^\infty \sigma(v)vf(v)\,\mathrm{d}v. $$ For the Maxwell equilibrium velocity distribution, $$ f_i(v_i) = \left(\frac{m_i}{2\pi k_\mathrm{B}T}\right)^{3/2}\exp\left(-\frac{m_iv_i^2}{2k_\mathrm{B}T}\right), $$ it can be shown [1, 2] that $$ \langle \sigma v \rangle = \frac{4\pi}{\sqrt{2\pi\mu}}\frac{1}{(k_\mathrm{B}T)^{3/2}} \int_0^\infty \sigma(E)E\exp\left(-\frac{E}{k_\mathrm{B}T}\right)\,\mathrm{d}E, $$ where $\mu = m_1m_2/(m_1+m_2)$ is the reduced mass of the colliding particles.

The code below requires the following data files:

The code:

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': 18})
rc('text', usetex=True)

# Plot resolution (dpi) and figure size (pixels)
DPI, FIG_WIDTH, FIG_HEIGHT = 72, 800, 600

# 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, '11B': 11.009305167,
          'p': 1.007276466620409}

# Energy grid, 1 – 1000 keV, evenly spaced in log-space.
Egrid = np.logspace(0, 3, 500)

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')

# p + 11B -> 3α
pB_xs = read_xsec('p_11B_-_3a.txt')

# T + T -> 4He + 2n
TT_xs = read_xsec('T_T_-_4He_n_n.txt')

# T + 3He -> 4He + n + p
THe_xs = read_xsec('T_3He_-_n_p_4He.txt')

# 3He + 3He -> 4He + 2p
HeHe_xs = read_xsec('3He_3He_-_p_p_4He.txt')


fig = plt.figure(figsize=(FIG_WIDTH/DPI, FIG_HEIGHT/DPI), dpi=DPI)
ax = fig.add_subplot(111)

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.loglog(Egrid, pB_xs, lw=2, label='$\mathrm{p-^{11}B}$')
ax.loglog(Egrid, TT_xs, lw=2, label='$\mathrm{T-T}$')
ax.loglog(Egrid, THe_xs, lw=2, label='$\mathrm{T-^3He}$')
ax.loglog(Egrid, HeHe_xs, lw=2, label='$\mathrm{^3He-^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(CM) /keV'
else:
    xlabel ='E /keV'
ax.set_xlabel(xlabel)

# A second x-axis for energies as temperatures in millions of K
ax2 = ax.twiny()
ax2.set_xscale('log')
ax2.set_xlim(1,1000)
xticks2 = np.array([15, 150, 1500, 5000])
ax2.set_xticks(xticks2 * kB/e * 1.e3)
ax2.set_xticklabels(xticks2)
ax2.set_xlabel('$T$ /million K')

ax.set_ylabel('$\sigma\;/\mathrm{m^2}$')
ax.set_ylim(1.e-32, 1.e-27)

ax.legend(loc='upper left')
plt.savefig('fusion-xsecs2.png', dpi=DPI)
plt.show()
  1. D. D. Clayton, "Principles of Stellar Evolution and Nucleosynthesis", University of Chicago Press (1983).
  2. S. Atzeni and J. Meyer-ter-Vehn, "The Physics of Inertial Fusion", OUP (2004).
Currently unrated

Comments

Comments are pre-moderated. Please be patient and your comment will appear soon.

There are currently no comments

New Comment

required

required (not published)

optional

required