diff --git a/documentation/source/physics-models/fusion_reactions/beam_reactions.md b/documentation/source/physics-models/fusion_reactions/beam_reactions.md index 7b12c755a8..9d77fb25d8 100644 --- a/documentation/source/physics-models/fusion_reactions/beam_reactions.md +++ b/documentation/source/physics-models/fusion_reactions/beam_reactions.md @@ -1,392 +1,860 @@ -# Neutral Beam Injection Fusion +# Beam Fusion Model -The main function called for calculating the fusion reactions produced by neutral beam injection is `beam_fusion()` -Due to the small contribution of fusion power from the neutral beams only D-T reactions are taken into account, as D-D additions to fusion power are deemed to be negligible. -The beam fusion calculations will only run if the calculated beam current is greater than 0. This is done by having a NBI heating and current drive configuration. +This page describes the neutral beam fusion model currently implemented in `beam_fusion()`. -The NBI parameters taken from the current drive module to be used in the beam fusion calculations are the beam current (`c_beam_total`), beam energy (`e_beam_kev`) and the tritium component of the beam (`f_beam_tritium`). +The model is a reduced beam-target treatment for neutral beam ions injected into the plasma. `PROCESS` first determines an effective beam current that already accounts for beam transport effects such as shine-through losses. The beam-fusion model then uses this net beam current to build a steady population of fast ions and estimates how many of those ions fuse with the background plasma. `beam_fusion()` does not explicitly model beam attenuation, orbit effects, or spatial beam evolution, and instead treats the fast-ion population in a volume-averaged (0D) sense. The model is calculated in the following steps: -Please see the [H&CD section](../../eng-models/heating_and_current_drive/heating-and-current-drive.md) of the docs for more info. - ------------------------- - -## Beam slowing down properties | `beam_fusion()` - -1. **Calculate the beam ion slowing down time given by**: - - $$ - \tau_{\text{slow}} = 1.99 \times 10^{19}\left(A_{\text{D}}\left(1.0-f_{\text{tritium-beam}}\right)+(A_{\text{T}}f_{\text{tritium-beam}})\right)\frac{\langle T_{\text{e}}\rangle^{3/2}}{\langle n_{\text{e}} \rangle \Lambda_{\text{ie}}} - $$ - -2. **Calculate the the beam critical energy** - - The alpha particles are born with an energy of 3.5 MeV and initially slow down mainly by collisions with electrons. At a critical energy $E_{\text{crit}}$ the rate of loss to the ions becomes equal to that to the electrons, and at lower energies the loss to the ions predominates.[^1] - - $$ - E_{\text{crit}} = 14.8 A T_{\text{e}}\left[\frac{1}{n_{\text{e} \ln{ \Lambda_{\text{e}}}}}\left[\Sigma \frac{n_j Z_j^2\ln{\Lambda_{\text{i}}}}{A_j}\right]\right]^{2/3} \ [\text{eV}] - $$ - - This can be approximated to: - - $$ - E_{\text{crit}} \approx 0.1AT_{10} \ \ [\text{MeV}] - $$ - - Though is currently implemented for deuterium as: - - $$ - E_{\text{crit,D}} = 14.8A_{\text{D}}T_{\text{e}}Z_{\text{eff}}^{2/3}\frac{\ln \Lambda_{\text{ie}}+4.0}{\Lambda_{\text{ie}}} - $$ - - The tritium critical energy is simply just scaled with the ratio of atomic mass numbers - - $$ - E_{\text{crit,T}} = E_{\text{crit,D}}\left(\frac{A_{\text{T}}}{A_{\text{D}}}\right) - $$ - - --------------------------- - - ### Derivation of beam slowing down rate and critical energy - - The rate of slowing down of a test particle of mass $M$, charge $Z\text{e}$ and energy $E$, due to Coulomb collisions with a background species off mass $m_{\text{j}}$, charge $Z_{\text{j}}\text{e}$, density $n_{\text{j}}$ and temperature $T_{\text{j}}$ is given by[^2]: - - $$ - \frac{\mathrm{dE}}{\mathrm{dt}}=\left[-\Phi\left(x_{\text{j}}\right)+x_{\text{j}}\left(1+\frac{m_{\text{j}}}{M} \Phi^{\prime}\left(x_{\text{j}}\right)\right)\right] \frac{4 \pi n_{\text{j}}}{m_{\text{j}} V}\left(\frac{Z z_{\text{j}} Z \text{e}^2}{4 \pi \varepsilon_0}\right)^2 \ln \Lambda_{\text{j}}, - $$ - - where $\Phi(x)$ is the error function, +1. computes a beam slowing-down time and a critical energy for deuterium and tritium beam ions, +2. estimates the steady-state hot beam ion densities from the beam source rate and slowing-down residence time, +3. evaluates a fast-ion pressure and a pressure-equivalent deposited beam energy, +4. computes effective beam-target fusion rate coefficients for hot deuterium and hot tritium beam ions, +5. forms total beam-target reaction rates and converts them to alpha power, +6. computes the neutral beam beta contribution from the hot beam population. - $$ - \Phi^{\prime}(x) = \frac{\mathrm{d\Phi}}{\mathrm{dx}}, \ \ V = \sqrt{\frac{2E}{M}}, \ \ V_{\text{j}} = \sqrt{\frac{2kT_{\text{j}}}{m_{\text{j}}}}, \ \ x_{\text{j}} = \frac{V}{V_{\text{j}}} - $$ +The main functions involved are: - and $\ln \Lambda_{\text{j}}$ is the usual Coulomb logarithm for the test particle and background species interactions. +- `beam_fusion()` +Top-level routine. Computes beam slowing-down properties, beam-target fusion alpha power, and neutral beam beta. - For the contribution to the slowing down due to interaction with the background ions, we can use a large argument expansion for the error function. This is because the fusion born ions have a velocity $V$, much greater than the thermal velocity $V_{\text{j}}$, of the background ions. The velocity of the fast ions is much less than the thermal velocity of electrons, however. For the electron contribution to the slowing down we use the small argument expansion for the error function. The net slowing down rate is then given by - $$ - \frac{\mathrm{d E}}{\mathrm{d t}}=-\frac{A Z^2 \sqrt{M}}{\sqrt{E}}-\frac{B Z^2 E}{M} - $$ - where the coefficients $A$ and $B$ are given by - $$ - \begin{aligned} - & A=\frac{4 \pi}{\sqrt{2}}\left(\frac{\text{e}^2}{4 \pi \varepsilon_0}\right)^2 \sum_j\left(\frac{n_{\text{j}} Z_{\text{j}}^2}{m_{\text{j}}} \ln \Lambda_j\right) \\ - & B=\frac{16 \sqrt{\pi}}{3 k T_{\text{e}}} \sqrt{\frac{m_{\text{e}}}{2 k T_{\text{e}}}}\left(\frac{\text{e}^2}{4 \pi \varepsilon_0}\right)^2 n_{\text{e}} \ln \Lambda_{\text{e}} - \end{aligned} - $$ +- `beam_slowing_down_state()` +Splits the beam into deuterium and tritium components, computes steady-state hot beam densities, critical speeds, fast-ion pressures, and the density-weighted deposited beam energy. - The sum over $\text{j}$ in $A$ is over the various ionic species. The quantities with subscript $\text{e}$ refer to the electrons. +- `fast_ion_pressure_integral()` +Returns the closed-form dimensionless integral factor used in the fast-ion pressure expression. - $\frac{\mathrm{d E}}{\mathrm{d t}}$ can be rewritten in the form, +- `beam_reaction_rate_coefficient()` +Evaluates an effective beam-target fusion rate coefficient for a slowing-down fast-ion population. - $$ - \frac{d E}{d t}=-\frac{2 E}{\tau_{\text{slow}}}\left[1+\left(\frac{E}{E}\right)^{3 / 2}\right] - $$ +- `_hot_beam_fusion_reaction_rate_integrand()` +Internal integrand used in the beam-target rate coefficient calculation. - The critical energy $E_{\text{c}}$ is given by +- `_beam_fusion_cross_section()` +Internal beam fusion cross-section fit used by the beam-target rate coefficient calculation. - $$ - E_c=\left[\frac{3 \sqrt{\pi}}{4} \frac{M^{3 / 2}}{n_{\text{e}} \sqrt{m}_{\text{e}}} \sum_j\left(\frac{n_{\text{j}} z_{\text{j}}^2}{m_{\text{j}}} \ln \Lambda_{\text{j}}\right) \frac{1}{\ln \Lambda_{\text{e}}}\right]^{2 / 3} k T_{\text{e}} - $$ +- `beam_target_reaction_rate()` +Converts beam density, target density, and beam-target reactivity into a total reaction rate. - Some authors take $\ln \Lambda_i=\ln \Lambda_{\text{e}}$ in the expression for $E_C$, but this is not a good approximation for fast ion slowing down in fusion plasmas since $\ln \Lambda_{\text{e}} \approx 17$, while $\ln \Lambda_i \approx 22$. The slowing down time, $\tau_{\text{slow}}$ is given by +- `alpha_power_beam()` +Converts the total beam-target reaction rate into alpha power. - $$ - \tau_{\text{slow}}=\frac{3(k T_{\text{e}})^{3 / 2}}{4 \sqrt{2 \pi m_{\text{e}}} Z^2}\left(\frac{4 \pi \varepsilon_0}{\text{e}^2}\right)^2 \frac{M}{n_{\text{e}} \ln \Lambda_{\text{e}}} - $$ +!!! info "Model assumption: beam fusion channels" - When the particle energy $E$ is above $E_C$ the contribution of the electrons to the slowing down is larger than that of the ions. The slowing down time $\tau_s$ is actually the time scale for $V$ to decrease due to electron drag, i.e. $\tau_{\text{slow}}=-V /(d V / d t) e^*$ + Due to the relatively small contribution of fusion power from neutral beams, + only D–T beam–target reactions are included in this model. - $\blacksquare$ + Beam-driven D–D reactions are neglected. - ------------------------- +The beam fusion calculations will only run if the calculated beam current is greater than 0, i.e. when an NBI heating and current drive configuration is active. Currently, the model only runs when the plasma is not ignited `i_plasma_ignited = 0`. -3. **Set the plasma tritium and ion densities** +The NBI parameters taken from the current drive and plasma state are the total beam current (`c_beam_total`), the beam energy (`e_beam_kev`), the tritium fraction in the beam (`f_beam_tritium`), the bulk deuterium and tritium fractions (`f_deuterium_plasma`, `f_tritium_plasma`), and the plasma state quantities required for slowing down. - $$ - \mathtt{deuterium\_density = nd_plasma_fuel_ions_vol_avg * f\_deuterium\_plasma} \\ - \mathtt{tritium\_density = nd_plasma_fuel_ions_vol_avg * f\_tritium\_plasma} - $$ - -4. **Calculate the beam alpha powers, density and deposited energy** - - [`beamcalc()`](#neutral-beam-alpha-power-and-ion-energy--beamcalc) is ran to find the alpha power from the beams, the beam densities and the total energy deposited into the plasma. - -5. **Set the returned alpha power** - - $$ - P_{\alpha,\text{beam}} = \mathtt{beamfus0} \times \left(P_{\alpha,\text{D-beam}} + P_{\alpha,\text{T-beam}}\right) - $$ - -6. **Calculate the neutral beam beta** - - $$ - \beta_{\text{beam}} = \mathtt{betbm0}\times \frac{2}{3}4.03\times10^{-22} \frac{n_{\text{beam}}E_{\text{hot,beam}}}{B_{\text{tot}}^2} - $$ - - The value of $E_{\text{hot,beam}}$ is the energy deposited by the fast beam ions into the plasma, NOT the initial energy of the beam. - ------------------------- - -## Neutral beam alpha powers, density and deposited energy | `beamcalc()` - -1. **Calculate the beam current densities** - - $$ - I_{\text{beam,D}} = I_{\text{beam}} \times \left(1-\mathtt{f\_tritium\_beam}\right) \\ - I_{\text{beam,T}} = I_{\text{beam}} \times \mathtt{f\_tritium\_beam} - $$ +Please see the [H&CD section](../../eng-models/heating_and_current_drive/heating-and-current-drive.md) of the docs for more info. -2. **Calculate the characteristic time taken for the beam energy to comparable to that of the thermal energy** +--- - The attenuation of the beam energy as it penetrates into the plasma is given by[^3]: +## Beam slowing down properties | `beam_fusion()` - $$ - E_{\text{beam}} = E_{\text{beam,0}} \left[e^{-\frac{3t}{\tau_{\text{slow}}}}-\left(\frac{E_{\text{crit}}}{E_{\text{beam,0}}}\right)^{\frac{3}{2}}\left(1-e^{-\frac{3t}{\tau_{\text{slow}}}}\right)\right]^{\frac{2}{3}} - $$ +This section explains the stages of calculation in the function `beam_fusion()`. + +### Calculate the beam ion slowing down time + +The beam slowing down time used in the model is implemented as[^deng1987]: + +$$ +\tau_{\text{slow}} = +1.99\times10^{19} +\left[ +\text{A}_{\text{D}}\left(1-\text{f}_{\text{beam,T}}\right) ++ +\text{A}_{\text{T}}\text{f}_{\text{beam,T}} +\right] +\frac{\langle \text{T}_{\text{e}}\rangle^{3/2}} +{\langle \text{n}_{\text{e}}\rangle \ln \Lambda_{\text{ie}}} +$$ + +where: + +- $A_{\text{D}}$ is the deuteron mass in amu, +- $A_{\text{T}}$ is the triton mass in amu, +- $f_{\text{beam,T}}$ is the tritium fraction in the injected beam, +- $\langle T_{\text{e}}\rangle$ is the density-weighted electron temperature in keV, +- $\langle n_{\text{e}}\rangle$ is the volume-averaged electron density in $\text{m}^{-3}$ +- $\ln\Lambda_{\text{ie}}$ is the ion-electron Coulomb logarithm. +- The numerical prefactor $1.99\times10^{19}$ arises from rewriting the classical slowing-down time expression in practical units. + +--- + +### Calculate the beam critical energy + +For an energetic ion slowing down in a plasma, there is a critical energy $E_{\text{crit}}$ at which the rate of energy loss to ions equals the rate of energy loss to electrons. Above this energy electron drag dominates, while below it ion drag dominates. + +In the current implementation, the deuterium critical energy [^sheffield1994] is + +$$ +\text{E}_{\text{crit,D}} = +14.8 +\text{A}_{\text{D}} +\text{T}_{\text{e}} +\text{Z}_{\text{eff,mw}}^{2/3} +\frac{\ln\Lambda_{\text{ie}}+4.0}{\ln\Lambda_{\text{ie}}} +\qquad [\text{keV}] +$$ + +where $Z_{\text{eff,mw}}$ is the mass-weighted effective plasma charge. + +The tritium critical energy is then scaled by the beam ion mass ratio: + +$$ +\text{E}_{\text{crit,T}} = +\text{E}_{\text{crit,D}} +\left(\frac{\text{A}_{\text{T}}}{\text{A}_{\text{D}}}\right) +$$ + +#### Derivation of beam slowing down rate and critical energy + +The rate of slowing down of a test particle of mass $M$, charge $Ze$ and energy $E$, due to Coulomb collisions with a background species of mass $m_j$, charge $Z_j e$, density $n_j$ and temperature $T_j$, is given by[^deng1987] + +$$ +\frac{\text{dE}}{\text{dt}} += +\left[ +-\Phi(\text{x}_{\text{j}}) ++ +\text{x}_{\text{j}}\left(1+\frac{\text{m}_{\text{j}}}{\text{M}}\,\Phi^{\prime}(\text{x}_{\text{j}})\right) +\right] +\frac{4\pi \text{n}_{\text{j}}}{\text{m}_{\text{j}} \text{V}} +\left(\frac{\text{Z} \text{Z}_{\text{j}} \text{e}^2}{4\pi\varepsilon_0}\right)^2 +\ln \Lambda_{\text{j}} +$$ + +where + +$$ +\Phi^{\prime}(\text{x}) = \frac{\text{d}\Phi}{\text{dx}} +$$ + +$$ +\text{V} = \sqrt{\frac{2\text{E}}{\text{M}}}, \qquad +\text{V}_{\text{j}} = \sqrt{\frac{2\text{kT}_{\text{j}}}{\text{m}_{\text{j}}}}, \qquad +\text{x}_{\text{j}} = \frac{\text{V}}{\text{V}_{\text{j}}} +$$ + +For fast ions in fusion plasmas, the ion contribution and electron contribution may be approximated separately, leading to the standard slowing-down form[^wesson2011]: + +$$ +\frac{\mathrm{d}\text{E}}{\mathrm{d}\text{t}} += +-\frac{\text{A} \text{Z}^2 \sqrt{\text{M}}}{\sqrt{\text{E}}} +- +\frac{\text{B} \text{Z}^2 \text{E}}{\text{M}} +$$ + +with coefficients + +$$ +\begin{aligned} +\text{A} &= +\frac{4\pi}{\sqrt{2}} +\left(\frac{\text{e}^2}{4\pi\varepsilon_0}\right)^2 +\sum_{\text{j}} +\left( +\frac{\text{n}_{\text{j}} \text{Z}_{\text{j}}^2}{\text{m}_{\text{j}}}\ln\Lambda_{\text{j}} +\right) \\ +\text{B} &= +\frac{16\sqrt{\pi}}{3\text{kT}_{\text{e}}} +\sqrt{\frac{\text{m}_{\text{e}}}{2\text{kT}_{\text{e}}}} +\left(\frac{\text{e}^2}{4\pi\varepsilon_0}\right)^2 +\text{n}_{\text{e}} \ln\Lambda_{\text{e}} +\end{aligned} +$$ + +This can be rewritten in the form + +$$ +\frac{\mathrm{d}\text{E}}{\mathrm{d}\text{t}} += +-\frac{2\text{E}}{\tau_{\text{slow}}} +\left[ +1+\left(\frac{\text{E}_{\text{c}}}{\text{E}}\right)^{3/2} +\right] +$$ + +where + +$$ +\text{E}_{\text{c}} += +\left[ +\frac{3\sqrt{\pi}}{4} +\frac{\text{M}^{3/2}}{\text{n}_{\text{e}}\sqrt{\text{m}_{\text{e}}}} +\sum_{\text{j}} +\left( +\frac{\text{n}_{\text{j}} \text{Z}_{\text{j}}^2}{\text{m}_{\text{j}}}\ln\Lambda_{\text{j}} +\right) +\frac{1}{\ln\Lambda_{\text{e}}} +\right]^{2/3} +\text{kT}_{\text{e}} +$$ + +and + +$$ +\tau_{\text{slow}} += +\frac{3(\text{kT}_{\text{e}})^{3/2}}{4\sqrt{2\pi \text{m}_{\text{e}}} \text{Z}^2} +\left(\frac{4\pi\varepsilon_0}{\text{e}^2}\right)^2 +\frac{\text{M}}{\text{n}_{\text{e}} \ln\Lambda_{\text{e}}} +$$ + +In this regime, $\tau_{\text{slow}}$ is the characteristic electron-drag slowing-down timescale. + +$\blacksquare$ + +--- + +#### Set the plasma deuterium and tritium ion densities + +The bulk target ion densities used in the beam-target reactions are + +$$ +\text{n}_{\text{D,plasma}} += +\text{n}_{\text{fuel}} +\text{f}_{\text{D,plasma}} +$$ + +$$ +\text{n}_{\text{T,plasma}} += +\text{n}_{\text{fuel}} +\text{f}_{\text{T,plasma}} +$$ + +where $n_{\text{fuel}}$ is the volume-averaged total fuel ion density. + +--- + +### Calculate the beam alpha powers, beam densities and deposited energy + +[`beam_slowing_down_state()`](#beam-slowing-down-state) is run to calculate the hot beam densities, critical speeds, and deposited beam energy. + +[`beam_reaction_rate_coefficient()`](#beam-reaction-rate-coefficient) is then run for the deuterium and tritium beam components to obtain effective beam-target rate coefficients. + +[`beam_target_reaction_rate()`](#beam-target-reaction-rate) and [`alpha_power_beam()`](#alpha-power-beam) are used to convert these into alpha power. + +### Set the returned alpha power + +The total neutral beam alpha power is + +$$ +\text{P}_{\alpha,\text{beam}} += +\mathtt{beamfus0} +\left( +\text{P}_{\alpha,\text{D-beam}} ++ +\text{P}_{\alpha,\text{T-beam}} +\right) +$$ + +### Calculate the neutral beam beta + +The neutral beam beta is computed from the hot beam density and the deposited beam energy: + +$$ +\beta_{\text{beam}} += +\mathtt{betbm0} +\times +4.03\times10^{-22} +\times +\frac{2}{3} +\frac{\text{n}_{\text{beam,hot}} \text{E}_{\text{beam,deposited}}} +{\text{B}_{\phi}^2 + \text{B}_{\theta}^2} +$$ + +where: + +- $n_{\text{beam,hot}}$ is the total hot beam ion density, +- $E_{\text{beam,deposited}}$ is the density-weighted deposited beam ion energy in keV, +- $B_{\phi}$ is the toroidal magnetic field on axis, +- $B_{\theta}$ is the poloidal magnetic field. + +The value of $E_{\text{beam,deposited}}$ is the pressure-equivalent deposited energy of the hot beam ions, **not** the initial beam injection energy. + +--- + +## Neutral beam alpha power, beam densities and deposited energy | `beam_slowing_down_state()` {#beam-slowing-down-state} + +### Calculate the beam current fractions + +The beam current is split into deuterium and tritium components: + +$$ +\text{I}_{\text{beam,D}} += +\text{I}_{\text{beam}} +\left(1-\text{f}_{\text{beam,T}}\right) +$$ + +$$ +\text{I}_{\text{beam,T}} += +\text{I}_{\text{beam}} +\text{f}_{\text{beam,T}} +$$ - Where the characteristic time taken for the beam energy to fall to that of the thermal energy: +### Calculate the characteristic slowing-down time to the thermal range - $$ - \tau_{\text{slow,D*}} = \frac{\tau_{\text{slow}}}{3}\ln\left(1+\left(\frac{E_{\text{beam,0}}}{E_{\text{crit,D}}}\right)^{\frac{3}{2}}\right) \\ - \tau_{\text{slow,T*}} = \frac{\tau_{\text{slow}}}{3}\ln\left(1+\left(\frac{E_{\text{beam,0}}}{E_{\text{crit,T}}}\right)^{\frac{3}{2}}\right) - $$ +Using the classical slowing-down model, the characteristic time for the beam ion energy to slow from the birth energy to the thermal range is implemented as: + +$$ +\tau_{\text{slow,D}}^{*} += +\frac{\tau_{\text{slow}}}{3} +\ln\left[ +1+ +\left( +\frac{\text{E}_{\text{beam}}}{\text{E}_{\text{crit,D}}} +\right)^{3/2} +\right] +$$ + +$$ +\tau_{\text{slow,T}}^{*} += +\frac{\tau_{\text{slow}}}{3} +\ln\left[ +1+ +\left( +\frac{\text{E}_{\text{beam}}}{\text{E}_{\text{crit,T}}} +\right)^{3/2} +\right] +$$ + +### Set the fast beam ion densities + +The steady-state hot beam ion densities are obtained from the balance between the beam particle source rate and the slowing-down residence time: + +$$ +\text{n}_{\text{beam,D}} += +\frac{\text{I}_{\text{beam,D}}\tau_{\text{slow,D}}^{*}} +{\text{e} \text{V}_{\text{plasma}}} +$$ + +$$ +\text{n}_{\text{beam,T}} += +\frac{\text{I}_{\text{beam,T}}\tau_{\text{slow,T}}^{*}} +{\text{e} \text{V}_{\text{plasma}}} +$$ + +The total hot beam ion density is then + +$$ +\text{n}_{\text{beam,hot}} += +\text{n}_{\text{beam,D}} + \text{n}_{\text{beam,T}} +$$ + +### Calculate the speeds of ions at the critical energy + +Assuming non-relativistic energies, the beam ion speeds at the critical energy are + +$$ +\text{v}_{\text{crit,D}} += +\sqrt{ +\frac{2 \text{e}_{\text{keV}} \text{E}_{\text{crit,D}}} +{\text{m}_{\text{u}} \text{A}_{\text{D}}} +} +$$ + +$$ +\text{v}_{\text{crit,T}} += +\sqrt{ +\frac{2 \text{e}_{\text{keV}} \text{E}_{\text{crit,T}}} +{\text{m}_{\text{u}} \text{A}_{\text{T}}} +} +$$ + +where: + +- $e_{\text{keV}}$ is the conversion from keV to joules, +- $m_u$ is the atomic mass unit. + +### Calculate the fast ion pressures + +First define the source rates per unit volume: + +$$ +\text{S}_{\text{D}} += +\frac{\text{I}_{\text{beam,D}}}{\text{e} \text{V}_{\text{plasma}}} +$$ + +$$ +\text{S}_{\text{T}} += +\frac{\text{I}_{\text{beam,T}}}{\text{e} \text{V}_{\text{plasma}}} +$$ + +The implemented pressure coefficients are then + +$$ +\text{C}_{\text{p},\text{D}} += +\frac{ +\text{A}_{\text{D}} \text{m}_{\text{u}} \tau_{\text{slow}} \text{v}_{\text{crit,D}}^2 \text{S}_{\text{D}} +} +{3 \text{e}_{\text{keV}}} +$$ + +$$ +\text{C}_{\text{p},\text{T}} += +\frac{ +\text{A}_{\text{T}} \text{m}_{\text{u}} \tau_{\text{slow}} \text{v}_{\text{crit,T}}^2 \text{S}_{\text{T}} +} +{3 \text{e}_{\text{keV}}} +$$ + +The fast ion pressures are + +$$ +\text{p}_{\text{D}} += +\text{C}_{\text{p},\text{D}} +\times +\underbrace{ +\left[ +\frac{\text{x}_{\text{c}}^2}{2} ++ +\frac{1}{6}\ln\left(\frac{\text{x}_{\text{c}}^2+2\text{x}_{\text{c}}+1}{\text{x}_{\text{c}}^2-\text{x}_{\text{c}}+1}\right) +- +\frac{1}{\sqrt{3}}\tan^{-1}\left(\frac{2\text{x}_{\text{c}}-1}{\sqrt{3}}\right) +- +\frac{1}{\sqrt{3}}\tan^{-1}\left(\frac{1}{\sqrt{3}}\right) +\right] +}_{\mathtt{fast\_ion\_pressure\_integral}(\text{E}_{\text{beam}},\text{E}_{\text{crit,D}})} +$$ + +$$ +\text{p}_{\text{T}} += +\text{C}_{\text{p},\text{T}} +\times +\underbrace{ +\left[ +\frac{\text{x}_{\text{c}}^2}{2} ++ +\frac{1}{6}\ln\left(\frac{\text{x}_{\text{c}}^2+2\text{x}_{\text{c}}+1}{\text{x}_{\text{c}}^2-\text{x}_{\text{c}}+1}\right) +- +\frac{1}{\sqrt{3}}\tan^{-1}\left(\frac{2\text{x}_{\text{c}}-1}{\sqrt{3}}\right) +- +\frac{1}{\sqrt{3}}\tan^{-1}\left(\frac{1}{\sqrt{3}}\right) +\right] +}_{\mathtt{fast\_ion\_pressure\_integral}(\text{E}_{\text{beam}},\text{E}_{\text{crit,T}})} +$$ + +with + +$$ +\text{x}_{\text{c}} = \sqrt{\frac{\text{E}_{\text{beam}}}{\text{E}_{\text{crit}}}} +$$ + +--- + +### Beam fast ion pressure integral | `fast_ion_pressure_integral()` + +This internal function returns the dimensionless pressure integral factor used in the fast-ion pressure expression. It takes the beam birth energy $E_{\text{beam}}$ and the critical energy $E_{\text{crit}}$ for the required beam ion species. + +#### Derivation + +The fast ions, because of their finite slowing down time, develop a finite pressure which must be supported by the magnetic field. + +To calculate this pressure, introduce the distribution function $g(E)$ of fast ions, where $g(E)$ is the number of ions per unit energy per unit volume, satisfying the steady-state kinetic equation + +$$ +\frac{\partial}{\partial \text{E}} +\left( +\text{g} \frac{\mathrm{d}\text{E}}{\mathrm{d}\text{t}} +\right) += +\text{S}(\text{E}) +$$ + +with slowing-down law + +$$ +\frac{\mathrm{d}\text{E}}{\mathrm{d}\text{t}} += +-\frac{2\text{E}}{\tau_{\text{s}}} +\left[ +1+\left(\frac{\text{E}_{\text{c}}}{\text{E}}\right)^{3/2} +\right] +$$ + +If the ions are born monoenergetically, + +$$ +\text{S}(\text{E})=\text{S}_0\delta(\text{E}-\text{E}_0) +$$ + +and with the boundary condition $g(E)=0$ for $E>E_0$, the steady-state distribution becomes + +$$ +\text{g}(\text{E})= +\begin{cases} +\dfrac{\text{S}_0\tau_{\text{s}}}{2\text{E}\left[1+\left(\text{E}_{\text{c}}/\text{E}\right)^{3/2}\right]}, & \text{E}<\text{E}_0 \\ +0, & \text{E}>\text{E}_0 +\end{cases} +$$ + +The fast ion pressure is then + +$$ +\text{p} += +\frac{2}{3} +\int_0^{\text{E}_0} \text{g}(\text{E}) \text{E} \, \mathrm{d}\text{E} +$$ + +which can be written as + +$$ +\text{p} += +\frac{\text{M} \text{S}_0 \tau_{\text{s}} \text{V}_{\text{c}}^2}{3} +\int_0^{\text{x}_{\text{c}}} +\frac{\text{x}^4}{1+\text{x}^3}\,\mathrm{d}\text{x} +$$ + +and evaluated analytically as + +$$ +\text{p} += +\frac{\text{M} \text{S}_0 \tau_{\text{s}} \text{V}_{\text{c}}^2}{3} +\left[ +\frac{\text{x}_{\text{c}}^2}{2} ++ +\frac{1}{6}\ln\left(\frac{\text{x}_{\text{c}}^2+2\text{x}_{\text{c}}+1}{\text{x}_{\text{c}}^2-\text{x}_{\text{c}}+1}\right) +- +\frac{1}{\sqrt{3}}\tan^{-1}\left(\frac{2\text{x}_{\text{c}}-1}{\sqrt{3}}\right) +- +\frac{1}{\sqrt{3}}\tan^{-1}\left(\frac{1}{\sqrt{3}}\right) +\right] +$$ - This is calculated for the deuterium and the tritium beam components +`fast_ion_pressure_integral()` returns the bracketed dimensionless factor only. -3. **Set the fast beam ion densities** +$\blacksquare$ - $$ - \langle n_{\text{beam}} \rangle_{\text{D}} = \frac{I_{\text{beam,D}}\tau_{\text{slow,D*}}}{V_{\text{plasma}}} \\ - \langle n_{\text{beam}} \rangle_{\text{T}} = \frac{I_{\text{beam,T}}\tau_{\text{slow,T*}}}{V_{\text{plasma}}} - $$ +--- - We can also set the total hot ion beam density as: +### Calculate the deposited fast ion energy from the pressure - $$ - \langle n_{\text{beam}} \rangle_{\text{total}} = \langle n_{\text{beam}} \rangle_{\text{D}} + \langle n_{\text{beam}} \rangle_{\text{T}} - $$ +The deposited beam ion energy is inferred from the pressure relation -4. **Calculate the speeds of ions when at the critical energy** +$$ +\text{p} = \frac{1}{3}\text{nmv}_{\text{rms}}^2 = \frac{2}{3}\text{n}\langle \text{E}\rangle +$$ - Assuming non-relativistic energies, we set the velocities of the deuterium and tritium particles when they have the critical energy: +so that the pressure-equivalent deposited energies are - $$ - v_{\text{crit,D}} = \sqrt{\left(2m_{\text{D}}E_{\text{crit,D}}\right)} \quad - v_{\text{crit,T}} = \sqrt{\left(2m_{\text{T}}E_{\text{crit,T}}\right)} - $$ +$$ +\text{E}_{\text{hot,D}} += +\frac{3}{2} +\frac{\text{p}_{\text{D}}}{\text{n}_{\text{beam,D}}} +$$ -5. **Calculate the fast ion pressures** - - The fast ion pressure is set as: +$$ +\text{E}_{\text{hot,T}} += +\frac{3}{2} +\frac{\text{p}_{\text{T}}}{\text{n}_{\text{beam,T}}} +$$ - $$ - p_{\text{D}}=\frac{m_{\text{D}} \tau_{\text{slow}} v_{\text{crit,D}}^2 I_{\text{beam,D}}}{3 V_{\text{plasma}}} \times \\ - \underbrace{\left[\frac{x_c^2}{2}+\frac{1}{6} \ln \left(\frac{x_c^2+2 x_c+1}{x_c^2-x_c+1}\right)-\frac{1}{\sqrt{3}} \tan ^{-1}\left(\frac{2 x_c-1}{\sqrt{3}}\right)-\frac{1}{\sqrt{3}} \tan \left(\frac{1}{\sqrt{3}}\right)\right]}_{\mathtt{\_fast\_ion\_pressure\_integral(E_{\text{beam}},E_{\text{crit,D}})}} \\ - p_{\text{T}}=\frac{m_{\text{T}} \tau_{\text{slow}} v_{\text{crit,T}}^2 I_{\text{beam,T}}}{3 V_{\text{plasma}}} \times \\ \quad - \underbrace{\left[\frac{x_c^2}{2}+\frac{1}{6} \ln \left(\frac{x_c^2+2 x_c+1}{x_c^2-x_c+1}\right)-\frac{1}{\sqrt{3}} \tan ^{-1}\left(\frac{2 x_c-1}{\sqrt{3}}\right)-\frac{1}{\sqrt{3}} \tan \left(\frac{1}{\sqrt{3}}\right)\right]}_{\mathtt{\_fast\_ion\_pressure\_integral(E_{\text{beam}},E_{\text{crit,T}})}} - $$ +with the convention that the deposited energy is set to zero if the corresponding hot beam density is zero. - --------------------------------- +The total deposited beam ion energy returned to `beam_fusion()` is the density-weighted average beam ion energy: - ### Beam fast ion pressure integral | `fast_ion_pressure_integral()` +$$ +\langle \text{E}_{\text{hot}} \rangle_{\text{beam}} += +\frac{ +\text{E}_{\text{hot,D}} \, \text{n}_{\text{beam,D}} ++ +\text{E}_{\text{hot,T}} \, \text{n}_{\text{beam,T}} +} +{\text{n}_{\text{beam,hot}}} +$$ - This internal function is for returning the main integral value for the fast ion pressure. It takes the initial beam energy ($E_{\text{beam}}$) and the critical energy ($E_{\text{crit}}$) for the required ion species. +$$ +\text{E}_{\text{hot,total}} = \langle \text{E}_{\text{hot}} \rangle_{\text{beam}} +$$ - This integral derivation is originally of the form derived by D.Baiquan et.al.[^2] +with the convention that this is set to zero if $\text{n}_{\text{beam,hot}} = 0$. This represents a density-weighted average over the deuterium and tritium fast-ion populations. - #### Derivation +### Return the slowing-down state - The fast ions, because of their finite slowing down time, develop a certain amount of pressure which has to be supported by the magnetic field. +`beam_slowing_down_state()` returns the key quantities describing the steady-state fast-ion population: - This is in addition to the pressure of accumulated thermal "ash" in the plasma. To calculate this pressure we introduce a kinetic equation for the slowing down particles and solve it for their distribution function. Integration of the distribution function then determines the fast ion pressure. +- deuterium beam ion density, +- tritium beam ion density, +- deuterium critical speed, +- tritium critical speed, +- total hot beam ion density, +- density-weighted deposited beam energy. - We introduce the distribution function $g$ , of fast ions; $g$ is defined as the number of ions per unit energy per unit spatial volume and satisfies the steady-state kinetic equation, +--- - $$ - \frac{\partial}{\partial E}\left(g \frac{d E}{d t}\right)=S(E), - $$ +## Beam fusion reaction rate coefficient | `beam_reaction_rate_coefficient()` {#beam-reaction-rate-coefficient} - where $S(E)$ is the source function in this "phase" space and $d E / d t$ is given by: +### Calculate the beam velocity - $$ - \frac{dE}{dT} = -\frac{2E}{\tau_s}\left[1+\left(\frac{E_c}{E}^{\frac{3}{2}}\right)\right] - $$ +The beam velocity is calculated from the input beam energy and beam ion mass: - This is the same form as above for deriving the critical energy and beam slow time in [`beam_fusion()`](#beam-slowing-down-properties--beam_fusion) +$$ +\text{E}_{\text{beam}}[\mathrm{J}] = +\text{E}_{\text{beam}}[\mathrm{keV}] +\left(1.602176634\times10^{-16}\,\mathrm{J\,keV^{-1}}\right) +$$ - This equation assumes the ions have a confinement time much longer than the slowing down time. For cases in which this is not true, an additional loss term would have to be introduced in the equation above. If we assume the ions are born monoenergetically, then: +$$ +\text{v}_{\text{beam}} = +\sqrt{ +\frac{ +2\,\text{E}_{\text{beam}}[\mathrm{J}] +}{ +\text{A}\,\text{m}_{\text{u}} +} +} +$$ - $$ - S(E)=S_0 \delta\left(E-E_0\right), - $$ +### Define the integral coefficient - where $\mathrm{S}_0$ is the number of ions born per unit time per unit volume and $\delta$ is the [Dirac delta function](https://en.wikipedia.org/wiki/Dirac_delta_function). +The implemented coefficient is - We impose the boundary condition that $\mathrm{g}=0$ for $\mathrm{E}>\mathrm{E}_0$. Our steady state kinetic distribution function can then be integrated to yield: +$$ +\frac{3\text{v}_{\text{crit}}} +{\ln\left(1+\left(\frac{\text{v}_{\text{beam}}}{\text{v}_{\text{crit}}}\right)^3\right)} +$$ - $$ - g(E)= \begin{cases}\frac{S_0 T_s}{2 E\left(1+\left(E_c / E\right)^{3 / 2}\right)}, & EE_0\end{cases} - $$ +### Perform the fusion rate integral - If we consider the appropriate limiting cases. +The slowing-down-weighted integral is evaluated as - The pressure $p$ of the fast ions is +$$ +\int_0^{\text{v}_{\text{beam}}/\text{v}_{\text{crit}}} +\frac{\text{u}^3}{1+\text{u}^3} +\sigma_{\text{bmfus}}(\text{E}_{\text{amu}}(\text{u})) +\,\mathrm{d}\text{u} +$$ - $$ - p=\frac{2}{3} \int_0^{E_0} g(E) E \ \ \mathrm{dE} - $$ +where $u = v/v_{\text{crit}}$. - which can be rewritten as +The quantity returned by `_beam_fusion_cross_section()` is in cm$^2$, so the integrated area is converted to m$^2$ before forming the final rate coefficient. - $$ - p=\frac{M S_0 \tau_s V_c^2}{3} \int_0^{X_c} \frac{x^4}{1+x^3} \quad \mathrm{dx} - $$ +--- - The integral above can also be evaluated analytically to yield: +### Hot beam fusion reaction rate integrand | `_hot_beam_fusion_reaction_rate_integrand()` - $$ - p=\frac{M \tau_s V_c^2 S_0}{3}\left[\frac{x_c^2}{2}+\frac{1}{6} \ln \left(\frac{x_c^2+2 x_c+1}{x_c^2-x_c+1}\right) \\ - -\frac{1}{\sqrt{3}} \tan ^{-1}\left(\frac{2 x_c-1}{\sqrt{3}}\right)-\frac{1}{\sqrt{3}} \tan \left(\frac{1}{\sqrt{3}}\right)\right] - $$ +This internal function evaluates the integrand used in the beam-target fusion rate coefficient. - For most applications to fusion born particles, the dominant term is the first term in the square brackets. - This function returns the terms in the square brackets. +The implemented integrand is - $\blacksquare$ +$$ +\frac{\text{u}^3}{1+\text{u}^3}\sigma_{\text{bmfus}}(\text{E}_{\text{amu}}(\text{u})) +$$ - -------------------------- +where: -6. **Calculate the deposited fast ion energy from the pressure** +- $u$ is the ratio of the instantaneous beam speed to the critical speed, +- $E_{\text{amu}}(u)$ is the instantaneous beam kinetic energy per amu, +- $\sigma_{\text{bmfus}}$ is returned by [`_beam_fusion_cross_section()`](#beam-fusion-cross-section). - This can be simply done by applying the Ideal gas law approximation that, $P = \frac{1}{3}nmv_{\text{rms}^2} = \frac{2}{3}n \langle E \rangle$ +The beam kinetic energy per amu used in the code is constructed from - $$ - E_{\text{hot,D}} = \frac{3}{2}p_{\text{D}} \langle n_{\text{beam}} \rangle_{\text{D}} \\ - E_{\text{hot,T}} = \frac{3}{2}p_{\text{T}} \langle n_{\text{beam}} \rangle_{\text{T}} - $$ +$$ +\text{E}_{\text{amu}} += +\frac{(\text{u} \text{v}_{\text{crit}})^2 \text{m}_{\text{u}}}{\text{e}_{\text{keV}}} +$$ - We can thus also define the total hot ion depoisted energy as: +#### Beam fusion cross section | `_beam_fusion_cross_section()` {#beam-fusion-cross-section} - $$ - E_{\text{hot,total}} = \frac{\left( E_{\text{hot,D}} \langle n_{\text{beam}} \rangle_{\text{D}}\right) + \left( E_{\text{hot,T}} \langle n_{\text{beam}} \rangle_{\text{T}}\right)}{\langle n_{\text{beam}} \rangle_{\text{total}}} - $$ +This internal function returns the beam fusion cross-section fit used by the beam-target rate coefficient calculation. Note: the provenance of this fit is currently unverified in the documentation. -7. **Calculate the hot ion species fusion rates** +The plasma ions are assumed to be stationary. - The D-T fusion reaction rate is calculated from the [`beam_reaction_rate()`](#beam-fusion-reaction-rate--beam_reaction_rate) function. +The implementation is - ------------------------ +$$ +\sigma_{\text{bm}}(\text{E}) = +\begin{cases} +1.0\times10^{-27}\ \text{cm}^2, & \text{E} < 10.0\ \text{keV} \\ +8.0\times10^{-26}\ \text{cm}^2, & \text{E} > 10^4\ \text{keV} \\ +\dfrac{ +1.0\times10^{-24} +\left[ +\dfrac{\text{a}_2}{1.0+(\text{a}_3 \text{E}-\text{a}_4)^2}+\text{a}_5 +\right] +}{ +\text{E}\left[\exp\left(\dfrac{\text{a}_1}{\sqrt{\text{E}}}\right)-1.0\right] +} +\ \text{cm}^2, & \text{otherwise} +\end{cases} +$$ - ### Beam fusion reaction rate | `beam_reaction_rate()` +where the beam energy used inside the fit is - 1. **Calculate beam velocity** +$$ +\text{E} += +\frac{1}{2} \text{A}_{\text{D}} \text{E}_{\text{amu}} +$$ - The beam velocity ($v_{\text{beam}}$) is calculated from the inputted beam energy and relative ion mass of the beam by simply re-arranging the kinetic energy equation +and the constants are - 2. **Define the integral coefficient** +- $a_1 = 45.95$ +- $a_2 = 5.02 \times 10^4$ +- $a_3 = 1.368 \times 10^{-2}$ +- $a_4 = 1.076$ +- $a_5 = 4.09 \times 10^2$ - $$ - \frac{3v_{\text{critical}}}{\ln\left(1+\frac{v_{\text{beam}}}{v_{\text{critical}}}\right)^{3}} - $$ +The exact provenance of this particular fit and its coefficients has not yet been fully traced in the present documentation. It is retained here to document the current implementation. - 3. **Perform the fusion rate integral** +--- - $$ - \int_0^{v_{\text{relative}}} \frac{u^3}{1+u^3}\sigma_{\text{bmfus}}(E_{\text{amu}}) - $$ - - -------------------- +### Multiply by the coefficient to get the full rate coefficient - #### Hot Beam Fusion Reaction Rate Integrand | `_hot_beam_fusion_reaction_rate_integrand()` +The final effective beam-target fusion rate coefficient is - This function computes the integrand for the hot beam fusion reaction rate based on the ratio of beam velocity to the critical velocity and the critical velocity for electron/ion slowing down of the beam ion. +$$ +\langle \sigma \text{v} \rangle_{\text{beam}} += +\frac{3\text{v}_{\text{crit}}} +{\ln\left(1+\left(\frac{\text{v}_{\text{beam}}}{\text{v}_{\text{crit}}}\right)^3\right)} +\int_0^{\text{v}_{\text{beam}}/\text{v}_{\text{crit}}} +\frac{\text{u}^3}{1+\text{u}^3} +\sigma_{\text{bmfus}}(\text{E}_{\text{amu}}(\text{u})) +\,\mathrm{d}\text{u} +$$ - The integrand function is: +--- - $$ - \int \frac{u^3}{1+u^3}\sigma_{\text{bmfus}}(E_{\text{amu}}) - $$ +## Beam-target fusion reaction rate | `beam_target_reaction_rate()` {#beam-target-reaction-rate} - Where $u$ is the inputted ratio of the beam to the critical velocity. - $E_{\text{amu}}$ represents the beam kinetic energy per atomic mass unit. +`beam_target_reaction_rate()` calculates the total beam-target fusion reaction rate from the hot beam ion density, thermal target ion density, effective beam-target rate coefficient, and plasma volume. This is consistent with simple beam-plasma reaction-rate models in which fast ions interact with a Maxwellian background plasma during the slowing-down process[^niikura1990]. - The calculated beam fusion cross section $\sigma_{\text{bmfus}}$ is calculated from [`_beam_fusion_cross_section()`](#beam-fusion-cross-section--_beam_fusion_cross_section) +$$ +\text{R}_{\text{beam-target}} += +\text{n}_{\text{beam}} +\text{n}_{\text{target}} +\langle \sigma \text{v} \rangle_{\text{beam}} +\text{V}_{\text{plasma}} +$$ - ------------------------ +where: - #### Beam fusion cross section | `_beam_fusion_cross_section()` +- $\text{n}_{\text{beam}}$ is the hot beam ion density, +- $\text{n}_{\text{target}}$ is the thermal target ion density, +- $\langle \sigma v \rangle_{\text{beam}}$ is the effective beam-target rate coefficient, +- $\text{V}_{\text{plasma}}$ is the plasma volume. - This internal function is used to find the beam cross section. - It sets limits on cross-section at low and high beam energies. The plasma ions are assumed to be stationary: +This returns the total reaction rate in $\text{s}^{-1}$. - $$ - \sigma_{\text{bm}}(E) = - \begin{cases} - 1.0 \times 10^{-27} \ \text{cm}^2 & \text{if } E < 10.0 \ \text{keV/amu} \\ - 8.0 \times 10^{-26} \ \text{cm}^2 & \text{if } E > 10^4 \ \text{keV/amu} \\ - \frac{1.0 \times 10^{-24} \cdot \left( \frac{a_2}{1.0 + (a_3 E - a_4)^2} + a_5 \right)}{E \left( \exp\left(\frac{a_1}{\sqrt{E}}\right) - 1.0 \right)} \ \text{cm}^2 & \text{otherwise} - \end{cases} - $$ +--- - where: +## Beam fusion alpha power | `alpha_power_beam()` {#alpha-power-beam} - - \( E \) is the beam energy in, $\text{keV/amu}$ +The beam-target alpha power is obtained from the total beam-target reaction rate by - - \( a_1, a_2, a_3, a_4, a_5 \) are constants. +$$ +\text{P}_{\alpha,\text{beam}} += +\text{R}_{\text{beam-target}} \text{E}_{\alpha} +$$ - The constants are defined as: +and is converted from W to MW in the implementation. - - \( a_1 = 45.95 \) - - \( a_2 = 5.02 \times 10^4 \) - - \( a_3 = 1.368 \times 10^{-2} \) - - \( a_4 = 1.076 \) - - \( a_5 = 4.09 \times 10^2 \) +This returns the alpha power in MW. - ---------------------- +--- - 4. **Multiply by the coefficient to get the full fusion rate** +## Full beam-target alpha power calculation in `beam_fusion()` - $$ - \frac{3v_{\text{critical}}}{\ln\left(1+\frac{v_{\text{beam}}}{v_{\text{critical}}}\right)^{3}}\int_0^{v_{\text{relative}}} \frac{u^3}{1+u^3}\sigma_{\text{bmfus}}(E_{\text{amu}}) - $$ +For the deuterium beam component reacting with thermal tritium: + +$$ +\text{R}_{\text{D-beam,DT}} += +\text{n}_{\text{beam,D}} +\text{n}_{\text{T,plasma}} +\langle \sigma \text{v} \rangle_{\text{beam,D}} +\text{V}_{\text{plasma}} +$$ - ------------------------- +$$ +\text{P}_{\alpha,\text{D-beam}} += +\text{R}_{\text{D-beam,DT}} \text{E}_{\alpha} +$$ -8. **Calculate the alpha power produced by the hot ion species** +For the tritium beam component reacting with thermal deuterium: - The function [`alpha_power_beam()`](#beam-fusion-alpha-power--alpha_power_beam) is ran to calculate the alpha power produced by the deuterium and tritium fast ions. +$$ +\text{R}_{\text{T-beam,DT}} += +\text{n}_{\text{beam,T}} +\text{n}_{\text{D,plasma}} +\langle \sigma \text{v} \rangle_{\text{beam,T}} +\text{V}_{\text{plasma}} +$$ - ----------------------- +$$ +\text{P}_{\alpha,\text{T-beam}} += +\text{R}_{\text{T-beam,DT}} \text{E}_{\alpha} +$$ - ### Beam fusion alpha power | `alpha_power_beam()` +The returned beam alpha power is then - 1. **Calculate reactivity ratio** +$$ +\text{P}_{\alpha,\text{beam}} += +\mathtt{beamfus0} +\left( +\text{P}_{\alpha,\text{D-beam}} ++ +\text{P}_{\alpha,\text{T-beam}} +\right) +$$ - The ratio between the profile averaged reactivity for D-T reactions and the reactivity for the D-T reactions if the plasma is assumed to be homogeneously at the volume averaged ion temperature ($T_{\text{i}}$) is calculated. +--- - $$ - f_{\text{DT}} = \frac{\langle\langle \sigma v \rangle\rangle_{\text{DT}}}{\langle \sigma v \rangle_{\text{DT}}} - $$ +## Key Constraints - 2. **Calculate the alpha fusion power** +### Hot beam ion density limit - The alpha powers from the deuterium and tritium beam components are calculated: +This constraint can be activated by stating `icc = 7` in the input file. - $$ - P_{\alpha,\text{D}} = \langle n_{\text{beam}} \rangle_{\text{D}} n_{\text{i}} \langle \sigma v \rangle_{\text{beam}} E_{\alpha} V_{\text{plasma}} f_{\text{DT}} \\ - P_{\alpha,\text{T}} = \langle n_{\text{beam}} \rangle_{\text{T}} n_{\text{i}} \langle \sigma v \rangle_{\text{beam}} E_{\alpha} V_{\text{plasma}} f_{\text{DT}} - $$ +The desired value of the hot ion beam density calculated from the code (`nd_beam_ions_out`) can be constrained using the input variable `f_nd_beam_electron`, which is the ratio of the beam density to the plasma electron density. It can be set as an iteration variable by setting `ixc = 7`. ------------------------------- +--- -## Key Constraints +## References -### Hot beam ion density limit +[^sheffield1994]: J. W. Sheffield, “The physics of magnetic fusion reactors,” *Rev. Mod. Phys.*, vol. 66, no. 3, pp. 1015–1103, Jul. 1994. -This constraint can be activated by stating `icc = 7` in the input file. +[^deng1987]: B. Deng and G. A. Emmert, “Fast ion pressure in fusion plasma,” *Nuclear Fusion and Plasma Physics*, vol. 9, no. 3, pp. 136–141, 1987. Available: -The desired value of the hot ion beam density calculated from the code (`nd_beam_ions_out`) can be constrained using the input variable, `f_nd_beam_electron`. Which is the ratio of the beam density to the plasma electron density. It can be set as an iteration variable by setting `ixc = 7`. +[^wesson2011]: J. Wesson, *Tokamaks*, 4th ed., Oxford Science Publications, 2011. -[^1]: J. W. Sheffield, “The physics of magnetic fusion reactors,” vol. 66, no. 3, pp. 1015–1103,Jul. 1994, doi: https://doi.org/10.1103/revmodphys.66.1015. -[^2]: Deng Baiquan and G. A. Emmert, “Fast ion pressure in fusion plasma,” Nuclear Fusion and Plasma Physics,vol. 9, no. 3, pp. 136–141, 2022, Available: https://fti.neep.wisc.edu/fti.neep.wisc.edu/pdf/fdm718.pdf -[^3]: Wesson, J. (2011) Tokamaks. 4th Edition, 2011 Oxford Science Publications,International Series of Monographs on Physics, Volume 149. \ No newline at end of file +[^niikura1990]: S. Niikura and M. Nagami, “Improvement of fusion reactivity and fusion power multiplication factor in the presence of fast ions,” *Fusion Engineering and Design*, vol. 12, pp. 467–480, 1990. diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py index 4fcf8c46f8..972f92d297 100644 --- a/process/models/physics/fusion_reactions.py +++ b/process/models/physics/fusion_reactions.py @@ -880,8 +880,7 @@ def set_fusion_powers( def beam_fusion( beamfus0: float, betbm0: float, - b_plasma_poloidal_average: float, - b_plasma_toroidal_on_axis: float, + b_plasma_total: float, c_beam_total: float, nd_plasma_electrons_vol_avg: float, nd_plasma_fuel_ions_vol_avg: float, @@ -890,69 +889,63 @@ def beam_fusion( f_deuterium_plasma: float, f_tritium_plasma: float, f_beam_tritium: float, - sigmav_dt_average: float, temp_plasma_electron_density_weighted_kev: float, - temp_plasma_ion_density_weighted_kev: float, vol_plasma: float, n_charge_plasma_effective_mass_weighted_vol_avg: float, -) -> tuple: - """Routine to calculate beam slowing down properties. +) -> tuple[float, float, float]: + """Calculate neutral beam slowing-down properties and beam-target fusion. - This function computes the neutral beam beta component, hot beam ion density, - and alpha power from hot neutral beam ions based on the provided plasma parameters. + This function computes the neutral beam beta contribution, hot beam ion density, + and alpha power generated from beam-target fusion reactions in the plasma. Parameters ---------- beamfus0: - Multiplier for beam-background fusion calculation. + Multiplier for beam target fusion calculation. betbm0: Leading coefficient for neutral beam beta fraction. - b_plasma_poloidal_average: - Poloidal field (T). - b_plasma_toroidal_on_axis: - Toroidal field on axis (T). - c_beam_total: - Neutral beam current (A). - nd_plasma_electrons_vol_avg: - Electron density (m^-3). - nd_plasma_fuel_ions_vol_avg: - Fuel ion density (m^-3). - ion_electron_coulomb_log: - Ion-electron coulomb logarithm. - e_beam_kev: - Neutral beam energy (keV). - f_deuterium_plasma: - Deuterium fraction of main plasma. - f_tritium_plasma: - Tritium fraction of main plasma. - f_beam_tritium: - Tritium fraction of neutral beam. - sigmav_dt_average: - Profile averaged for D-T (m^3/s). - temp_plasma_electron_density_weighted_kev: - Density-weighted electron temperature (keV). - temp_plasma_ion_density_weighted_kev: - Density-weighted ion temperature (keV). - vol_plasma: + b_plasma_poloidal_average : + Poloidal magnetic field (T). + b_plasma_toroidal_on_axis : + Toroidal magnetic field on axis (T). + c_beam_total : + Total neutral beam current (A). + nd_plasma_electrons_vol_avg : + Volume-averaged electron density (m^-3). + nd_plasma_fuel_ions_vol_avg : + Volume-averaged fuel ion density (m^-3). + ion_electron_coulomb_log : + Ion-electron Coulomb logarithm. + e_beam_kev : + Neutral beam injection energy (keV). + f_deuterium_plasma : + Deuterium fraction in the bulk plasma. + f_tritium_plasma : + Tritium fraction in the bulk plasma. + f_beam_tritium : + Tritium fraction in the injected beam. + temp_plasma_electron_density_weighted_kev : + Density weighted average electron temperature (keV). + vol_plasma : Plasma volume (m^3). - n_charge_plasma_effective_mass_weighted_vol_avg: - Mass weighted plasma effective charge. + n_charge_plasma_effective_mass_weighted_vol_avg : + Mass-weighted effective plasma charge. Returns ------- - : - tuple: A tuple containing the following elements: - - beta_beam (float): Neutral beam beta component. - - nd_beam_ions_out (float): Hot beam ion density (m^-3). - - p_beam_alpha_mw (float): Alpha power from hot neutral beam ions (MW). + beta_beam : + Neutral beam beta contribution. + n_beam : + Hot beam ion density (m^-3). + p_alpha_beam : + Alpha power from beam-target fusion (MW). Notes ----- - - The function uses the Bosch-Hale parametrization to compute the reactivity. - - The critical energy for electron/ion slowing down of the beam ion is calculated - for both deuterium and tritium neutral beams. - - The function integrates the hot beam fusion reaction rate integrand over the - range of beam velocities up to the critical velocity. + - The function uses the Bosch-Hale parametrization to evaluate fusion reactivity. + - Beam ion slowing-down properties are computed using Coulomb collision physics. + - Fusion power is obtained by integrating the beam-target reaction rate over + the fast-ion distribution. References ---------- @@ -960,15 +953,16 @@ def beam_fusion( Nuclear Fusion, vol. 32, no. 4, pp. 611-631, Apr. 1992, doi: https://doi.org/10.1088/0029-5515/32/4/i07. - - J. W. Sheffield, “The physics of magnetic fusion reactors,” vol. 66, no. 3, pp. 1015-1103, - Jul. 1994, doi: https://doi.org/10.1103/revmodphys.66.1015. + - J. W. Sheffield, “The physics of magnetic fusion reactors,” + Rev. Mod. Phys., vol. 66, no. 3, pp. 1015-1103, Jul. 1994, Jul. 1994, + doi: https://doi.org/10.1103/revmodphys.66.1015. - - Deng Baiquan and G. A. Emmert, “Fast ion pressure in fusion plasma,” Nuclear Fusion and Plasma Physics, - vol. 9, no. 3, pp. 136-141, 2022, Available: https://fti.neep.wisc.edu/fti.neep.wisc.edu/pdf/fdm718.pdf - ‌ + - Deng Baiquan and G. A. Emmert, “Fast ion pressure in fusion plasma,” + Nuclear Fusion and Plasma Physics, vol. 9, no. 3, pp. 136-141, 2022. + Available: https://fti.neep.wisc.edu/fti.neep.wisc.edu/pdf/fdm718.pdf """ # Beam ion slowing down time given by Deng Baiquan and G. A. Emmert 1987 - beam_slow_time = ( + t_beam_slow = ( 1.99e19 * ( constants.M_DEUTERON_AMU * (1.0 - f_beam_tritium) @@ -994,100 +988,142 @@ def beam_fusion( ) # Deuterium and tritium ion densities - deuterium_density = nd_plasma_fuel_ions_vol_avg * f_deuterium_plasma - tritium_density = nd_plasma_fuel_ions_vol_avg * f_tritium_plasma - - ( - deuterium_beam_alpha_power, - tritium_beam_alpha_power, - hot_beam_density, - beam_deposited_energy, - ) = beamcalc( - deuterium_density, - tritium_density, + nd_plasma_deuterium = nd_plasma_fuel_ions_vol_avg * f_deuterium_plasma + nd_plasma_tritium = nd_plasma_fuel_ions_vol_avg * f_tritium_plasma + + beam_state = beam_slowing_down_state( e_beam_kev, critical_energy_deuterium, critical_energy_tritium, - beam_slow_time, + t_beam_slow, f_beam_tritium, c_beam_total, - temp_plasma_ion_density_weighted_kev, vol_plasma, - sigmav_dt_average, + ) + + sigv_beam_deuterium = beam_reaction_rate_coefficient( + constants.M_DEUTERON_AMU, + beam_state.deuterium_critical_energy_speed, + e_beam_kev, + ) + + sigv_beam_tritium = beam_reaction_rate_coefficient( + constants.M_TRITON_AMU, + beam_state.tritium_critical_energy_speed, + e_beam_kev, + ) + + reaction_rate_beam_deuterium_dt = beam_target_reaction_rate( + beam_state.deuterium_beam_density, + nd_plasma_tritium, + sigv_beam_deuterium, + vol_plasma, + ) + p_alpha_beam_deuterium_dt_mw = alpha_power_beam( + reaction_rate_beam_deuterium_dt, + ) + + reaction_rate_beam_tritium_dt = beam_target_reaction_rate( + beam_state.tritium_beam_density, + nd_plasma_deuterium, + sigv_beam_tritium, + vol_plasma, + ) + p_alpha_beam_tritium_dt_mw = alpha_power_beam( + reaction_rate_beam_tritium_dt, ) # Neutral beam alpha power - p_beam_alpha_mw = beamfus0 * (deuterium_beam_alpha_power + tritium_beam_alpha_power) + p_beam_alpha_mw = beamfus0 * ( + p_alpha_beam_deuterium_dt_mw + p_alpha_beam_tritium_dt_mw + ) # Neutral beam beta beta_beam = ( betbm0 * 4.03e-22 * (2 / 3) - * hot_beam_density - * beam_deposited_energy - / (b_plasma_toroidal_on_axis**2 + b_plasma_poloidal_average**2) + * beam_state.nd_beam_hot + * beam_state.e_beam_deposited_kev + / (b_plasma_total**2) ) - return beta_beam, hot_beam_density, p_beam_alpha_mw + return beta_beam, beam_state.nd_beam_hot, p_beam_alpha_mw + + +@dataclass(frozen=True) +class BeamSlowingDownState: + """Beam slowing-down state and hot-ion properties. + + Attributes + ---------- + deuterium_beam_density : + Steady-state deuterium beam ion density (m^-3). + tritium_beam_density : + Steady-state tritium beam ion density (m^-3). + deuterium_critical_energy_speed : + Critical speed of deuterium beam ions corresponding to the + critical energy (m/s). + tritium_critical_energy_speed : + Critical speed of tritium beam ions corresponding to the + critical energy (m/s). + nd_beam_hot : + Total hot beam ion density (m^-3). + e_beam_deposited_kev : + Density-weighted average deposited beam ion energy (keV). + """ + + deuterium_beam_density: float + tritium_beam_density: float + deuterium_critical_energy_speed: float + tritium_critical_energy_speed: float + nd_beam_hot: float + e_beam_deposited_kev: float -def beamcalc( - nd: float, - nt: float, +def beam_slowing_down_state( e_beam_kev: float, critical_energy_deuterium: float, critical_energy_tritium: float, - beam_slow_time: float, + t_beam_slow: float, f_beam_tritium: float, c_beam_total: float, - temp_plasma_ion_vol_avg_kev: float, vol_plasma: float, - svdt: float, -) -> tuple[float, float, float, float]: - """Calculate neutral beam alpha power and ion energy. +) -> BeamSlowingDownState: + """Calculate beam slowing-down state and hot ion properties. - This function computes the alpha power generated from the interaction between - hot beam ions and thermal ions in the plasma, as well as the hot beam ion density - and average hot beam ion energy. + Estimates the steady-state hot D and T beam-ion densities from beam + source rate times slowing-down residence time, computes species critical + velocities from critical energies, evaluates a slowing-down-distribution + fast-ion pressure for each species, and converts that pressure into a + pressure-equivalent mean fast-ion energy. Parameters ---------- - nd : - Thermal deuterium density (m^-3). - nt : - Thermal tritium density (m^-3). e_beam_kev : Beam energy (keV). critical_energy_deuterium : - Critical energy for electron/ion slowing down of the beam ion (deuterium neutral beam) (keV). + Critical energy for electron/ion slowing down of deuterium beam ions (keV). critical_energy_tritium : - Critical energy for beam slowing down (tritium neutral beam) (keV). - beam_slow_time : - Beam ion slowing down time on electrons (s). + Critical energy for electron/ion slowing down of tritium beam ions (keV). + t_beam_slow : + Beam ion slowing-down time on electrons (s). f_beam_tritium : - Beam tritium fraction (0.0 = deuterium beam). + Tritium fraction in the injected beam (0.0 = pure deuterium beam). c_beam_total : - Beam current (A). - temp_plasma_ion_vol_avg_kev : - Thermal ion temperature (keV). + Total neutral beam current (A). vol_plasma : Plasma volume (m^3). - svdt : - Profile averaged for D-T (m^3/s). Returns ------- : - tuple[float, float, float, float]: A tuple containing the following elements: - - Alpha power from deuterium beam-background fusion (MW). - - Alpha power from tritium beam-background fusion (MW). - - Hot beam ion density (m^-3). - - Average hot beam ion energy (keV). + Beam slowing-down state containing deuterium and tritium hot beam-ion + densities, deuterium and tritium critical speeds, total hot beam-ion + density, and average deposited beam-ion energy. Notes ----- - - The function uses the Bosch-Hale parametrization to compute the reactivity. - The critical energy for electron/ion slowing down of the beam ion is calculated for both deuterium and tritium neutral beams. - The function integrates the hot beam fusion reaction rate integrand over the @@ -1119,13 +1155,13 @@ def beamcalc( beam_energy_ratio_deuterium = e_beam_kev / critical_energy_deuterium # Calculate the characterstic time for the deuterium ions to slow down to the thermal energy, eg E = 0. - characteristic_deuterium_beam_slow_time = ( - beam_slow_time / 3.0 * np.log(1.0 + (beam_energy_ratio_deuterium) ** 1.5) + characteristic_deuterium_t_beam_slow = ( + t_beam_slow / 3.0 * np.log(1.0 + (beam_energy_ratio_deuterium) ** 1.5) ) deuterium_beam_density = ( beam_current_deuterium - * characteristic_deuterium_beam_slow_time + * characteristic_deuterium_t_beam_slow / (constants.ELECTRON_CHARGE * vol_plasma) ) @@ -1134,17 +1170,17 @@ def beamcalc( # Calculate the characterstic time for the tritium to slow down to the thermal energy, eg E = 0. # Wesson, J. (2011) Tokamaks. - characteristic_tritium_beam_slow_time = ( - beam_slow_time / 3.0 * np.log(1.0 + (beam_energy_ratio_tritium) ** 1.5) + characteristic_tritium_t_beam_slow = ( + t_beam_slow / 3.0 * np.log(1.0 + (beam_energy_ratio_tritium) ** 1.5) ) tritium_beam_density = ( beam_current_tritium - * characteristic_tritium_beam_slow_time + * characteristic_tritium_t_beam_slow / (constants.ELECTRON_CHARGE * vol_plasma) ) - hot_beam_density = deuterium_beam_density + tritium_beam_density + nd_beam_hot = deuterium_beam_density + tritium_beam_density # Find the speed of the deuterium particle when it has the critical energy. # Re-arrange kinetic energy equation to find speed. Non-relativistic. @@ -1175,7 +1211,7 @@ def beamcalc( pressure_coeff_deuterium = ( constants.M_DEUTERON_AMU * constants.ATOMIC_MASS_UNIT - * beam_slow_time + * t_beam_slow * deuterium_critical_energy_speed**2 * source_deuterium / (constants.KILOELECTRON_VOLT * 3.0) @@ -1183,7 +1219,7 @@ def beamcalc( pressure_coeff_tritium = ( constants.M_TRITON_AMU * constants.ATOMIC_MASS_UNIT - * beam_slow_time + * t_beam_slow * tritium_critical_energy_speed**2 * source_tritium / (constants.KILOELECTRON_VOLT * 3.0) @@ -1200,64 +1236,55 @@ def beamcalc( # Beam deposited energy # Find the energy from the ideal gas pressure, P=1/3 * nmv^2 = 2/3 * n - deuterium_deposited_energy = 1.5 * deuterium_pressure / deuterium_beam_density - tritium_deposited_energy = 1.5 * tritium_pressure / tritium_beam_density - - total_depsoited_energy = ( - (deuterium_beam_density * deuterium_deposited_energy) - + (tritium_beam_density * tritium_deposited_energy) - ) / hot_beam_density - - hot_deuterium_rate = 1e-4 * beam_reaction_rate( - constants.M_DEUTERON_AMU, deuterium_critical_energy_speed, e_beam_kev + deuterium_deposited_energy = ( + 1.5 * deuterium_pressure / deuterium_beam_density + if deuterium_beam_density > 0.0 + else 0.0 ) - - hot_tritium_rate = 1e-4 * beam_reaction_rate( - constants.M_TRITON_AMU, tritium_critical_energy_speed, e_beam_kev + tritium_deposited_energy = ( + 1.5 * tritium_pressure / tritium_beam_density + if tritium_beam_density > 0.0 + else 0.0 ) - deuterium_beam_alpha_power = alpha_power_beam( - deuterium_beam_density, - nt, - hot_deuterium_rate, - vol_plasma, - temp_plasma_ion_vol_avg_kev, - svdt, - ) - tritium_beam_alpha_power = alpha_power_beam( - tritium_beam_density, - nd, - hot_tritium_rate, - vol_plasma, - temp_plasma_ion_vol_avg_kev, - svdt, + e_beam_deposited_kev = ( + ( + deuterium_beam_density * deuterium_deposited_energy + + tritium_beam_density * tritium_deposited_energy + ) + / nd_beam_hot + if nd_beam_hot > 0.0 + else 0.0 ) - return ( - deuterium_beam_alpha_power, - tritium_beam_alpha_power, - hot_beam_density, - total_depsoited_energy, + return BeamSlowingDownState( + deuterium_beam_density=deuterium_beam_density, + tritium_beam_density=tritium_beam_density, + deuterium_critical_energy_speed=deuterium_critical_energy_speed, + tritium_critical_energy_speed=tritium_critical_energy_speed, + nd_beam_hot=nd_beam_hot, + e_beam_deposited_kev=e_beam_deposited_kev, ) def fast_ion_pressure_integral(e_beam_kev: float, critical_energy: float) -> float: - """Calculate the fraction of initial beam energy given to the ions. + """Evaluate the dimensionless fast-ion pressure integral. - This function computes the fraction of initial beam energy given to the ions. based on the neutral beam energy - and the critical energy for electron/ion slowing down of the beam ion. + This function returns the closed-form dimensionless integral factor + appearing in the analytic fast-ion pressure expression for a + monoenergetic source with a classical slowing-down distribution. Parameters ---------- e_beam_kev : - Neutral beam energy (keV). + Beam birth energy (keV). critical_energy : - Critical energy for electron/ion slowing down of the beam ion (keV). + Critical energy for electron/ion slowing down (keV). Returns ------- : - float: Fraction of initial beam energy given to the ions. + Dimensionless pressure integral factor. Notes ----- @@ -1287,130 +1314,151 @@ def fast_ion_pressure_integral(e_beam_kev: float, critical_energy: float) -> flo return t1 + t2 - t3 - t4 -def alpha_power_beam( - beam_ion_desnity: float, - plasma_ion_desnity: float, - sigv: float, +def beam_target_reaction_rate( + nd_beam_ion: float, + nd_target_ion: float, + sigv_beam: float, vol_plasma: float, - temp_plasma_ion_vol_avg_kev: float, - sigmav_dt: float, ) -> float: - """Calculate alpha power from beam-background fusion. + """Calculate the total beam-target fusion reaction rate. - This function computes the alpha power generated from the interaction between - hot beam ions and thermal ions in the plasma. + This function computes the total number of beam-target fusion reactions + occurring per second in the plasma. + + The total beam-target reaction rate is given by: + reaction_rate = nd_beam_ion * nd_target_ion * sigv_beam * vol_plasma + + Unit analysis: + (m^-3) * (m^-3) * (m^3/s) * (m^3) = s^-1 Parameters ---------- - beam_ion_desnity : + nd_beam_ion : Hot beam ion density (m^-3). - plasma_ion_desnity : - Thermal ion density (m^-3). - sigv : - Hot beam fusion reaction rate (m^3/s). + nd_target_ion : + Thermal target ion density (m^-3). + sigv_beam : + Beam-target fusion reactivity (m^3/s). vol_plasma : Plasma volume (m^3). - temp_plasma_ion_vol_avg_kev : - Thermal ion temperature (keV). - sigmav_dt : - Profile averaged for D-T (m^3/s). Returns ------- : - float: Alpha power from beam-background fusion (MW). + Total beam-target fusion reaction rate (s^-1). Notes ----- - - The function uses the Bosch-Hale parametrization to compute the reactivity. - - The ratio of the profile-averaged to the reactivity at the given - thermal ion temperature is used to scale the alpha power. + Any spatial/profile effects should be introduced through the + supplied nd_target_ion. + """ + return nd_beam_ion * nd_target_ion * sigv_beam * vol_plasma - References + +def alpha_power_beam( + beam_target_reaction_rate: float, +) -> float: + """Calculate alpha power from beam-target fusion. + + This function converts the total beam-target fusion reaction rate into + alpha particle power. + + The alpha power is given by: + P_alpha = beam_target_reaction_rate * DT_ALPHA_ENERGY + + Unit analysis: + (s^-1) * (J) = W + + The result is converted from watts to megawatts (MW). + + Parameters ---------- - - H.-S. Bosch and G. M. Hale, “Improved formulas for fusion cross-sections and thermal reactivities,” - Nuclear Fusion, vol. 32, no. 4, pp. 611-631, Apr. 1992, - doi: https://doi.org/10.1088/0029-5515/32/4/i07. - """ - # Calculate the reactivity ratio - ratio = ( - sigmav_dt - / bosch_hale_reactivity( - np.array([temp_plasma_ion_vol_avg_kev]), - BoschHaleConstants(**REACTION_CONSTANTS_DT), - ).item() - ) + beam_target_reaction_rate : + Total beam-target fusion reaction rate (s^-1). - # Calculate and return the alpha power - return ( - beam_ion_desnity - * plasma_ion_desnity - * sigv - * (constants.DT_ALPHA_ENERGY / 1e6) - * vol_plasma - * ratio - ) + Returns + ------- + : + Alpha power from beam-target fusion (MW). + """ + return beam_target_reaction_rate * constants.DT_ALPHA_ENERGY / 1.0e6 -def beam_reaction_rate( - relative_mass_ion: float, critical_velocity: float, beam_energy_keV: float +def beam_reaction_rate_coefficient( + relative_mass_ion: float, critical_velocity: float, beam_energy_kev: float ) -> float: - """Calculate the hot beam fusion reaction rate. + """Calculate an effective beam-target fusion rate coefficient. - This function computes the fusion reaction rate for hot beam ions - using the critical velocity for electron/ion slowing down and the - neutral beam energy. + This function evaluates an effective fusion rate coefficient for a + slowing-down fast-ion population interacting with a thermal target + plasma. The beam ions are assumed to follow a classical slowing-down + distribution characterised by the critical velocity, and the fusion + contribution is obtained by integrating a beam-energy-dependent fusion + cross-section over the slowing-down distribution. Parameters ---------- relative_mass_ion : - Relative atomic mass of the ion (e.g., approx 2.0 for D, 3.0 for T). + Relative atomic mass of the ion (e.g. approx 2.0 for D, 3.0 for T). critical_velocity : Critical velocity for electron/ion slowing down of the beam ion [m/s]. - beam_energy_keV : + beam_energy_kev : Neutral beam energy [keV]. Returns ------- : - float: Hot beam fusion reaction rate (m^3/s). + Hot beam fusion rate coefficient [m^3/s]. Notes ----- - - The function integrates the hot beam fusion reaction rate integrand - over the range of beam velocities up to the critical velocity. - - The integration is performed using the quad function from scipy.integrate. + - The function returns a reactivity rate coefficient, not the reaction rate itself. + - The integration variable is the beam speed normalised to the critical speed. + - The underlying beam fusion cross-section is evaluated in cm^2 and converted + internally to m^2 so that this function returns an SI rate coefficient. + - The integration is performed using scipy.integrate.quad. + References + ---------- + - S. Niikura and M. Nagami, "Improvement of fusion reactivity and fusion + power multiplication factor in the presence of fast ions," Fusion + Engineering and Design, vol. 12, pp. 467-480, 1990. """ # Find the speed of the beam particle when it has the critical energy. # Re-arrange kinetic energy equation to find speed. Non-relativistic. beam_velocity = np.sqrt( - (beam_energy_keV * constants.KILOELECTRON_VOLT) + (beam_energy_kev * constants.KILOELECTRON_VOLT) * 2.0 / (relative_mass_ion * constants.ATOMIC_MASS_UNIT) ) relative_velocity = beam_velocity / critical_velocity - integral_coefficient = 3.0 * critical_velocity / np.log(1.0 + (relative_velocity**3)) + integral_coefficient = 3.0 * critical_velocity / np.log(1.0 + relative_velocity**3) - fusion_integral = integrate.quad( + fusion_integral_cm2 = integrate.quad( _hot_beam_fusion_reaction_rate_integrand, 0.0, relative_velocity, args=(critical_velocity,), )[0] - return integral_coefficient * fusion_integral + # _hot_beam_fusion_cross_section returns cm^2, so the integrated quantity is in + # cm^2. Convert to m^2 here so the final rate coefficient is returned in m^3/s. + fusion_integral_m2 = fusion_integral_cm2 * 1.0e-4 + + return integral_coefficient * fusion_integral_m2 def _hot_beam_fusion_reaction_rate_integrand( velocity_ratio: float, critical_velocity: float ) -> float: - """Integrand function for the hot beam fusion reaction rate. + """Evaluate the integrand for the beam-target fusion rate coefficient. - This function computes the integrand for the hot beam fusion reaction rate - based on the ratio of beam velocity to the critical velocity and the critical - velocity for electron/ion slowing down of the beam ion. + This function returns the slowing-down-weighted integrand used to + compute the effective beam-target fusion rate coefficient for a fast-ion + population. The integration variable is the beam speed normalised to the + critical speed, and the integrand combines a slowing-down weighting with + the beam fusion cross-section evaluated at the corresponding beam energy. Parameters ---------- @@ -1423,14 +1471,22 @@ def _hot_beam_fusion_reaction_rate_integrand( Returns ------- : - float: Value of the integrand for the hot beam fusion reaction rate. + Value of the integrand for the hot beam fusion reaction rate. Notes ----- - - The function uses the ratio of the beam velocity to the critical velocity - to compute the integrand. - - The integrand involves the fusion reaction cross-section and the critical - velocity for electron/ion slowing down of the beam ion. + - The factor `velocity_ratio**3 / (1.0 + velocity_ratio**3)` represents + the slowing-down weighting in the normalised velocity variable. + - The beam energy passed to `_beam_fusion_cross_section` is constructed + from the instantaneous beam speed and expressed in keV/amu. + - This function is an internal helper used by + `beam_reaction_rate_coefficient`. + + References + ---------- + - S. Niikura and M. Nagami, "Improvement of fusion reactivity and fusion + power multiplication factor in the presence of fast ions," Fusion + Engineering and Design, vol. 12, pp. 467-480, 1990. """ intgeral_term = (velocity_ratio**3) / (1.0 + velocity_ratio**3) @@ -1448,7 +1504,7 @@ def _hot_beam_fusion_reaction_rate_integrand( def _beam_fusion_cross_section(vrelsq: float) -> float: - """Calculate the fusion reaction cross-section. + """Calculate the beam fusion cross-section from beam energy. This function computes the fusion reaction cross-section based on the square of the speed of the beam ion (keV/amu). The functional form of @@ -1468,11 +1524,12 @@ def _beam_fusion_cross_section(vrelsq: float) -> float: Notes ----- - The cross-section is limited at low and high beam energies. - For beam kinetic energy less than 10 keV, the cross-section is set to 1.0e-27 cm^2. - For beam kinetic energy greater than 10,000 keV, the cross-section is set to 8.0e-26 cm^2. - The cross-section is calculated using a functional form with parameters a1 to a5. - + - The cross-section is limited at low and high beam energies. + - For beam kinetic energy less than 10 keV, the cross-section is set to 1.0e-27 cm^2. + - For beam kinetic energy greater than 10,000 keV, the cross-section is set to 8.0e-26 cm^2. + - The cross-section is calculated using a functional form with parameters a1 to a5. + - The exact provenance of this particular fit and its coefficients has + not yet been traced in the present code documentation. """ a1 = 45.95 a2 = 5.02e4 diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 8f1e10ccc0..1c3b0c6286 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -568,7 +568,6 @@ def run(self): # Calculate neutral beam slowing down effects # If ignited, then ignore beam fusion effects - if (current_drive_variables.c_beam_total != 0.0e0) and ( # noqa: RUF069 physics_variables.i_plasma_ignited == 0 ): @@ -579,8 +578,7 @@ def run(self): ) = reactions.beam_fusion( physics_variables.beamfus0, physics_variables.betbm0, - physics_variables.b_plasma_surface_poloidal_average, - physics_variables.b_plasma_toroidal_on_axis, + physics_variables.b_plasma_total, current_drive_variables.c_beam_total, physics_variables.nd_plasma_electrons_vol_avg, physics_variables.nd_plasma_fuel_ions_vol_avg, @@ -589,9 +587,7 @@ def run(self): physics_variables.f_plasma_fuel_deuterium, physics_variables.f_plasma_fuel_tritium, current_drive_variables.f_beam_tritium, - physics_variables.sigmav_dt_average, physics_variables.temp_plasma_electron_density_weighted_kev, - physics_variables.temp_plasma_ion_density_weighted_kev, physics_variables.vol_plasma, physics_variables.n_charge_plasma_effective_mass_weighted_vol_avg, ) diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py index 7761b57b31..e158edd95c 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -2040,8 +2040,7 @@ def st_phys(self, output): ) = reactions.beam_fusion( physics_variables.beamfus0, physics_variables.betbm0, - physics_variables.b_plasma_surface_poloidal_average, - physics_variables.b_plasma_toroidal_on_axis, + physics_variables.b_plasma_total, current_drive_variables.c_beam_total, physics_variables.nd_plasma_electrons_vol_avg, physics_variables.nd_plasma_fuel_ions_vol_avg, @@ -2050,9 +2049,7 @@ def st_phys(self, output): physics_variables.f_plasma_fuel_deuterium, physics_variables.f_plasma_fuel_tritium, current_drive_variables.f_beam_tritium, - physics_variables.sigmav_dt_average, physics_variables.temp_plasma_electron_density_weighted_kev, - physics_variables.temp_plasma_ion_density_weighted_kev, physics_variables.vol_plasma, physics_variables.n_charge_plasma_effective_mass_weighted_vol_avg, ) diff --git a/tests/unit/models/physics/test_fusion_reactions.py b/tests/unit/models/physics/test_fusion_reactions.py index e34b44e990..f85557a20e 100644 --- a/tests/unit/models/physics/test_fusion_reactions.py +++ b/tests/unit/models/physics/test_fusion_reactions.py @@ -5,6 +5,7 @@ import numpy as np import pytest +from process.core import constants from process.data_structure import physics_variables as pv from process.models.physics import fusion_reactions as reactions @@ -232,8 +233,7 @@ def test_beam_fusion(): beta_beam, nd_beam_ions_out, p_beam_alpha_mw = reactions.beam_fusion( 1.0, 1.5, - 0.85, - 5.3, + 5.367727, 130, 7.8e19, 6.6e19, @@ -242,8 +242,6 @@ def test_beam_fusion(): 0.5, 0.5, 1e-06, - 2.8e-22, - 13.5, 13.5, 1888.0, 0.425, @@ -251,33 +249,28 @@ def test_beam_fusion(): assert beta_beam == pytest.approx(0.0026264022466211366) assert nd_beam_ions_out == pytest.approx(4.2133504058678246e17) - assert p_beam_alpha_mw == pytest.approx(11.593221085189192) + assert p_beam_alpha_mw == pytest.approx(9.271206216041564) -def test_beamcalc(): - ( - deuterium_beam_alpha_power, - tritium_beam_alpha_power, - hot_beam_density, - beam_deposited_energy, - ) = reactions.beamcalc( - 3.3e19, - 3.3e19, +def test_beam_slowing_down_state(): + beam_state = reactions.beam_slowing_down_state( 1000.0, 276.7, 415.0, 1.42, 1e-06, 130, - 13.5, 1888.0, - 2.8e-22, ) - assert deuterium_beam_alpha_power == pytest.approx(11.561197668655383) - assert tritium_beam_alpha_power == pytest.approx(1.0434445041616093e-05) - assert hot_beam_density == pytest.approx(4.1968331737565126e17) - assert beam_deposited_energy == pytest.approx(445.05787301616635) + assert beam_state.deuterium_beam_density == pytest.approx(4.1968331737565126e17) + assert beam_state.tritium_beam_density == pytest.approx(316553077182.4059, rel=1e-6) + assert beam_state.deuterium_critical_energy_speed == pytest.approx( + 5.1495e6, rel=1e-4 + ) + assert beam_state.tritium_critical_energy_speed == pytest.approx(5.1534e6, rel=1e-4) + assert beam_state.nd_beam_hot == pytest.approx(4.1968331737565126e17) + assert beam_state.e_beam_deposited_kev == pytest.approx(445.05787301616635) def test__fast_ion_pressure_integral(): @@ -286,15 +279,28 @@ def test__fast_ion_pressure_integral(): assert pressure_integral == pytest.approx(1.1061397270783706) -def test_alpha_power_beam(): - alpha_power_beam = reactions.alpha_power_beam( - 316000000000, 3.3e19, 7.5e-22, 1888.0, 13.5, 2.8e-22 +def test_beam_target_reaction_rate(): + reaction_rate = reactions.beam_target_reaction_rate( + nd_beam_ion=3.16e11, + nd_target_ion=3.3e19, + sigv_beam=7.5e-22, + vol_plasma=1888.0, ) - assert alpha_power_beam == pytest.approx(1.047572705194316e-05) + assert reaction_rate == pytest.approx(1.4766048e13) -def test_beam_reaction_rate(): - beam_reaction_rate = reactions.beam_reaction_rate(3.01550071597, 5140000.0, 1000.0) +def test_alpha_power_beam(): + beam_target_rate = 1.0e13 # s^-1 + result = reactions.alpha_power_beam(beam_target_rate) + expected = beam_target_rate * constants.DT_ALPHA_ENERGY / 1.0e6 + + assert result == pytest.approx(expected) + + +def test_beam_reaction_rate_coefficient(): + beam_reaction_rate = reactions.beam_reaction_rate_coefficient( + 3.01550071597, 5140000.0, 1000.0 + ) - assert beam_reaction_rate == pytest.approx(7.465047902975452e-18) + assert beam_reaction_rate == pytest.approx(7.465047902975452e-22)