Particle balance constraint updated

documentation/source/physics-models/plasma_fuelling.md

@@ -0,0 +1,244 @@
+# Plasma Fuelling | `PlasmaFuelling()`
+
+## Particle balance
+
+The control of fuelling is governed by 4 key particle flux equations for each of the primary fuel species and the thermal helium $\alpha$ ash.
+
+$$
+\frac{dN_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{T,beam}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
+$$
+
+$$
+\frac{dN_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{D,beam}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
+$$
+
+$$
+\frac{dN_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}} -\Gamma_{\text{D+3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
+$$
+
+$$
+\frac{dN_{\alpha,\text{thermal}}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha,\text{thermal}}}{\tau_{\alpha}^*}
+$$
+
+In a steady state equilibrium all 4 of these equations should balance, therefore:
+
+$$
+\frac{dN_{\text{D}}}{dt} = \frac{dN_{\text{T}}}{dt} = \frac{dN_{\text{3He}}}{dt} = \frac{dN_{\alpha,\text{thermal}}}{dt} = 0
+$$
+
+Here, $\eta_{\text{fuelling}}$ is the fuelling efficiency, which quantifies how effectively fuel injected into the vacuum vessel reaches the plasma core. The value of $\eta_{\text{fuelling}}$ depends on the fuelling method. Typical values are around 0.01--0.1 for low-field-side gas puffing, 0.1--0.2 for supersonic gas injection, and 0.5--0.9 for pellet injection, which can approach unity under favourable conditions. $\Gamma_{\text{fuelling}}$ is the fuel injection rate into the vacuum vessel. The product $\eta_{\text{fuelling}} \Gamma_{\text{fuelling}}$ therefore represents the effective fuelling rate, i.e. the fraction of the injected fuel that successfully penetrates into the plasma and becomes available for fusion reactions. $\Gamma_{\text{T,beam}}$ and $\Gamma_{\text{D,beam}}$ respectively are the particle source rates from neutral beam systems (if present).
+
+
+The fuelling fractional compositions is given by $f$
+
+ - $N$ is the total amount of ions in the plasma.
+
+ - $\tau_{\text{fuel}}^*$ is the recycling corrected fuel particle confinement time given by:
+
+    $\tau_{\text{fuel}}^* = (\tau_p) /  (1-R)$
+
+The (effective) exhaust efficiency is is given by , $\eta_{\text{eff}} = 1- R$ 
+
+The factor $\frac{R}{1-R}$ is the mean number of recycling events back into the burning region experieced by a particle before it is pumped away.
+
+The definition of the recycling coefficient $R = 1- \frac{\Gamma_{\text{pumps}}}{\Gamma_{\text{out}}}$, where $\Gamma_{\text{pumps}}$ is the number of particles exhausted by the pumps per second and $\Gamma_{\text{out}}$ is the number of particles per second transported radially outwards across the separatrix.
+
+
+
+Where $\tau_p$ is the particle confinement time which we can assume is approximately equal to the energy confinement time ($\tau_p = \tau_E$). 
+
+!!! warning "Relation between $\tau_p$ and $\tau_E$"
+
+    In this model the "raw" particle confinement time ($\tau_p$) is set to always match the energy confinement time ($\tau_E$). The only variation of this in terms of the "effective" particle confinement time $\tau_{p}^*$ is via the recycling coefficient ($R$). If needed new coeffcients may be added to scale ($\tau_p$) before recyling corrections if needed.
+
+!!! note "Quantifying $R$"
+
+    The recycling coefficient $R$, defined as the fraction of particles crossing the LCFS that return to the plasma, can depend on numerous factors—including vessel pumping speed, neutral pressure in the private‑divertor region, impurity seeding levels, and the detailed properties of the SOL. Among these parameters, $R$ is the least certain and the most difficult to quantify. In next‑step devices, the SOL temperature is expected to be high, so particles reflected from the vessel walls are mostly ionized within the SOL and are removed by pumping before they can effectively refuel the burning plasma. As a result, the recycling coefficient is anticipated to be lower than in present‑day tokamaks, where $R$ can often approach unity. An additional uncertainty is the extent of neutral penetration at the plasma edge, which influences both the pedestal density and the density profile, and therefore also affects $R$[^1].
+
+
+--------------
+
+### Tritium Flow Rate | `calculate_plasma_tritium_flow_rate()`
+
+$$
+\frac{dN_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{T,beam}} + \Gamma_{\text{D+D} \rightarrow \text{T}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{T}}^*}
+$$
+
+---------------
+
+### Deuterium Flow Rate | `calculate_plasma_deuterium_flow_rate()`
+
+$$
+\frac{dN_{\text{D}}}{dt} = f_{\text{fuelling,D}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}}+\Gamma_{\text{D,beam}} -2 \Gamma_{\text{D+D}}- \Gamma_{\text{D+3He}} - \Gamma_{\text{D+T}} - \frac{N_{\text{T}}}{\tau_{\text{D}}^*}
+$$
+
+---------------
+
+### Helium-3 Flow Rate | `calculate_plasma_helium3_flow_rate()`
+
+$$
+\frac{dN_{\text{3He}}}{dt} = f_{\text{fuelling,3He}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \Gamma_{\text{D+D} \rightarrow \text{3He}}-\Gamma_{\text{D+3He}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*}
+$$
+
+---------------
+
+### Thermal Alpha Particle Flow Rate | `calculate_plasma_alphas_thermal_flow_rate()`
+
+$$
+\frac{dN_{\alpha,\text{thermal}}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha,\text{thermal}}}{\tau_{\alpha}^*}
+$$
+
+-----------------
+
+## Fuel Burnup Fration
+
+In a steady state tokamak, the burnup fraction ($f_b$) is explicitly defined by the following rate equation:
+
+$$
+f_b = \frac{\Gamma_{\text{fusion}}}{\sum\Gamma_{\text{fuel}}}
+$$
+
+where $\Gamma_{\text{fusion}}$ is the fusion reaction rate per second $[\text{s}^{-1}]$, and $\sum\Gamma_{\text{fuel}}$ is the sum of all fuel injection rates into the vacuum vessel by any means in particles per second $[\text{s}^{-1}]$.
+
+-----------------
+
+### Total Fuel Burnup Fraction | `calculate_fuel_burnup_fraction()`
+
+For the total burnup fraction we state:
+
+$$
+\overbrace{f_b}^{\texttt{f_plasma_fuel_burnup}} = \frac{2\left(\Gamma_{\text{D+D}}+\Gamma_{\text{D+T}}+\Gamma_{\text{D+3He}}\right)}{\Gamma_{\text{fuel}}+\Gamma_{\text{D,beam}}+\Gamma_{\text{T,beam}}}
+$$
+
+Here the factor of 2 is included as each fusion reaction removes 2 particles but our fuelling rate looks at indivudal particles injected.
+
+-------------------
+
+### Tritium Burnup Fraction | `calculate_tritium_burnup_fraction()`
+
+For just the tritium burnup fraction we state:
+
+$$
+\overbrace{f_b}^{\texttt{f_plasma_tritium_burnup}} = \frac{\left(\Gamma_{\text{D+T}}\right)}{\Gamma_{\text{fuel}}f_{\text{fuelling,T}}+\Gamma_{\text{T,beam}}} 
+$$
+
+------------------
+
+### Deuterium Burnup Fraction | `calculate_deuterium_burnup_fraction()`
+
+For just the deuterium burnup fraction we state:
+
+$$
+\overbrace{f_b}^{\texttt{f_plasma_deuterium_burnup}} = \frac{2\left(\Gamma_{\text{D+D}}+\Gamma_{\text{D+3He}}\right)}{\Gamma_{\text{fuel}}f_{\text{fuelling,D}}+\Gamma_{\text{D,beam}}}
+$$
+
+------------------
+
+## Key Constraints
+
+The implementation of the equations above as equality constraints is required in order to achieve a plasma steady state equilibrium where there is not a net rate of change of species number. In any given run several of these constraints and iterations variables will almost always be on.Below is the generic input used for a D-T only plasma where we have allowed all of the key solution variables to be iteration variables: 
+
+```python
+* Tritium particle balance
+icc = 93
+
+* Deuterium particle balance
+icc = 94
+
+* Alpha particle balance
+icc = 96
+
+* Fuelling composition consistency
+icc = 97
+
+
+* Particle recycling fraction
+ixc = 178
+f_plasma_particles_lcfs_recycled = 0.9
+
+* Plasma fuelling efficiecy
+ixc = 179
+eta_plasma_fuelling = 0.7
+
+* Injected VV fuelling rate
+ixc = 180
+molflow_plasma_fuelling_vv_injected = 5e21
+boundl(180) = 1e20
+
+* Deuterium fuelling fraction
+ixc = 181
+f_molflow_plasma_fuelling_deuterium = 0.5
+boundl(181) = 0.4
+
+* Tritium fuelling fraction
+ixc = 182
+f_molflow_plasma_fuelling_tritium = 0.5
+boundl(182) = 0.4
+```
+
+Note that the Helium-3 consistency equation (`icc= 95`) is not active as we are not actively fuelling and have no stated Helium-3 in the plasma composition. 
+
+If solving this as an evaluation problem it would be typical to fix the values of the reycling fraction, $R$ (`f_plasma_particles_lcfs_recycled`) and the fuelling efficiency $\eta_{\text{fuelling}}$ (`eta_plasma_fuelling`) as these represent higher order assumptions already about the plant design that `PROCESS` cannot model. This realistically leaves the total fuelling rate $\Gamma_{\text{fuel}}$ (`molflow_plasma_fuelling_vv_injected`) and the fuelling fractions $f_{\text{fuelling,D}}$,$f_{\text{fuelling,T}}$ (`f_molflow_plasma_fuelling_deuterium,f_molflow_plasma_fuelling_tritium`) to solve this consistency equations. This is ideal as the fuelling rates and fractions are ideal optimisation paramters as they can be directly controlled from the machine control room.
+
+------------
+
+### Tritium Flow Consistency
+
+This constraint can be activated by stating `icc = 93` in the input file.
+
+This constraint ensures that the change in tritium particles as a function of time is zero. It ensures the output of `calculate_plasma_tritium_flow_rate()` is zero
+
+**It is recommended to have this constraint on as it is a plasma equilibrium solution model**
+
+-----------------
+
+### Deuterium Flow Consistency
+
+This constraint can be activated by stating `icc = 94` in the input file.
+
+This constraint ensures that the change in deuterium particles as a function of time is zero. It ensures the output of `calculate_plasma_deuterium_flow_rate()` is zero
+
+**It is recommended to have this constraint on as it is a plasma equilibrium solution model**
+
+----------------
+
+### Helium-3 Flow Consistency
+
+This constraint can be activated by stating `icc = 95` in the input file.
+
+This constraint ensures that the change in helium-3 particles as a function of time is zero. It ensures the output of `calculate_plasma_helium3_flow_rate()` is zero
+
+**It is recommended to have this constraint on as it is a plasma equilibrium solution model**
+
+----------------
+
+### Thermal Alpha Particle Flow Consistency
+
+This constraint can be activated by stating `icc = 96` in the input file.
+
+This constraint ensures that the change in thermal alpha particles as a function of time is zero. It ensures the output of `calculate_plasma_alphas_thermal_flow_rate()` is zero
+
+**It is recommended to have this constraint on as it is a plasma equilibrium solution model**
+
+------------------
+
+### Fuelling Proportion Consistency
+
+This constraint can be activated by stating `icc = 97` in the input file.
+
+This ensures that all 3 injected fuelling fractions sum up to 1:
+
+$$
+f_{\text{fuelling,D}} + f_{\text{fuelling,T}} + f_{\text{fuelling,3He}} = 1.0
+$$
+
+**It is recommended to have this constraint on as it is a plasma consistency model**
+
+-----------------
+
+
+[^1]: G. L. Jackson, V. S. Chan, and R. D. Stambaugh, “An Analytic Expression for the Tritium Burnup Fraction in Burning-Plasma Devices,” Fusion Science and Technology, vol. 64, no. 1, pp. 8–12, Jul. 2013, doi: https://doi.org/10.13182/fst13-a17042.
+
+[^2]: J. F. Artaud et al., “Metis: a fast integrated tokamak modelling tool for scenario design,” Nuclear Fusion, vol. 58, no. 10, pp. 105001–105001, Aug. 2018, doi: https://doi.org/10.1088/1741-4326/aad5b1.
+
+[^3]: D. Reiter, H. Kever, G. H. Wolf, M. Baelmans, R. Behrisch, and R. Schneider, “Helium removal from tokamaks,” Plasma Physics and Controlled Fusion, vol. 33, no. 13, pp. 1579–1600, Nov. 1991, doi: https://doi.org/10.1088/0741-3335/33/13/008.

mkdocs.yml

@@ -53,6 +53,7 @@ nav:
                   - Overview: physics-models/plasma_beta/plasma_beta.md
                   - Fast Alpha: physics-models/plasma_beta/plasma_alpha_beta_contribution.md
               - Density Limit: physics-models/plasma_density.md
+              - Fuelling: physics-models/plasma_fuelling.md
               - Composition & Impurities: physics-models/plasma_composition.md
               - Radiation: physics-models/plasma_radiation.md
               - Plasma Current:

process/core/constants.py

@@ -9,6 +9,12 @@
 MFILE = 13
 """Machine-optimised output file unit"""
 
+AVOGADRO_NUMBER = 6.02214076e23
+"""Avogadro's number [1/mol]
+Reference: National Institute of Standards and Technology (NIST)
+https://physics.nist.gov/cgi-bin/cuu/Value?na|search_for=avogadro
+"""
+
 ELECTRON_CHARGE = 1.602176634e-19
 """Electron / elementary charge [C]
 Reference: National Institute of Standards and Technology (NIST)

process/core/input.py

@@ -169,7 +169,6 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "beta_vol_avg_min": InputVariable("physics", float, range=(0.0, 1.0)),
     "betbm0": InputVariable("physics", float, range=(0.0, 10.0)),
     "b_plasma_toroidal_on_axis": InputVariable("physics", float, range=(0.0, 30.0)),
-    "burnup_in": InputVariable("physics", float, range=(0.0, 1.0)),
     "radius_plasma_core_norm": InputVariable(
         "impurity_radiation", float, range=(0.0, 1.0)
     ),
@@ -1182,6 +1181,22 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]:
     "f_len_sol_power_decay_inboard_outboard": InputVariable(
         "physics", float, range=(0.01, 2.0)
     ),
+    "f_plasma_particles_lcfs_recycled": InputVariable(
+        "physics", float, range=(0.0, 1.0)
+    ),
+    "eta_plasma_fuelling": InputVariable("physics", float, range=(0.0, 1.0)),
+    "molflow_plasma_fuelling_vv_injected": InputVariable(
+        "physics", float, range=(1e18, 1e24)
+    ),
+    "f_molflow_plasma_fuelling_deuterium": InputVariable(
+        "physics", float, range=(0.0, 1.0)
+    ),
+    "f_molflow_plasma_fuelling_tritium": InputVariable(
+        "physics", float, range=(0.0, 1.0)
+    ),
+    "f_molflow_plasma_fuelling_helium3": InputVariable(
+        "physics", float, range=(0.0, 1.0)
+    ),
 }