diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md new file mode 100644 index 0000000000..9b2423ea28 --- /dev/null +++ b/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. \ No newline at end of file diff --git a/mkdocs.yml b/mkdocs.yml index 6cfa762788..cbdaf4def1 100644 --- a/mkdocs.yml +++ b/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: diff --git a/process/core/constants.py b/process/core/constants.py index 8583458d35..d6b8699de9 100644 --- a/process/core/constants.py +++ b/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) diff --git a/process/core/input.py b/process/core/input.py index fb282b26ec..7f8c515546 100644 --- a/process/core/input.py +++ b/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) + ), } diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index bd34580d43..ee9c1e50ce 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -68,6 +68,7 @@ ) from process.models.physics.density_limit import DensityLimitModel from process.models.physics.exhaust import calculate_brunner_divertor_power_splits +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.impurity_radiation import read_impurity_file from process.models.physics.l_h_transition import PlasmaConfinementTransitionModel from process.models.physics.physics import ( @@ -2942,7 +2943,7 @@ def plot_main_plasma_information( f"| D: {mfile.get('f_plasma_fuel_deuterium', scan=scan):.2f} | T: {mfile.get('f_plasma_fuel_tritium', scan=scan):.2f} | 3He: {mfile.get('f_plasma_fuel_helium3', scan=scan):.2f} |\n\n" f"Fusion Power, $P_{{\\text{{fus}}}}:$ {mfile.get('p_fusion_total_mw', scan=scan):,.2f} MW\n" f"D-T Power, $P_{{\\text{{fus,DT}}}}:$ {mfile.get('p_dt_total_mw', scan=scan):,.2f} MW\n" - f"D-D Power, $P_{{\\text{{fus,DD}}}}:$ {mfile.get('p_dd_total_mw', scan=scan):,.2f} MW\n" + f"D-D Power, $P_{{\\text{{fus,DD}}}}:$ {mfile.get('p_dd_total_mw', scan=scan):,.4f} MW\n" f"D-3He Power, $P_{{\\text{{fus,D3He}}}}:$ {mfile.get('p_dhe3_total_mw', scan=scan):,.2f} MW\n" f"Alpha Power, $P_{{\\alpha}}:$ {mfile.get('p_alpha_total_mw', scan=scan):,.2f} MW" ) @@ -2966,9 +2967,9 @@ def plot_main_plasma_information( f" - Average mass of all plasma ions: {mfile.get('m_ions_total_amu', scan=scan):.3f} amu\n" f"Fuel mass: {mfile.get('m_plasma_fuel_ions', scan=scan) * 1000:.4f} g\n" f" - Average mass of all fuel ions: {mfile.get('m_fuel_amu', scan=scan):.3f} amu\n\n" - f"Fueling rate: {mfile.get('molflow_plasma_fuelling_required', scan=scan):.3e} nucleus-pairs/s\n" - f"Fuel burn-up rate: {mfile.get('rndfuel', scan=scan):.3e} reactions/s \n" - f"Burn-up fraction: {mfile.get('burnup', scan=scan):.4f} \n" + f"Fueling rate: {mfile.get('molflow_plasma_fuelling_vv_injected', scan=scan):.3e} nucleus-pairs/s\n" + f"Fuel burn-up rate: {mfile.get('fusrat_total', scan=scan):.3e} reactions/s \n" + f"Burn-up fraction: {mfile.get('f_plasma_fuel_burnup', scan=scan):.4f} \n" ) axis.text( @@ -12113,8 +12114,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): # Add plasma volume, areas and shaping information textstr_general = ( f"Total fusion rate: {mfile.get('fusrat_total', scan=scan):.4e} reactions/s\n" - f"Total fusion rate density: {mfile.get('fusden_total', scan=scan):.4e} reactions/m3/s\n" - f"Plasma fusion rate density: {mfile.get('fusden_plasma', scan=scan):.4e} reactions/m3/s\n" + f"Total fusion rate density: {mfile.get('fusden_total', scan=scan):.4e} reactions/m$^3$/s\n" + f"Plasma fusion rate density: {mfile.get('fusden_plasma', scan=scan):.4e} reactions/m$^3$/s" ) axis.text( @@ -12136,9 +12137,11 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): # ============================================================================ textstr_dt = ( - f"Total fusion power: {mfile.get('p_dt_total_mw', scan=scan):,.2f} MW\n" - f"Plasma fusion power: {mfile.get('p_plasma_dt_mw', scan=scan):,.2f} MW \n" - f"Beam fusion power: {mfile.get('p_beam_dt_mw', scan=scan):,.2f} MW\n" + f"Total fusion power: {mfile.get('p_dt_total_mw', scan=scan):,.4f} MW\n" + f"Total fusion rate: {mfile.get('fusrat_dt_total', scan=scan):.4e} reactions/s\n" + f"Plasma fusion power: {mfile.get('p_plasma_dt_mw', scan=scan):,.4f} MW \n" + f"Plasma fusion rate: {mfile.get('fusrat_plasma_dt', scan=scan):.4e} reactions/s\n" + f"Beam fusion power: {mfile.get('p_beam_dt_mw', scan=scan):,.4f} MW" ) axis.text( @@ -12158,7 +12161,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): ) axis.text( - 0.24, + 0.285, 0.8, "$\\text{D - T}$", fontsize=20, @@ -12169,13 +12172,16 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): # ================================================= textstr_dd = ( - f"Total fusion power: {mfile.get('p_dd_total_mw', scan=scan):,.2f} MW\n" - f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f} \n" + f"Total fusion power: {mfile.get('p_dd_total_mw', scan=scan):,.4f} MW\n" + f"Tritium branching ratio: {mfile.get('f_dd_branching_trit', scan=scan):.4f} \n\n" + f"D+D -> T fusion rate: {mfile.get('fusrat_plasma_dd_triton', scan=scan):.4e} reactions/s \n" + f"D+D -> 3He fusion rate: {mfile.get('fusrat_plasma_dd_helion', scan=scan):.4e} reactions/s \n" + f"Total D-D fusion rate: {mfile.get('fusrat_plasma_dd_total', scan=scan):.4e} reactions/s" ) axis.text( 0.05, - 0.65, + 0.625, textstr_dd, fontsize=9, verticalalignment="bottom", @@ -12190,8 +12196,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): ) axis.text( - 0.22, - 0.685, + 0.31, + 0.69, "$\\text{D - D}$", fontsize=20, verticalalignment="top", @@ -12200,8 +12206,10 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): # ================================================= - textstr_dhe3 = f"Total fusion power: {mfile.get('p_dhe3_total_mw', scan=scan):,.2f} MW \n\n" - + textstr_dhe3 = ( + f"Total fusion power: {mfile.get('p_dhe3_total_mw', scan=scan):,.4f} MW \n" + f"D+3He fusion rate: {mfile.get('fusrat_plasma_dhe3', scan=scan):.4e} reactions/s \n" + ) axis.text( 0.05, 0.55, @@ -12219,8 +12227,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): ) axis.text( - 0.21, - 0.59, + 0.285, + 0.58, "$\\text{D - 3He}$", fontsize=20, verticalalignment="top", @@ -12230,15 +12238,15 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): # ================================================= textstr_alpha = ( - f"Total power: {mfile.get('p_alpha_total_mw', scan=scan):.2f} MW\n" - f"Plasma power: {mfile.get('p_plasma_alpha_mw', scan=scan):.2f} MW\n" - f"Beam power: {mfile.get('p_beam_alpha_mw', scan=scan):.2f} MW\n\n" - f"Rate density total: {mfile.get('fusden_alpha_total', scan=scan):.4e} particles/m3/sec\n" - f"Rate density, plasma: {mfile.get('fusden_plasma_alpha', scan=scan):.4e} particles/m3/sec\n\n" - f"Total power density: {mfile.get('pden_alpha_total_mw', scan=scan):.4e} MW/m3\n" - f"Plasma power density: {mfile.get('pden_plasma_alpha_mw', scan=scan):.4e} MW/m3\n\n" - f"Power per unit volume transferred to electrons: {mfile.get('f_pden_alpha_electron_mw', scan=scan):.4e} MW/m3\n" - f"Power per unit volume transferred to ions: {mfile.get('f_pden_alpha_ions_mw', scan=scan):.4e} MW/m3\n\n" + f"Total power: {mfile.get('p_alpha_total_mw', scan=scan):.4f} MW\n" + f"Plasma power: {mfile.get('p_plasma_alpha_mw', scan=scan):.4f} MW\n" + f"Beam power: {mfile.get('p_beam_alpha_mw', scan=scan):.4f} MW\n\n" + f"Rate density total: {mfile.get('fusden_alpha_total', scan=scan):.4e} particles/m$^3$/s\n" + f"Rate density, plasma: {mfile.get('fusden_plasma_alpha', scan=scan):.4e} particles/m$^3$/s\n\n" + f"Total power density: {mfile.get('pden_alpha_total_mw', scan=scan):.4e} MW/m$^3$\n" + f"Plasma power density: {mfile.get('pden_plasma_alpha_mw', scan=scan):.4e} MW/m$^3$\n\n" + f"Power per unit volume transferred to electrons: {mfile.get('f_pden_alpha_electron_mw', scan=scan):.4e} MW/m$^3$\n" + f"Power per unit volume transferred to ions: {mfile.get('f_pden_alpha_ions_mw', scan=scan):.4e} MW/m$^3$" ) axis.text( @@ -12269,11 +12277,12 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): # ================================================= textstr_neutron = ( - f"Total power: {mfile.get('p_neutron_total_mw', scan=scan):,.2f} MW\n" - f"Plasma power: {mfile.get('p_plasma_neutron_mw', scan=scan):,.2f} MW\n" - f"Beam power: {mfile.get('p_beam_neutron_mw', scan=scan):,.2f} MW\n\n" - f"Total power density: {mfile.get('pden_neutron_total_mw', scan=scan):,.4e} MW/m3\n" - f"Plasma power density: {mfile.get('pden_plasma_neutron_mw', scan=scan):,.4e} MW/m3\n" + f"Total power: {mfile.get('p_neutron_total_mw', scan=scan):,.4f} MW\n" + f"Plasma power: {mfile.get('p_plasma_neutron_mw', scan=scan):,.4f} MW\n" + f"Beam power: {mfile.get('p_beam_neutron_mw', scan=scan):,.4f} MW\n\n" + f"Total power density: {mfile.get('pden_neutron_total_mw', scan=scan):,.4e} MW/m$^3$\n" + f"Plasma power density: {mfile.get('pden_plasma_neutron_mw', scan=scan):,.4e} MW/m$^3$\n\n" + f"Neutron production rate: {mfile.get('fusrat_neutron_production_total', scan=scan):.4e} particles/s" ) axis.text( @@ -12293,8 +12302,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: MFile, scan: int): ) axis.text( - 0.25, - 0.2, + 0.26, + 0.21, "$n$", fontsize=20, verticalalignment="top", @@ -16474,6 +16483,325 @@ def plot_pf_dimensions( axis.set_aspect("equal", adjustable="box") +def plot_tritium_flow_contour(axis: plt.Axes, mfile: MFile, scan: int): + """Plot contour of tritium flow rate vs recycling and fuelling rate.""" + recycling_range = np.linspace(0.01, 0.99, 20) + fuelling_range = np.linspace(0.01, 1.0, 20) + tritium_flow = np.zeros((len(recycling_range), len(fuelling_range))) + + for i, recycling in enumerate(recycling_range): + for j, fuelling in enumerate(fuelling_range): + plasma_fuelling = PlasmaFuelling() + tritium_flow[i, j] = plasma_fuelling.calculate_plasma_tritium_flow_rate( + f_molflow_plasma_fuelling_tritium=mfile.get( + "f_molflow_plasma_fuelling_tritium", scan=scan + ), + eta_plasma_fuelling=fuelling, + molflow_plasma_fuelling_vv_injected=mfile.get( + "molflow_plasma_fuelling_vv_injected", scan=scan + ), + molflow_beam_tritium_vv_injected=mfile.get( + "molflow_beam_tritium_vv_injected", scan=scan + ), + fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), + fusrat_plasma_dd_triton=mfile.get("fusrat_plasma_dd_triton", scan=scan), + t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), + f_plasma_particles_lcfs_recycled=recycling, + nd_plasma_fuel_ions_vol_avg=mfile.get( + "nd_plasma_fuel_ions_vol_avg", scan=scan + ), + vol_plasma=mfile.get("vol_plasma", scan=scan), + f_plasma_fuel_tritium=mfile.get("f_plasma_fuel_tritium", scan=scan), + ) + + contour = axis.contourf( + fuelling_range, + recycling_range, + tritium_flow, + levels=15, + cmap="seismic", + norm=plt.matplotlib.colors.CenteredNorm(vcenter=0), + ) + + axis.contour( + fuelling_range, + recycling_range, + tritium_flow, + levels=[0], + colors="black", + linewidths=2, + ) + + # Plot star for mfile values + recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) + fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan) + axis.plot( + fuelling_mfile, + recycling_mfile, + marker="*", + markersize=15, + color="yellow", + markeredgecolor="black", + markeredgewidth=1.5, + ) + + axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") + axis.set_ylabel("Recycling Fraction [$R$]") + axis.set_title("Plasma Tritium Flow Rate (particles/s)", pad=20) + axis.minorticks_on() + axis.grid(True, which="major", linestyle="-", alpha=0.7) + axis.grid(True, which="minor", linestyle=":", alpha=0.4) + cbar = plt.colorbar(contour, ax=axis, label="Tritium Flow Rate") + cbar.ax.axhline(y=0, color="black", linewidth=2) + + +def plot_deuterium_flow_contour(axis: plt.Axes, mfile: MFile, scan: int): + """Plot contour of deuterium flow rate vs recycling and fuelling rate.""" + recycling_range = np.linspace(0.01, 0.99, 20) + fuelling_range = np.linspace(0.01, 1.0, 20) + deuterium_flow = np.zeros((len(recycling_range), len(fuelling_range))) + + for i, recycling in enumerate(recycling_range): + for j, fuelling in enumerate(fuelling_range): + deuterium_flow[i, j] = PlasmaFuelling().calculate_plasma_deuterium_flow_rate( + f_molflow_plasma_fuelling_deuterium=mfile.get( + "f_molflow_plasma_fuelling_deuterium", scan=scan + ), + eta_plasma_fuelling=fuelling, + molflow_plasma_fuelling_vv_injected=mfile.get( + "molflow_plasma_fuelling_vv_injected", scan=scan + ), + molflow_beam_deuterium_vv_injected=mfile.get( + "molflow_beam_deuterium_vv_injected", scan=scan + ), + fusrat_dt_total=mfile.get("fusrat_dt_total", scan=scan), + fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan), + fusrat_plasma_dd_total=mfile.get("fusrat_plasma_dd_total", scan=scan), + t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), + f_plasma_particles_lcfs_recycled=recycling, + nd_plasma_fuel_ions_vol_avg=mfile.get( + "nd_plasma_fuel_ions_vol_avg", scan=scan + ), + vol_plasma=mfile.get("vol_plasma", scan=scan), + f_plasma_fuel_deuterium=mfile.get("f_plasma_fuel_deuterium", scan=scan), + ) + + contour = axis.contourf( + fuelling_range, + recycling_range, + deuterium_flow, + levels=15, + cmap="seismic", + norm=plt.matplotlib.colors.CenteredNorm(vcenter=0), + ) + axis.contour( + fuelling_range, + recycling_range, + deuterium_flow, + levels=[0], + colors="black", + linewidths=2, + ) + + # Plot star for mfile values + recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) + fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan) + axis.plot( + fuelling_mfile, + recycling_mfile, + marker="*", + markersize=15, + color="yellow", + markeredgecolor="black", + markeredgewidth=1.5, + ) + + axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") + axis.set_title("Plasma Deuterium Flow Rate (particles/s)", pad=20) + axis.minorticks_on() + axis.grid(True, which="major", linestyle="-", alpha=0.7) + axis.grid(True, which="minor", linestyle=":", alpha=0.4) + cbar = plt.colorbar(contour, ax=axis, label="Deuterium Flow Rate") + cbar.ax.axhline(y=0, color="black", linewidth=2) + + +def plot_alpha_flow_contour(axis: plt.Axes, mfile: MFile, scan: int): + """Plot contour of alpha particle flow rate vs recycling and fuelling rate.""" + fusion_dt_range = np.linspace(1e19, 5e21, 20) + f_t_alpha_energy_confinement_range = np.linspace(2.0, 10.0, 20) + alpha_flow = np.zeros(( + len(fusion_dt_range), + len(f_t_alpha_energy_confinement_range), + )) + + for i, fusion_dt in enumerate(fusion_dt_range): + for j, f_t_alpha_energy_confinement in enumerate( + f_t_alpha_energy_confinement_range + ): + alpha_flow[i, j] = ( + PlasmaFuelling().calculate_plasma_alphas_thermal_flow_rate( + fusrat_dt_total=fusion_dt, + fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan), + t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), + nd_plasma_alphas_thermal_vol_avg=mfile.get( + "nd_plasma_alphas_thermal_vol_avg", scan=scan + ), + vol_plasma=mfile.get("vol_plasma", scan=scan), + f_t_alpha_energy_confinement=f_t_alpha_energy_confinement, + ) + ) + + contour = axis.contourf( + f_t_alpha_energy_confinement_range, + fusion_dt_range, + alpha_flow, + levels=15, + cmap="seismic", + norm=plt.matplotlib.colors.CenteredNorm(vcenter=0), + ) + axis.contour( + f_t_alpha_energy_confinement_range, + fusion_dt_range, + alpha_flow, + levels=[0], + colors="black", + linewidths=2, + ) + + # Plot star for mfile values + fusion_dt_mfile = mfile.get("fusrat_dt_total", scan=scan) + f_t_alpha_mfile = mfile.get("f_t_alpha_energy_confinement", scan=scan) + axis.plot( + f_t_alpha_mfile, + fusion_dt_mfile, + marker="*", + markersize=15, + color="yellow", + markeredgecolor="black", + markeredgewidth=1.5, + ) + + axis.set_xlabel( + "Alpha to Energy Confinement Time Ratio ($f_{\\alpha, \\text{energy confinement}}$)" + ) + axis.set_ylabel("Fusion DT Rate [$\\text{particles/s}$]") + axis.set_title("Plasma Alpha Particle Flow Rate (particles/s)", pad=20) + axis.minorticks_on() + axis.grid(True, which="major", linestyle="-", alpha=0.7) + axis.grid(True, which="minor", linestyle=":", alpha=0.4) + cbar = plt.colorbar(contour, ax=axis, label="Alpha Particle Flow Rate") + cbar.ax.axhline(y=0, color="black", linewidth=2) + + +def plot_helium3_flow_contour(axis: plt.Axes, mfile: MFile, scan: int): + """Plot contour of helium-3 flow rate vs recycling and fuelling rate.""" + recycling_range = np.linspace(0.01, 0.99, 20) + fuelling_range = np.linspace(0.01, 1.0, 20) + helium3_flow = np.zeros((len(recycling_range), len(fuelling_range))) + + for i, recycling in enumerate(recycling_range): + for j, fuelling in enumerate(fuelling_range): + helium3_flow[i, j] = PlasmaFuelling().calculate_plasma_helium3_flow_rate( + f_molflow_plasma_fuelling_helium3=mfile.get( + "f_molflow_plasma_fuelling_helium3", scan=scan + ), + eta_plasma_fuelling=fuelling, + molflow_plasma_fuelling_vv_injected=mfile.get( + "molflow_plasma_fuelling_vv_injected", scan=scan + ), + fusrat_plasma_dhe3=mfile.get("fusrat_plasma_dhe3", scan=scan), + fusrat_plasma_dd_helion=mfile.get("fusrat_plasma_dd_helion", scan=scan), + t_energy_confinement=mfile.get("t_energy_confinement", scan=scan), + f_plasma_particles_lcfs_recycled=recycling, + nd_plasma_fuel_ions_vol_avg=mfile.get( + "nd_plasma_fuel_ions_vol_avg", scan=scan + ), + vol_plasma=mfile.get("vol_plasma", scan=scan), + f_plasma_fuel_helium3=mfile.get("f_plasma_fuel_helium3", scan=scan), + ) + + contour = axis.contourf( + fuelling_range, + recycling_range, + helium3_flow, + levels=15, + cmap="seismic", + norm=plt.matplotlib.colors.CenteredNorm(vcenter=0), + ) + axis.contour( + fuelling_range, + recycling_range, + helium3_flow, + levels=[0], + colors="black", + linewidths=2, + ) + + # Plot star for mfile values + recycling_mfile = mfile.get("f_plasma_particles_lcfs_recycled", scan=scan) + fuelling_mfile = mfile.get("eta_plasma_fuelling", scan=scan) + axis.plot( + fuelling_mfile, + recycling_mfile, + marker="*", + markersize=15, + color="yellow", + markeredgecolor="black", + ) + axis.set_xlabel("Fuelling Rate Efficiency ($\\eta_{\\text{fuelling}}$)") + cbar = plt.colorbar(contour, ax=axis, label="Helium-3 Flow Rate") + cbar.ax.axhline(y=0, color="black", linewidth=2) + axis.set_title("Plasma Helium-3 Flow Rate (particles/s)", pad=20) + axis.minorticks_on() + axis.grid(True, which="major", linestyle="-", alpha=0.7) + axis.grid(True, which="minor", linestyle=":", alpha=0.4) + + +def plot_fuelling_info(fig: plt.Figure, mfile: MFile, scan: int): + """Plot fuelling information.""" + msg = ( + f"$\\mathbf{{Plasma \\ Fuelling \\ Information:}}$\n\n" + f"Fuelling rate into vessel:" + f"{mfile.get('molflow_plasma_fuelling_vv_injected', scan=scan):.4e} particles/s\n" + f"Fuelling rate into vessel: " + f"{mfile.get('molflow_plasma_fuelling_vv_injected_moles', scan=scan):.4e} moles/s\n" + f"Additional deuterium beam fuelling: " + f"{mfile.get('molflow_beam_deuterium_vv_injected', scan=scan):.4e} particles/s\n" + f"Additional tritium beam fuelling: " + f"{mfile.get('molflow_beam_tritium_vv_injected', scan=scan):.4e} particles/s\n\n" + f"Fuelling loss that doesn't enter the plasma: " + f"{mfile.get('molflow_plasma_fuelling_loss', scan=scan):.4e} particles/s\n" + f"Fuelling loss that doesn't enter the plasma: " + f"{mfile.get('molflow_plasma_fuelling_loss_moles', scan=scan):.4e} moles/s\n" + f"Fuelling Rate Efficiency ($\\eta_{{\\text{{fuelling}}}}$): " + f"{mfile.get('eta_plasma_fuelling', scan=scan):.4f}\n" + f"Recycling Fraction ($R$): " + f"{mfile.get('f_plasma_particles_lcfs_recycled', scan=scan):.4f}\n\n" + f"Fraction of Tritium Fuelling: " + f"{mfile.get('f_molflow_plasma_fuelling_tritium', scan=scan):.4f}\n" + f"Fraction of Deuterium Fuelling: " + f"{mfile.get('f_molflow_plasma_fuelling_deuterium', scan=scan):.4f}\n" + f"Fraction of 3-Helium Fuelling: " + f"{mfile.get('f_molflow_plasma_fuelling_helium3', scan=scan):.4f}\n\n" + f"Total Fuel Burnup Fraction: " + f"{mfile.get('f_plasma_fuel_burnup', scan=scan):.4f}\n" + f"Tritium Burnup Fraction: " + f"{mfile.get('f_plasma_tritium_burnup', scan=scan):.4f}\n" + f"Deuterium Burnup Fraction: " + f"{mfile.get('f_plasma_deuterium_burnup', scan=scan):.4f}" + ) + fig.text( + 0.75, + 0.25, + msg, + ha="center", + va="center", + transform=fig.transFigure, + fontsize=9, + bbox={"boxstyle": "round", "facecolor": "wheat", "alpha": 1.0}, + ) + + def main_plot( m_file: MFile, scan: int, @@ -16647,6 +16975,21 @@ def _add_page(name: str | None = None): if m_file.get("i_plasma_shape", scan=scan) == PlasmaShapeModelType.SAUTER: plot_fusion_rate_contours(pages["rx_1_2"], pages["rx_3_4"], m_file, scan) + plot_tritium_flow_contour( + axis=_add_page("fuelling").add_subplot(231), mfile=m_file, scan=scan + ) + plot_deuterium_flow_contour( + axis=pages["fuelling"].add_subplot(232), mfile=m_file, scan=scan + ) + plot_helium3_flow_contour( + axis=pages["fuelling"].add_subplot(233), mfile=m_file, scan=scan + ) + plot_alpha_flow_contour( + axis=pages["fuelling"].add_subplot(223), mfile=m_file, scan=scan + ) + + plot_fuelling_info(pages["fuelling"], m_file, scan) + if i_shape != PlasmaShapeModelType.SAUTER: msg = ( "Fusion-rate contour plots require a closed (Sauter) plasma boundary " diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 72c11014fd..6826204d2c 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -14,6 +14,7 @@ from process.data_structure.build_variables import TFCSRadialConfiguration from process.models.physics.density_limit import DensityLimitModel from process.models.physics.exhaust import PlasmaExhaust +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.physics import ( BetaComponentLimits, PlasmaBeta, @@ -1955,6 +1956,155 @@ def constraint_equation_92(constraint_registration, data): ) +@ConstraintManager.register_constraint(93, "particles/s", "=") +def constraint_equation_93(constraint_registration, data): + """ + Tritium particle balance consistency equation. + + This constraint ensures that the tritium particle balance is maintained in the + plasma. + + The numerator represents the total tritium sources, including fuelling and fusion + reactions, while the denominator represents the total tritium sinks, which are + related to the plasma volume, fuel ion density, and confinement time. + + The equation ensures that the rate of tritium production equals the rate of tritium + loss, maintaining a steady-state condition for tritium in the plasma. + """ + source = PlasmaFuelling.calculate_plasma_tritium_source_rate( + f_molflow_plasma_fuelling_tritium=data.physics.f_molflow_plasma_fuelling_tritium, + eta_plasma_fuelling=data.physics.eta_plasma_fuelling, + molflow_plasma_fuelling_vv_injected=data.physics.molflow_plasma_fuelling_vv_injected, + fusrat_plasma_dd_triton=data.physics.fusrat_plasma_dd_triton, + molflow_beam_tritium_vv_injected=data.physics.molflow_beam_tritium_vv_injected, + ) + sink = PlasmaFuelling.calculate_plasma_tritium_loss_rate( + fusrat_dt_total=data.physics.fusrat_dt_total, + t_energy_confinement=data.physics.t_energy_confinement, + f_plasma_particles_lcfs_recycled=data.physics.f_plasma_particles_lcfs_recycled, + nd_plasma_fuel_ions_vol_avg=data.physics.nd_plasma_fuel_ions_vol_avg, + vol_plasma=data.physics.vol_plasma, + f_plasma_fuel_tritium=data.physics.f_plasma_fuel_tritium, + ) + + return eq(source, -sink, constraint_registration) + + +@ConstraintManager.register_constraint(94, "particles/s", "=") +def constraint_equation_94(constraint_registration, data): + """ + Deuterium particle balance consistency equation. + + This constraint ensures that the deuterium particle balance is maintained in the + plasma. + + The numerator represents the total deuterium sources, including fuelling and fusion + reactions, while the denominator represents the total deuterium sinks, which are + related to the plasma volume, fuel ion density, and confinement time. + + The equation ensures that the rate of deuterium production equals the rate of + deuterium loss, maintaining a steady-state condition for deuterium in the plasma. + """ + source = PlasmaFuelling.calculate_plasma_deuterium_source_rate( + f_molflow_plasma_fuelling_deuterium=data.physics.f_molflow_plasma_fuelling_deuterium, + eta_plasma_fuelling=data.physics.eta_plasma_fuelling, + molflow_plasma_fuelling_vv_injected=data.physics.molflow_plasma_fuelling_vv_injected, + molflow_beam_deuterium_vv_injected=data.physics.molflow_beam_deuterium_vv_injected, + ) + sink = PlasmaFuelling.calculate_plasma_deuterium_loss_rate( + fusrat_dt_total=data.physics.fusrat_dt_total, + fusrat_plasma_dd_total=data.physics.fusrat_plasma_dd_total, + fusrat_plasma_dhe3=data.physics.fusrat_plasma_dhe3, + t_energy_confinement=data.physics.t_energy_confinement, + f_plasma_particles_lcfs_recycled=data.physics.f_plasma_particles_lcfs_recycled, + nd_plasma_fuel_ions_vol_avg=data.physics.nd_plasma_fuel_ions_vol_avg, + vol_plasma=data.physics.vol_plasma, + f_plasma_fuel_deuterium=data.physics.f_plasma_fuel_deuterium, + ) + + return eq(source, -sink, constraint_registration) + + +@ConstraintManager.register_constraint(95, "particles/s", "=") +def constraint_equation_95(constraint_registration, data): + """ + Helium-3 particle balance consistency equation. + + This constraint ensures that the helium-3 particle balance is maintained in the + plasma. + + The numerator represents the total helium-3 sources, including fuelling and fusion + reactions, while the denominator represents the total helium-3 sinks, which are + related to the plasma volume, fuel ion density, and confinement time. + + The equation ensures that the rate of helium-3 production equals the rate of helium-3 + loss, maintaining a steady-state condition for helium-3 in the plasma. + """ + source = PlasmaFuelling.calculate_plasma_helium3_source_rate( + f_molflow_plasma_fuelling_helium3=data.physics.f_molflow_plasma_fuelling_helium3, + eta_plasma_fuelling=data.physics.eta_plasma_fuelling, + molflow_plasma_fuelling_vv_injected=data.physics.molflow_plasma_fuelling_vv_injected, + fusrat_plasma_dd_helion=data.physics.fusrat_plasma_dd_helion, + ) + sink = PlasmaFuelling.calculate_plasma_helium3_loss_rate( + fusrat_plasma_dhe3=data.physics.fusrat_plasma_dhe3, + t_energy_confinement=data.physics.t_energy_confinement, + f_plasma_particles_lcfs_recycled=data.physics.f_plasma_particles_lcfs_recycled, + nd_plasma_fuel_ions_vol_avg=data.physics.nd_plasma_fuel_ions_vol_avg, + vol_plasma=data.physics.vol_plasma, + f_plasma_fuel_helium3=data.physics.f_plasma_fuel_helium3, + ) + + return eq(source, -sink, constraint_registration) + + +@ConstraintManager.register_constraint(96, "particles/s", "=") +def constraint_equation_96(constraint_registration, data): + """ + Thermal alpha particle balance consistency equation. + + This constraint ensures that the thermal alpha particle balance is maintained in the + plasma. + + The numerator represents the total thermal alpha particle sources, + while the denominator represents the total thermal alpha particle sinks, which are + related to the plasma volume, fuel ion density, and confinement time. + + The equation ensures that the rate of thermal alpha particle production equals the + rate of thermal alpha particle loss, maintaining a steady-state condition for + thermal alpha particles in the plasma. + """ + source = PlasmaFuelling.calculate_plasma_alphas_thermal_source_rate( + fusrat_dt_total=data.physics.fusrat_dt_total, + fusrat_plasma_dhe3=data.physics.fusrat_plasma_dhe3, + ) + sink = PlasmaFuelling.calculate_plasma_alphas_thermal_loss_rate( + nd_plasma_alphas_thermal_vol_avg=data.physics.nd_plasma_alphas_thermal_vol_avg, + vol_plasma=data.physics.vol_plasma, + t_energy_confinement=data.physics.t_energy_confinement, + f_t_alpha_energy_confinement=data.physics.f_t_alpha_energy_confinement, + ) + + return eq(source, -sink, constraint_registration) + + +@ConstraintManager.register_constraint(97, "", "=") +def constraint_equation_97(constraint_registration, data): + """Equation for checking the fuelling composition is consistent. + + f_molflow_plasma_fuelling_deuterium: fraction of deuterium ions + f_molflow_plasma_fuelling_tritium: fraction of tritium ions + f_molflow_plasma_fuelling_helium3: fraction of helium-3 ions + """ + return eq( + data.physics.f_molflow_plasma_fuelling_deuterium + + data.physics.f_molflow_plasma_fuelling_tritium + + data.physics.f_molflow_plasma_fuelling_helium3, + 1.0, + constraint_registration, + ) + + def constraint_eqns(m: int, ieqn: int, data: DataStructure): """Evaluates the constraints given the current state of PROCESS. diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index d1df7c6feb..65c25b7f4b 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -239,6 +239,32 @@ class IterationVariable: 175: IterationVariable("kappa", "physics", 0.00, 10.00), 176: IterationVariable("f_st_coil_aspect", "stellarator", 0.70, 1.30), 177: IterationVariable("f_a_tf_turn_cable_space_extra_void", "tfcoil", 0.01, 1.0), + 178: IterationVariable("f_plasma_particles_lcfs_recycled", "physics", 0.01, 1.0), + 179: IterationVariable("eta_plasma_fuelling", "physics", 0.01, 1.0), + 180: IterationVariable( + "molflow_plasma_fuelling_vv_injected", + "physics", + 1e18, + 1e22, + ), + 181: IterationVariable( + "f_molflow_plasma_fuelling_deuterium", + "physics", + 0.0, + 1.0, + ), + 182: IterationVariable( + "f_molflow_plasma_fuelling_tritium", + "physics", + 0.0, + 1.0, + ), + 183: IterationVariable( + "f_molflow_plasma_fuelling_helium3", + "physics", + 0.0, + 1.0, + ), } diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 43b0b78314..1c306b767b 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -292,6 +292,11 @@ class NumericsData: "CS achievable stress load cycles lower limit ", "ECRH ignitability ", # Stellarator constraint "Fuel composition consistency ", + "Tritium particle balance consistency ", + "Deuterium particle balance consistency ", + "Helium-3 particle balance consistency ", + "Alpha particle balance consistency ", + "Fuelling composition consistency ", ] ) """Labels describing constraint equations (corresponding itvs)