name: BH_ReactionRates.ipynb
date: 11 Jan 14
auth: Zach Hartwig hartwig@psfc.mit.edu
This IPython notebook plots the Bosch-Hale volumetric reaction rate parameterizations found in H.-S. and G.M. Hale, Nuclear Fusion 32 (1992) p611. Reaction energies are in [keV]; cross sections are in [m$^2$]. The B-H parameterizations have been computed for four reactions: $\mathrm{^{2}H(d,n)^{3}He}$; $\mathrm{^{2}H(d,p)^{3}H}$; $\mathrm{^{3}H(d,n)^{4}He}$; and $\mathrm{^{3}He(d,p)^{4}He}$. Each of these reactions is treated here.
The final figure produced is used in Chapter 9 of the MFE Formulary. Additionally, the code used in this notebook can easily be extended or adapted for your own calculations since the B-H coefficients and parameters are alreadyentered correctly.
In [2]:
import numpy as np
import matplotlib.pyplot as plt
In [3]:
# Assign indices for each of the four reactions
i_2Hdn = 0 # 2H(d,n)3He "D-D"
i_2Hdp = 1 # 2H(d,p)3H "D-D"
i_3Hdn = 2 # 3H(d,n)4He "D-T"
i_3Hedp = 3 # 3He(d,p)4He "D-He3"
reactions = (i_2Hdn, i_2Hdp, i_3Hdn, i_3Hedp)
# Assign numerial values for all the B-H parameters. Notation
# ("A","B","B_G") corresponds exactly as found in the original
# B-H reference for clarity.
# Allocate the Bosch-Hale coefficients for each fusion reaction
B_G_2Hdn = 31.3970
B_G_2Hdp = 31.3970
B_G_3Hdn = 34.3827
B_G_3Hedp = 68.7508
B_G = (B_G_2Hdn, B_G_2Hdp, B_G_3Hdn, B_G_3Hedp)
muc2_2Hdn = 937814
muc2_2Hdp = 937814
muc2_3Hdn = 1124656
muc2_3Hedp = 1124572
muc2 = (muc2_2Hdn, muc2_2Hdp, muc2_3Hdn, muc2_3Hedp)
# C_1 C_2 C_3 C_4 C_5 C_6 C_7
C_2Hdn = np.array([5.43360e-12, 5.85778e-3, 7.68222e-3, 0.0 , -2.96400e-6, 0.0 , 0.0])
C_2Hdp = np.array([5.65718e-12, 3.41267e-3, 1.99167e-3, 0.0 , 1.05060e-5, 0.0 , 0.0])
C_3Hdn = np.array([1.17302e-9 , 1.51361e-2, 7.51886e-2, 4.60643e-3, 1.35000e-2, -1.06750e-4, 1.36600e-5])
C_3Hedp = np.array([5.51036e-10, 6.41918e-3, -2.02896e-3, -1.91080e-5, 1.35776e-4, 0.0 , 0.0])
C = (C_2Hdn, C_2Hdp, C_3Hdn, C_3Hedp)
# Assign temperature ranges for each reaction [E_min, E_max] [keV]
Tmp_2Hdn = np.array([0.5, 4900]) #
Tmp_2Hdp = np.array([0.5, 5000])
Tmp_3Hdn = np.array([0.5, 550])
Tmp_3Hedp = np.array([0.3, 900])
Tmp = (Tmp_2Hdn, Tmp_2Hdp, Tmp_3Hdn, Tmp_3Hedp)
In [4]:
# Compute the B-H volumetric reaction rates. The rates will be arrays
# corresponding to rate [1 / m3 / s] as functions of temperature [keV]
# List to hold the reaction rate arrays for each reaction
sigma_v = []
# Iterate over each of the four reactions
for r in reactions:
T = np.logspace(-1, 2, 200) # [keV]
#T = np.arange(1,20,1.0)
theta_p0 = T * (C[r][1] + T * (C[r][3] + T * C[r][5]))
theta_p1 = 1 + T * (C[r][2] + T * (C[r][4] + T * C[r][6]))
theta = T / (1 - theta_p0 / theta_p1)
xi = (B_G[r]**2 / 4 / theta)**(1.0 / 3)
sv_p0 = sqrt(xi / (muc2[r] * T**3))
sv_p1 = exp(-3 * xi)
sv = C[r][0] * theta * sv_p0 * sv_p1 # [cm^-3 s^-1]
# Transform units from [cm^-3 s^-1] to [m^-3 s^-1]
cm3_2_m3 = 1e-6
sv = sv * cm3_2_m3 # [cm^-3 s^-1]
# Append this reaction cross section to container
sigma_v.append(sv) # [m2]
In [5]:
# Draw all four reaction rates on the same log-log plot over
# reasonably useful temperature ranges
fig, ax = plt.subplots()
plt.hold(True)
ax.loglog(T, sigma_v[0], '-bo', linewidth=2, markersize=14, markevery=30, label='$\mathrm{^2\!H(d,n)^3\!He}$')
ax.loglog(T, sigma_v[1], '-k^', linewidth=2, markersize=12, markevery=30, label='$\mathrm{^2\!H(d,p)^3\!H}$')
ax.loglog(T, sigma_v[2], '-rs', linewidth=2, markersize=11, markevery=30, label='$\mathrm{^3\!H(d,n)^4\!He}$')
ax.loglog(T, sigma_v[3], '-g*', linewidth=2, markersize=15, markevery=30, label= '$\mathrm{^3\!He(d,p)^4\!He}$')
ax.hold(False)
ax.legend(loc='lower right', fontsize=24, numpoints=1)
font = {'family' : 'sans-serif',
'style' : 'normal',
'weight' : 'normal',
'size' : 28}
matplotlib.rc('font', **font)
ax.grid(axis='both')
plt.title('Bosch-Hale Volumetric Reaction Rates', fontsize=24)
plt.xlabel('Temperature [keV]')
ax.xaxis.set_tick_params(which='major', width=1, length=10)
ax.xaxis.set_tick_params(which='minor', width=1, length=5)
ax.yaxis.set_tick_params(which='major', width=1, length=10)
ax.yaxis.set_tick_params(which='minor', width=1, length=5)
plt.ylabel('Volumetric reaction rate [ m$^{-3}$ s$^{-1}$]')
plt.yticks([1e-32, 1e-30, 1e-28, 1e-26, 1e-24, 1e-22, 1e-20])
pylab.xlim([0.3, 100])
pylab.ylim([1e-32,1e-20])
gcf().subplots_adjust(top=0.92, bottom=0.2, left=0.2, right=0.95)
fig.set_size_inches(10,8)
plt.show()
In [6]:
# Save the figure as a vector-graphic EPS file
fig.savefig('BH_ReactionRates.eps')
In [7]:
sigma_v[0]
Out[7]:
In [ ]: