3385 Update the plasma composition docs to include all new models

documentation/proc-pages/development/standards.md

@@ -533,7 +533,7 @@ subroutine constraint_eqn_001(args)
   !! - \( T_i \) -- density weighted average ion temperature [keV]
   !! - \( B_{tot} \) -- total toroidal + poloidal field [T]
 
-  use physics_variables, only: beta_fast_alpha, beta_beam, dene, ten, dnitot, tin, btot, beta
+  use physics_variables, only: beta_fast_alpha, beta_beam, dene, ten, nd_ions_total, tin, btot, beta
   use constants, only: electron_charge,rmu0
 
   implicit none
@@ -542,7 +542,7 @@ subroutine constraint_eqn_001(args)
   !! constraint derived type
 
     args%cc = 1.0D0 - (beta_fast_alpha + beta_beam + &
-      2.0D3*rmu0*electron_charge * (dene*ten + dnitot*tin)/btot**2 )/beta
+      2.0D3*rmu0*electron_charge * (dene*ten + nd_ions_total*tin)/btot**2 )/beta
     args%con = beta * (1.0D0 - args%cc)
     args%err = beta * args%cc
     args%symbol = '='

documentation/proc-pages/eng-models/heating_and_current_drive/NBI/culham_nb.md

@@ -66,7 +66,7 @@ plus correction terms outlined in Culham Report AEA FUS 172.
 
 | Input       | Description                          |
 | :---------- | :----------------------------------- |
-| $\mathtt{abeam}$      | beam ion mass (amu)   |
+| $\mathtt{m_beam_amu}$      | beam ion mass (amu)   |
 | $\mathtt{alphan}$, $\alpha_n$       | density profile factor   |
 | $\mathtt{alphat}$, $\alpha_T$       | temperature profile factor  |
 |  $\mathtt{aspect}$, $A$      |   aspect ratio                            |
@@ -108,7 +108,7 @@ Critical energy ($\text{MeV}$) (power to electrons = power to ions) (IPDG89)
 N.B. ten is in keV
 
 $$
-\mathtt{ecrit} = 0.01 \times \mathtt{abeam} \times \mathtt{ten}
+\mathtt{ecrit} = 0.01 \times \mathtt{m_beam_amu} \times \mathtt{ten}
 $$
 
 Beam energy in MeV
@@ -128,7 +128,7 @@ $$
 $$
 
 $$
-\mathtt{yj} = \frac{0.8 \times Z_{\text{eff}}}{\mathtt{abeam}}
+\mathtt{yj} = \frac{0.8 \times Z_{\text{eff}}}{\mathtt{m_beam_amu}}
 $$
 
 Fitting function J0(x,y)

documentation/proc-pages/eng-models/heating_and_current_drive/NBI/iter_nb.md

@@ -59,7 +59,7 @@ ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990
 
 | Input   | Description                                               |
 |---------|-----------------------------------------------------------|
-| $\mathtt{abeam}$, $m_{\text{u,ion}}$   | Beam ion mass ($\text{amu}$)                                       |
+| $\mathtt{m_beam_amu}$, $m_{\text{u,ion}}$   | Beam ion mass ($\text{amu}$)                                       |
 | $\mathtt{alphan}$  | Density profile factor                                    |
 | $\mathtt{alphat}$  | Temperature profile factor                                |
 | $\mathtt{aspect}$, $A$  | Aspect ratio                                              |

documentation/proc-pages/eng-models/heating_and_current_drive/NBI/nbi_overview.md

@@ -60,6 +60,9 @@ with $E, n_e, T$ expressed in units of keV/u, $\text{cm}^3$ and keV, respectivel
 !!! info "Info" 
     For the full table of values for $A_{ijk}$ & $B_{ijk}^{(z)}$\  please see the accompanying paper[^1] or `current_drive.py`
 
+!!! warning "Model applicability with different impurities"
+    Be cautious that the model only incorporates the presence of carbon, oxygen and iron impurities. Using other impurities will mean that the cross-section is not calculated appropriately.
+
 For a plasma having an arbitrary mix of $N$ different types of impurities with densities $n$, and charges $Z_q (q = 1, ..., N)$, the beam stopping cross-section can be represented as the weighted sum of the stopping cross- sections for $N$ reference single-impurity plasmas. In each of these reference plasmas, the electron density and the proton density (including that of deuterium and tritium ions) are the same as in a true plasma. The impurity density, however, is increased in order to satisfy quasi-neutrality. The weighting function is the electron density $n_qZ_q$ associated with the amu impurity (in the true plasma), divided by the sum of these densities. The result is: 
 
 $$

documentation/proc-pages/physics-models/error.txt

@@ -470,8 +470,8 @@ density and temperature profiles\footnote{H. Lux, R. Kemp, D.J. Ward, M.
   \textbar{} Des. \textbf{101}, 42-51 (2015)}
 
 The impurity number density fractions relative to the electron density
-are set using input array \texttt{fimp(1,...,nimp)}, where
-\texttt{nimp\ =\ 14} is the number of impurity species in the radiation
+are set using input array \texttt{fimp(1,...,n_impurities)}, where
+\texttt{n_impurities\ =\ 14} is the number of impurity species in the radiation
 model. The available impurities are as follows:
 
 \begin{longtable}[]{@{}cl@{}}
@@ -500,7 +500,7 @@ As stated above, the number density fractions for hydrogen (all
 isotopes) and helium need not be set, as they are calculated by the code
 to ensure that plasma quasi-neutrality is maintained, and taking into
 account the fuel ratios \texttt{f_deuterium}, \texttt{f_tritium} and
-\texttt{f_helium3}, and the alpha particle fraction \texttt{ralpne} which may
+\texttt{f_helium3}, and the alpha particle fraction \texttt{f_nd_alpha_electron} which may
 be input by the user.
 
 The impurity fraction of one of the elements listed in array
@@ -509,7 +509,7 @@ The impurity fraction of one of the elements listed in array
 The synchrotron radiation power\footnote{Albajar, Nuclear Fusion
   \textbf{41} (2001) 665} \footnote{Fidone, Giruzzi and Granata, Nuclear
   Fusion \textbf{41} (2001) 1755} is assumed to originate from the
-plasma core. The wall reflection factor \texttt{ssync} may be set by the
+plasma core. The wall reflection factor \texttt{f_sync_reflect} may be set by the
 user.
 
 By changing the input parameter \texttt{coreradius}, the user may set
@@ -994,7 +994,7 @@ is derived directly from the energy confinement scaling law.
 \texttt{iradloss\ =\ 0} -- Total power lost is scaling power plus
 radiation
 
-\texttt{pscaling\ +\ pradpv\ =\ f_alpha_plasma*alpha_power_density\ +\ charged_power_density\ +\ pden_plasma_ohmic_mw\ +\ pinjmw/vol_plasma}
+\texttt{pscaling\ +\ pden_plasma_rad_mw\ =\ f_alpha_plasma*alpha_power_density\ +\ charged_power_density\ +\ pden_plasma_ohmic_mw\ +\ pinjmw/vol_plasma}
 
 \texttt{iradloss\ =\ 1} -- Total power lost is scaling power plus core
 radiation only

documentation/proc-pages/physics-models/fusion_reactions/beam_reactions.md

@@ -103,8 +103,8 @@ Please see the [H&CD section](../../eng-models/heating_and_current_drive/heating
 3. **Set the plasma tritium and ion densities**
 
     $$
-    \mathtt{deuterium\_density = deni * f\_deuterium\_plasma} \\
-    \mathtt{tritium\_density = deni * f\_tritium\_plasma}
+    \mathtt{deuterium\_density = nd_fuel_ions * f\_deuterium\_plasma} \\
+    \mathtt{tritium\_density = nd_fuel_ions * f\_tritium\_plasma}
     $$
 
 4. **Calculate the beam alpha powers, density and deposited energy**
@@ -385,7 +385,7 @@ Please see the [H&CD section](../../eng-models/heating_and_current_drive/heating
 
 This constraint can be activated by stating `icc = 7` in the input file.
 
-The desired value of the hot ion beam density calculated from the code (`beam_density_out`) can be constrained using the input variable, `rnbeam`. 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`.
+The desired value of the hot ion beam density calculated from the code (`beam_density_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`.
 
 [^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  

documentation/proc-pages/physics-models/fusion_reactions/plasma_reactions.md

@@ -56,7 +56,7 @@ profiles
 
 The fractional composition of the 'fuel' ions ($\text{D}$, $\text{T}$ and $^3\text{He}$) is
 controlled using the three variables `f_deuterium`, `f_tritium` and `f_helium3`, respectively.
-More information about setting seeded impurities and simulating first wall sputtering can be found in the [impurities and radiation section](../plasma_radiation_impurities.md)
+More information about setting seeded impurities and simulating first wall sputtering can be found in the [composition and impurities section](../plasma_composition.md)
 
 !!! note "Reactions not calculated"
 

documentation/proc-pages/physics-models/plasma_composition.md

@@ -0,0 +1,237 @@
+# Plasma composition and impurities.
+
+Within `PROCESS` we always assume the plasma is quasi-neutral eg:
+
+$$
+n_{\text{e}} = \underbrace{Z_{\text{fuel}}n_{\text{i}}}_{\text{Fuel Ions}} + \underbrace{2n_{\text{e}}f_{\alpha}}_{\text{Alpha particles}} + \underbrace{n_{\text{e}}f_{\text{beam}}}_{\text{Neutral beams}} + \underbrace{\sum_j Z_j n_{\text{e}} f_j}_{\text{Impurities}}
+$$
+
+* Since only deuterium and tritium can be placed into the beams the charge coefficient on the $n_{\text{e}}f_{\text{beam}}$ is just 1.
+
+--------------------
+
+## Setting the impurity composition
+
+
+
+The impurity number density fractions relative to the electron density are constant and are set 
+using input array `fimp(1,...,14)`. The available species along with their `fimp()` index and iteration variable number are as follows:
+
+
+- `fimp(1)`: Hydrogen isotopes (fraction calculated by code)
+- `fimp(2)`: Helium (fraction calculated by code)
+- `fimp(3)`: Beryllium, (`ixc = 125`)
+- `fimp(4)`: Carbon, (`ixc = 126`)
+- `fimp(5)`: Nitrogen, (`ixc = 127`)
+- `fimp(6)`: Oxygen, (`ixc = 128`)
+- `fimp(7)`: Neon, (`ixc = 129`)
+- `fimp(8)`: Silicon, (`ixc = 130`)
+- `fimp(9)`: Argon, (`ixc = 131`)
+- `fimp(10)`: Iron, (`ixc = 132`)
+- `fimp(11)`: Nickel, (`ixc = 133`)
+- `fimp(12)`: Krypton, (`ixc = 134`)
+- `fimp(13)`: Xenon, (`ixc = 135`)
+- `fimp(14)`: Tungsten, (`ixc = 136`)
+
+As stated above, the number density fractions for hydrogen (all isotopes) and
+helium should not be set, as they are calculated by the code. This is to ensure 
+plasma quasi-neutrality taking into account the fuel ratios
+`f_deuterium`, `f_tritium` and `f_helium3`, and the alpha particle fraction `f_nd_alpha_electron` which may 
+be input by the user or selected as an iteration variable.
+
+The impurity fraction of any one of the elements listed in array `fimp` (other than hydrogen 
+isotopes and helium) may be used as an iteration variable.
+
+**The impurity fraction to be varied can be set simply with `fimp(i) = <value>`, where `i` is the corresponding number value for the desired impurity in the table above.**
+
+----------------
+
+## Plasma Composition Calculation | `plasma_composition()`
+
+This function sets the fractional makeups and determines the relative density and charges of different plasma species.
+
+All of the plasma composites are normally given as a fraction of the volume averaged plasma electron density.
+
+1. **Alpha Ash Portion Calculation**
+
+    - Calculate the number density of alpha particles (`nd_alphas`) using the electron density (`dene`) and the alpha particle to electron ratio (`f_nd_alpha_electron`).
+        - `f_nd_alpha_electron` can be set as an iteration variable (`ixc = 109`) or set directly.
+
+    $$
+    n_{\alpha} = \mathtt{f\_nd\_alpha\_electron}\times n_{\text{e}}
+    $$
+
+
+2. **Protons Calculation**
+    - The calculation of proton density (`nd_protons`) depends on whether the alpha rate density has been calculated. This should only happen in the first function call as the rates are calculated later on in the code.
+
+    - If the alpha rate density is not yet calculated, use a rough estimate.
+
+    $$
+    \texttt{nd_protons} | n_{\text{p}} = \\
+    \text{max}\left[\texttt{f_nd_protium_electrons} \times n_{\text{e}}, n_{\alpha}\times \left(f_{\text{3He}} + 0.001\right)\right]
+    $$
+
+    - Otherwise, use the calculated proton rate density.
+
+    $$
+    \texttt{nd_protons} | n_{\text{p}} = \\
+    \text{max}\left[\texttt{f_nd_protium_electrons} \times n_{\text{e}}, n_{\alpha}\times \frac{r_{\text{p}}}{r_{\alpha,\text{total}}}\right]
+    $$
+
+    where $r_{\text{p}}$ is the rate of proton production and $r_{\alpha,\text{total}}$ is the rate of total alpha particle production, which includes beam fusion (if present).
+
+3. **Beam Hot Ion Component**
+
+    - Calculate the number density of beam ions (`nd_beam_ions`), using the electron density (`dene`) and the beam ion to electron ratio (`f_nd_beam_electron`). If the plasma is ignited, set it to zero.
+        - `f_nd_beam_electron` can be set as an iteration variable (`ixc = 7`) or set directly.
+
+    $$
+    \mathtt{nd\_beam\_ions} | n_{\text{beam}} = \mathtt{f\_nd\_beam\_electron} \times n_{\text{e}}
+    $$
+
+4. **Sum of charge number density for all impurity ions**
+
+    - Sum the product of charge number (`Zi`) and number density for all impurity ions with charge greater than helium. 
+
+    $$
+    \mathtt{znimp} = \sum_j Z_j n_{\text{e}} f_j
+    $$
+
+5. **Fuel Portion - Conserve Charge Neutrality**
+    - Calculate the fuel portion (`znfuel`) by conserving charge neutrality.
+
+    $$
+   \mathtt{znfuel} | \underbrace{Z_{\text{fuel}}n_{\text{i}}}_{\text{Fuel Ions}} =  n_{\text{e}} - \underbrace{2n_{\text{e}}f_{\alpha}}_{\text{Alpha particles}} - \underbrace{n_{\text{e}}f_{\text{beam}}}_{\text{Neutral beams}} - \underbrace{\sum_j Z_j n_{\text{e}} f_j}_{\text{Impurities}}
+    $$
+
+6. **Fuel Ion Density Calculation**
+
+    $$
+    \mathtt{nd\_fuel\_ions} | n_{\text{i}} = \frac{\mathtt{znfuel}}{1+f_{\text{3He}}}
+    $$
+
+    - Calculate the fuel ion density (`nd_fuel_ions`).
+
+7. **Set Hydrogen and Helium Impurity Fractions**
+
+    - Set the impurity fractions for hydrogen and helium species.
+
+    $$
+    \frac{n_{\text{H}}}{n_{\text{e}}} = n_{\text{protons}} + \left(f_{\text{deuterium}}+f_{\text{tritium}}\right)n_{\text{i}}+ n_{\text{beam}}
+    $$
+
+    $$
+    \frac{n_{\text{He}}}{n_{\text{e}}} = f_{\text{3He}}n_{\text{i}}+\mathtt{f\_nd\_alpha\_electron}
+    $$
+
+8. **Total Impurity Density Calculation**
+
+    - Calculate the total impurity density (`nd_impurities`).
+
+    $$
+    \mathtt{nd\_impurities} | n_{\text{impurities}} = \sum_j n_{\text{e}} f_j
+    $$
+
+9. **Total Ion Density Calculation**
+
+    - Calculate the total ion density (`nd_ions_total`).
+
+    $$
+    \mathtt{nd\_ions\_total} | n_{\text{i,total}} = n_{\text{i}} + n_{\alpha}+n_{\text{protons}}+ n_{\text{beam}}+n_{\text{impurities}}
+    $$
+
+10. **Set Impurity Fraction Variables**
+
+    - Set global impurity fraction variables for other routines.
+
+    $$
+    \mathtt{rncne} = \frac{n_{\text{C}}}{n_{\text{e}}}, \quad \mathtt{rnone} = \frac{n_{\text{O}}}{n_{\text{e}}} \quad \mathtt{rnfene} = \frac{n_{\text{Fe}}+ n_{\text{Ar}}}{n_{\text{e}}}
+    $$
+
+    The variable above are set as global physics variables to be used in the [`sigbeam()`](../eng-models/heating_and_current_drive/NBI/nbi_overview.md#beam-stopping-cross-section--sigbeam) routine, which calculates the beam stopping cross-section.
+
+11. **Effective Charge Calculation**
+
+    - Calculate the effective charge (`zeff`), which is given by:
+
+    $$
+    Z_{\text{eff}} = \frac{\sum_j Z^2_j n_{\text{i}}}{\sum_j Z_j n_{\text{i}}}
+    $$
+
+    As we assume the plasma is quai-neutral then:
+
+    $$
+    \sum_j Z_j n_{\text{i}} = n_{\text{e}}
+    $$
+
+    Thus we can write:
+
+    $$
+    \mathtt{zeff} | Z_{\text{eff}} = \frac{\sum_j Z^2_j n_{\text{e}} f_j}{n_{\text{e}}}
+    $$
+
+    More info can be found [here](https://wiki.fusion.ciemat.es/wiki/Effective_charge_state).
+
+
+12. **Fraction of Alpha Energy to Ions and Electrons**
+
+    - Calculate the fraction of alpha energy going to electrons and ions.
+
+    $$
+    \mathtt{f\_alpha\_electron} = 0.88155\exp{\left[-\langle T_{\text{e}} \rangle \frac{(1+\alpha_n)(1+\alpha_T)}{67.4036(1+\alpha_T+\alpha_n)}\right]}
+    $$
+
+    $$
+    \mathtt{f\_alpha\_ion} = (1.0 - \mathtt{f\_alpha\_electron})
+    $$
+
+13. **Average Atomic Masses of Injected Fuel Species**
+
+    - Calculate the average atomic masses of injected fuel species.
+
+    $$
+    \mathtt{m\_fuel\_amu} | m_{\text{fuel,amu}} = \left(m_{\text{D}}f_{\text{D}}\right) + \left(m_{\text{T}}f_{\text{T}}\right) + \left(m_{\text{3He}}f_{\text{3He}}\right)
+    $$
+
+14. **Average Atomic Masses of Injected Fuel Species in Neutral Beams**
+
+    - Calculate the average atomic masses of injected fuel species in the neutral beams.
+
+    $$
+    \mathtt{m\_beam\_amu} | m_{\text{beam,amu}} = \left(m_{\text{D}}(1-f_{\text{T,beam}}\right)) + \left(m_{\text{T}}f_{\text{T,beam}}\right)
+    $$
+
+15. **Density Weighted Mass Calculation**
+
+    - Calculate the density-weighted mass of ions.
+
+    $$
+    \mathtt{m\_ions\_total\_amu} = \\
+    \left[\frac{\left(m_{\text{fuel}}n_{\text{i}}\right) + \left(m_{\alpha}n_{\alpha}\right) + \left(m_{\text{p}}n_{p}\right) + \left(m_{\text{beam}}n_{\text{beam}}\right) + \sum_j m_{j} n_{\text{e}} f_j}{n_{\text{i,total}}}\right]
+    $$
+
+16. **Mass Weighted Plasma Effective Charge**
+
+    - Calculate the mass-weighted plasma effective charge (`zeffai`). Similar to the calculation of the effective charge except each element is divided by its mass
+
+    $$
+    \mathtt{zeffai} | Z_{\text{eff,m}} = \frac{\sum_j \frac{Z^2_j n_{\text{e}} f_j}{m_{\text{j}}}}{n_{\text{e}}}
+    $$
+
+---------------
+
+
+## Key Constraints
+
+
+### Reinke criterion, divertor impurity fraction lower limit
+
+This constraint can be activated by stating `icc = 78` in the input file.
+
+The minimum impurity fraction required from the Reinke module can be set with, `fzmin`
+
+The scaling value `freinke` can be varied also.
+
+
+[^1]: H. Lux, R. Kemp, D.J. Ward, M. Sertoli, Impurity radiation in DEMO systems modelling, Fusion Engineering and Design, Volume 101, 2015, Pages 42-51, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2015.10.002.

documentation/proc-pages/physics-models/plasma_confinement.md

@@ -78,7 +78,7 @@ derived directly from the energy confinement scaling law.
 
 `iradloss = 0` -- Total power lost is scaling power plus radiation:
 
-`pscaling + pradpv = f_alpha_plasma*alpha_power_density_total + charged_power_density + pden_plasma_ohmic_mw + pinjmw/vol_plasma`
+`pscaling + pden_plasma_rad_mw = f_alpha_plasma*alpha_power_density_total + charged_power_density + pden_plasma_ohmic_mw + pinjmw/vol_plasma`
 
 
 `iradloss = 1` -- Total power lost is scaling power plus radiation from a region defined as the "core":