WarpX
BoschHaleFusionCrossSection.H
Go to the documentation of this file.
1 /* Copyright 2022 Remi Lehe
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 
8 #ifndef WARPX_BOSCH_HALE_FUSION_CROSS_SECTION_H
9 #define WARPX_BOSCH_HALE_FUSION_CROSS_SECTION_H
10 
12 #include "Utils/WarpXConst.H"
13 
14 #include <AMReX_REAL.H>
15 
16 #include <cmath>
17 
29 amrex::ParticleReal BoschHaleFusionCrossSection (
30  const amrex::ParticleReal& E_kin_star,
31  const NuclearFusionType& fusion_type,
32  const amrex::ParticleReal& m1,
33  const amrex::ParticleReal& m2 )
34 {
35  using namespace amrex::literals;
36 
37  constexpr amrex::ParticleReal joule_to_keV = 1.e-3_prt/PhysConst::q_e;
38  const amrex::ParticleReal E_keV = E_kin_star*joule_to_keV;
39 
40  // If kinetic energy is 0, return a 0 cross section and avoid later division by 0.
41  if (E_keV == 0._prt) {return 0._prt;}
42 
43  // Compute the Gamow constant B_G (in keV^{1/2})
44  // (See Eq. 3 in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
45  const amrex::ParticleReal m_reduced = m1 / (1._prt + m1/m2);
46  // The formula for `B_G` below assumes that both reactants have Z=1
47  // When one of the reactants is helium (Z=2), this formula is corrected further below.
48  amrex::ParticleReal B_G = MathConst::pi * PhysConst::alpha
49  * std::sqrt( 2._prt*m_reduced*PhysConst::c*PhysConst::c * joule_to_keV );
51  // Take into account the fact that Z=2 for one of the reactant (helium) in this case
52  B_G *= 2;
53  }
54 
55  // Compute astrophysical_factor
56  // (See Eq. 9 and Table IV in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
57  amrex::ParticleReal A1=0_prt, A2=0_prt, A3=0_prt, A4=0_prt, A5=0_prt, B1=0_prt, B2=0_prt, B3=0_prt, B4=0_prt;
59  A1 = 6.927e4_prt;
60  A2 = 7.454e8_prt;
61  A3 = 2.050e6_prt;
62  A4 = 5.2002e4_prt;
63  A5 = 0_prt;
64  B1 = 6.38e1_prt;
65  B2 = -9.95e-1_prt;
66  B3 = 6.981e-5_prt;
67  B4 = 1.728e-4_prt;
68  }
70  A1 = 5.5576e4_prt;
71  A2 = 2.1054e2_prt;
72  A3 = -3.2638e-2_prt;
73  A4 = 1.4987e-6_prt;
74  A5 = 1.8181e-10_prt;
75  B1 = 0_prt;
76  B2 = 0_prt;
77  B3 = 0_prt;
78  B4 = 0_prt;
79  }
81  A1 = 5.3701e4_prt;
82  A2 = 3.3027e2_prt;
83  A3 = -1.2706e-1_prt;
84  A4 = 2.9327e-5_prt;
85  A5 = -2.5151e-9_prt;
86  B1 = 0_prt;
87  B2 = 0_prt;
88  B3 = 0_prt;
89  B4 = 0_prt;
90  }
91  else if (fusion_type == NuclearFusionType::DeuteriumHeliumToProtonHelium) {
92  A1 = 5.7501e6_prt;
93  A2 = 2.5226e3_prt;
94  A3 = 4.5566e1_prt;
95  A4 = 0_prt;
96  A5 = 0_prt;
97  B1 = -3.1995e-3_prt;
98  B2 = -8.5530e-6_prt;
99  B3 = 5.9014e-8_prt;
100  B4 = 0_prt;
101  }
102 
103  const amrex::ParticleReal astrophysical_factor =
104  (A1 + E_keV*(A2 + E_keV*(A3 + E_keV*(A4 + E_keV*A5)))) /
105  (1_prt + E_keV*(B1 + E_keV*(B2 + E_keV*(B3 + E_keV*B4))));
106 
107  // Compute cross-section in SI units
108  // (See Eq. 8 in H.-S. Bosch and G.M. Hale 1992 Nucl. Fusion 32 611)
109  constexpr amrex::ParticleReal millibarn_to_sqm = 1.e-31_prt;
110  return millibarn_to_sqm * astrophysical_factor/E_keV * std::exp(-B_G/std::sqrt(E_keV));
111 }
112 
113 #endif // WARPX_BOSCH_HALE_FUSION_CROSS_SECTION_H
#define AMREX_INLINE
#define AMREX_GPU_HOST_DEVICE
NuclearFusionType
Definition: BinaryCollisionUtils.H:26
@ DeuteriumDeuteriumToProtonTritium
@ DeuteriumDeuteriumToNeutronHelium
AMREX_GPU_HOST_DEVICE AMREX_INLINE amrex::ParticleReal BoschHaleFusionCrossSection(const amrex::ParticleReal &E_kin_star, const NuclearFusionType &fusion_type, const amrex::ParticleReal &m1, const amrex::ParticleReal &m2)
Computes the fusion cross section, using the analytical fits given in H.-S. Bosch and G....
Definition: BoschHaleFusionCrossSection.H:29
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
static constexpr auto q_e
elementary charge [C]
Definition: constant.H:50