diff --git a/documentation/source/physics-models/plasma_fuelling.md b/documentation/source/physics-models/plasma_fuelling.md new file mode 100644 index 000000000..86dfd90ea --- /dev/null +++ b/documentation/source/physics-models/plasma_fuelling.md @@ -0,0 +1,154 @@ +# Plasma Fuelling | `PlasmaFuelling()` + +The control of fuelling is governed by 4 key particle flux equations for each of the primary fuel species and the helium ash, $\alpha$. + +$$ +\frac{dn_{\text{T}}}{dt} = f_{\text{fuelling,T}}\eta_{\text{fuelling}}\Gamma_{\text{fuel}} + \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}} -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}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*} +$$ + +$$ +\frac{dn_{\alpha}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha}}{\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}}{dt} = 0 +$$ + +Here $\eta_{\text{fuelling}}$ is the fuelling efficiecny which represents the method of injecting fuel into the plasma. Gas puffing on the low field side is probably around 0.01-0.1, supersonic gas is 0.1 and 0.2 and using pellets can get you close to unity with 0.5-0.9. $\Gamma_{\text{fuelling}}$ is the fuel injection rate into the vacuum vessel, so $\eta_{\text{fuelling}} \Gamma_{\text{fuelling}}$ together presents the fraction of injected fuel that actually makes it into the plasma core to fuse. + + +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$). + +!!! 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]. + + + + +## METIS Alpha Confinement + +$$ +\tau_{\alpha} = f_{\alpha}\tau_{\text{E}}\frac{R}{1-R}\tau_{\text{ne}} +$$ + +This is the model currently in METIS[^2] and is found in [^3] + + +-------------- + +## 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{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}} -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}} - \frac{N_{\text{T}}}{\tau_{\text{3He}}^*} +$$ + +--------------- + +## Alpha Particle Flow Rate | `calculate_plasma_alphas_flow_rate()` + +$$ +\frac{dn_{\alpha}}{dt} = \Gamma_{\text{D+3He}} + \Gamma_{\text{D+T}} - \frac{N_{\alpha}}{\tau_{\alpha}^*} +$$ + +----------------- + +## Key Constraints + +### 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 consistency model** + +---------------- + +### 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 consistency model** + +----------------- + +### Helium-3 Flow Consistency + +### Alpha Particle Flow Consistency + +This constraint can be activated by stating `icc = 95` in the input file. + +This constraint ensures that the change in alpha particles as a function of time is zero. It ensures the output of `calculate_plasma_alphas_flow_rate(()` is zero + +**It is recommended to have this constraint on as it is a plasma consistency model** + +------------------ + +### Fuelling Proportion Consistency + +This constraint can be activated by stating `icc = 96` 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 6ccdc0a65..bd75f3996 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -52,6 +52,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 100fe4d65..abf35cfc0 100644 --- a/process/core/constants.py +++ b/process/core/constants.py @@ -7,6 +7,13 @@ 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 7cb7b7657..94dfb526d 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -1208,6 +1208,24 @@ def __post_init__(self): "molflow_vac_pumps": InputVariable( data_structure.vacuum_variables, float, range=(0.0, 1e30) ), + "f_plasma_particles_lcfs_recycled": InputVariable( + data_structure.physics_variables, float, range=(0.0, 1.0) + ), + "eta_plasma_fuelling": InputVariable( + data_structure.physics_variables, float, range=(0.0, 1.0) + ), + "molflow_plasma_fuelling_vv_injected": InputVariable( + data_structure.physics_variables, float, range=(1e18, 1e24) + ), + "f_molflow_plasma_fuelling_deuterium": InputVariable( + data_structure.physics_variables, float, range=(0.0, 1.0) + ), + "f_molflow_plasma_fuelling_tritium": InputVariable( + data_structure.physics_variables, float, range=(0.0, 1.0) + ), + "f_molflow_plasma_fuelling_helium3": InputVariable( + data_structure.physics_variables, float, range=(0.0, 1.0) + ), "pflux_plant_floor_electric": InputVariable( data_structure.heat_transport_variables, float, range=(0.0, 1000.0) ), @@ -1464,9 +1482,6 @@ def __post_init__(self): "t_plasma_energy_confinement_max": InputVariable( data_structure.physics_variables, float, range=(0.1, 100.0) ), - "tauratio": InputVariable( - data_structure.physics_variables, float, range=(0.1, 100.0) - ), "n_beam_decay_lengths_core_required": InputVariable( data_structure.current_drive_variables, float, range=(0.0, 10.0) ), diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index 9165eb484..4916a3569 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -67,6 +67,7 @@ ElectronBernstein, ElectronCyclotron, ) +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.plasma_current import PlasmaCurrent, PlasmaCurrentModel @@ -2978,7 +2979,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" ) @@ -3002,9 +3003,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( @@ -11299,8 +11300,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.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( @@ -11322,9 +11323,11 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.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( @@ -11344,7 +11347,7 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): ) axis.text( - 0.24, + 0.285, 0.8, "$\\text{D - T}$", fontsize=20, @@ -11355,13 +11358,16 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.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", @@ -11376,8 +11382,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): ) axis.text( - 0.22, - 0.685, + 0.31, + 0.69, "$\\text{D - D}$", fontsize=20, verticalalignment="top", @@ -11386,11 +11392,14 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.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, + 0.525, textstr_dhe3, fontsize=9, verticalalignment="bottom", @@ -11405,8 +11414,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): ) axis.text( - 0.21, - 0.59, + 0.285, + 0.56, "$\\text{D - 3He}$", fontsize=20, verticalalignment="top", @@ -11416,15 +11425,15 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.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$/sec\n" + f"Rate density, plasma: {mfile.get('fusden_plasma_alpha', scan=scan):.4e} particles/m$^3$/sec\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( @@ -11455,11 +11464,12 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.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( @@ -11479,8 +11489,8 @@ def plot_fusion_rate_profiles(axis: plt.Axes, fig, mfile: mf.MFile, scan: int): ) axis.text( - 0.25, - 0.2, + 0.26, + 0.21, "$n$", fontsize=20, verticalalignment="top", @@ -13706,16 +13716,33 @@ def main_plot( figs[9].text(0.5, 0.5, msg, ha="center", va="center", wrap=True, fontsize=12) figs[10].text(0.5, 0.5, msg, ha="center", va="center", wrap=True, fontsize=12) - plot_plasma_pressure_profiles(figs[11].add_subplot(222), m_file, scan) - plot_plasma_pressure_gradient_profiles(figs[11].add_subplot(224), m_file, scan) + PlasmaFuelling().plot_tritium_flow_contour( + axis=figs[11].add_subplot(231), mfile=m_file, scan=scan + ) + PlasmaFuelling().plot_deuterium_flow_contour( + axis=figs[11].add_subplot(232), mfile=m_file, scan=scan + ) + PlasmaFuelling().plot_helium3_flow_contour( + axis=figs[11].add_subplot(233), mfile=m_file, scan=scan + ) + PlasmaFuelling().plot_alpha_flow_contour( + axis=figs[11].add_subplot(223), mfile=m_file, scan=scan + ) + + PlasmaFuelling().plot_fuelling_info(figs[11], m_file, scan) + + figs[11].subplots_adjust(wspace=0.3, hspace=0.35) + + plot_plasma_pressure_profiles(figs[12].add_subplot(222), m_file, scan) + plot_plasma_pressure_gradient_profiles(figs[12].add_subplot(224), m_file, scan) # Currently only works with Sauter geometry as plasma has a closed surface if i_shape == 1: plot_plasma_poloidal_pressure_contours( - figs[11].add_subplot(121, aspect="equal"), m_file, scan + figs[12].add_subplot(121, aspect="equal"), m_file, scan ) else: - ax = figs[11].add_subplot(131, aspect="equal") + ax = figs[12].add_subplot(131, aspect="equal") msg = ( "Plasma poloidal pressure contours require a closed (Sauter) plasma boundary " "(i_plasma_shape == 1). " @@ -13734,30 +13761,30 @@ def main_plot( ) ax.axis("off") - plot_magnetic_fields_in_plasma(figs[12].add_subplot(122), m_file, scan) - plot_beta_profiles(figs[12].add_subplot(221), m_file, scan) + plot_magnetic_fields_in_plasma(figs[13].add_subplot(122), m_file, scan) + plot_beta_profiles(figs[13].add_subplot(221), m_file, scan) - plot_ebw_ecrh_coupling_graph(figs[13].add_subplot(111), m_file, scan) + plot_ebw_ecrh_coupling_graph(figs[14].add_subplot(111), m_file, scan) - plot_bootstrap_comparison(figs[14].add_subplot(221), m_file, scan) - PlasmaCurrent.plot_plasma_current_comparison(figs[14].add_subplot(224), m_file, scan) - plot_h_threshold_comparison(figs[15].add_subplot(224), m_file, scan) - plot_density_limit_comparison(figs[15].add_subplot(221), m_file, scan) - plot_confinement_time_comparison(figs[16].add_subplot(224), m_file, scan) + plot_bootstrap_comparison(figs[15].add_subplot(221), m_file, scan) + PlasmaCurrent.plot_plasma_current_comparison(figs[15].add_subplot(224), m_file, scan) + plot_h_threshold_comparison(figs[16].add_subplot(224), m_file, scan) + plot_density_limit_comparison(figs[16].add_subplot(221), m_file, scan) + plot_confinement_time_comparison(figs[17].add_subplot(224), m_file, scan) - plot_debye_length_profile(figs[17].add_subplot(232), m_file, scan) - plot_velocity_profile(figs[17].add_subplot(233), m_file, scan) - plot_plasma_coloumb_logarithms(figs[17].add_subplot(231), m_file, scan) + plot_debye_length_profile(figs[18].add_subplot(232), m_file, scan) + plot_velocity_profile(figs[18].add_subplot(233), m_file, scan) + plot_plasma_coloumb_logarithms(figs[18].add_subplot(231), m_file, scan) - plot_electron_frequency_profile(figs[17].add_subplot(212), m_file, scan) + plot_electron_frequency_profile(figs[18].add_subplot(212), m_file, scan) - plot_ion_frequency_profile(figs[18].add_subplot(311), m_file, scan) + plot_ion_frequency_profile(figs[19].add_subplot(311), m_file, scan) - plot_larmor_radius_profile(figs[18].add_subplot(313), m_file, scan) + plot_larmor_radius_profile(figs[19].add_subplot(313), m_file, scan) # Plot poloidal cross-section poloidal_cross_section( - figs[19].add_subplot(121, aspect="equal"), + figs[20].add_subplot(121, aspect="equal"), m_file, scan, demo_ranges, @@ -13767,7 +13794,7 @@ def main_plot( # Plot toroidal cross-section toroidal_cross_section( - figs[19].add_subplot(122, aspect="equal"), + figs[20].add_subplot(122, aspect="equal"), m_file, scan, demo_ranges, @@ -13775,27 +13802,27 @@ def main_plot( ) # Plot color key - ax17 = figs[19].add_subplot(222) + ax17 = figs[20].add_subplot(222) ax17.set_position([0.5, 0.5, 0.5, 0.5]) color_key(ax17, m_file, scan, colour_scheme) plot_full_machine_poloidal_cross_section( - figs[20].add_subplot(111, aspect="equal"), + figs[21].add_subplot(111, aspect="equal"), m_file, scan, radial_build, colour_scheme, ) - ax18 = figs[21].add_subplot(211) + ax18 = figs[22].add_subplot(211) ax18.set_position([0.1, 0.33, 0.8, 0.6]) plot_radial_build(ax18, m_file, colour_scheme) # Make each axes smaller vertically to leave room for the legend - ax185 = figs[22].add_subplot(211) + ax185 = figs[23].add_subplot(211) ax185.set_position([0.1, 0.61, 0.8, 0.32]) - ax18b = figs[22].add_subplot(212) + ax18b = figs[23].add_subplot(212) ax18b.set_position([0.1, 0.13, 0.8, 0.32]) plot_upper_vertical_build(ax185, m_file, colour_scheme) plot_lower_vertical_build(ax18b, m_file, colour_scheme) @@ -13803,62 +13830,62 @@ def main_plot( # Can only plot WP and turn structure if superconducting coil at the moment if m_file.get("i_tf_sup", scan=scan) == 1: # TF coil with WP - ax19 = figs[23].add_subplot(221, aspect="equal") + ax19 = figs[24].add_subplot(221, aspect="equal") ax19.set_position([ 0.025, 0.45, 0.5, 0.5, ]) # Half height, a bit wider, top left - plot_superconducting_tf_wp(ax19, m_file, scan, figs[23]) + plot_superconducting_tf_wp(ax19, m_file, scan, figs[24]) # TF coil turn structure - ax20 = figs[24].add_subplot(325, aspect="equal") + ax20 = figs[25].add_subplot(325, aspect="equal") ax20.set_position([0.025, 0.5, 0.4, 0.4]) - plot_tf_cable_in_conduit_turn(ax20, figs[24], m_file, scan) - plot_205 = figs[24].add_subplot(223, aspect="equal") + plot_tf_cable_in_conduit_turn(ax20, figs[25], m_file, scan) + plot_205 = figs[25].add_subplot(223, aspect="equal") plot_205.set_position([0.075, 0.1, 0.3, 0.3]) - plot_cable_in_conduit_cable(plot_205, figs[24], m_file, scan) + plot_cable_in_conduit_cable(plot_205, figs[25], m_file, scan) else: - ax19 = figs[23].add_subplot(211, aspect="equal") + ax19 = figs[24].add_subplot(211, aspect="equal") ax19.set_position([0.06, 0.55, 0.675, 0.4]) - plot_resistive_tf_wp(ax19, m_file, scan, figs[23]) - plot_resistive_tf_info(ax19, m_file, scan, figs[23]) + plot_resistive_tf_wp(ax19, m_file, scan, figs[24]) + plot_resistive_tf_info(ax19, m_file, scan, figs[24]) plot_tf_coil_structure( - figs[25].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme + figs[26].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme ) - plot_plasma_outboard_toroidal_ripple_map(figs[26], m_file, scan) + plot_plasma_outboard_toroidal_ripple_map(figs[27], m_file, scan) - plot_tf_stress(figs[27].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) + plot_tf_stress(figs[28].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) - plot_current_profiles_over_time(figs[28].add_subplot(111), m_file, scan) + plot_current_profiles_over_time(figs[29].add_subplot(111), m_file, scan) plot_cs_coil_structure( - figs[29].add_subplot(121, aspect="equal"), figs[29], m_file, scan + figs[30].add_subplot(121, aspect="equal"), figs[30], m_file, scan ) plot_cs_turn_structure( - figs[29].add_subplot(224, aspect="equal"), figs[29], m_file, scan + figs[30].add_subplot(224, aspect="equal"), figs[30], m_file, scan ) plot_first_wall_top_down_cross_section( - figs[30].add_subplot(221, aspect="equal"), m_file, scan + figs[31].add_subplot(221, aspect="equal"), m_file, scan ) - plot_first_wall_poloidal_cross_section(figs[30].add_subplot(122), m_file, scan) - plot_fw_90_deg_pipe_bend(figs[30].add_subplot(337), m_file, scan) + plot_first_wall_poloidal_cross_section(figs[31].add_subplot(122), m_file, scan) + plot_fw_90_deg_pipe_bend(figs[31].add_subplot(337), m_file, scan) - plot_blkt_pipe_bends(figs[31], m_file, scan) - ax_blanket = figs[31].add_subplot(122, aspect="equal") - plot_blkt_structure(ax_blanket, figs[31], m_file, scan, radial_build, colour_scheme) + plot_blkt_pipe_bends(figs[32], m_file, scan) + ax_blanket = figs[32].add_subplot(122, aspect="equal") + plot_blkt_structure(ax_blanket, figs[32], m_file, scan, radial_build, colour_scheme) plot_main_power_flow( - figs[32].add_subplot(111, aspect="equal"), m_file, scan, figs[32] + figs[33].add_subplot(111, aspect="equal"), m_file, scan, figs[33] ) - ax24 = figs[33].add_subplot(111) + ax24 = figs[34].add_subplot(111) # set_position([left, bottom, width, height]) -> height ~ 0.66 => ~2/3 of page height ax24.set_position([0.08, 0.35, 0.84, 0.57]) - plot_system_power_profiles_over_time(ax24, m_file, scan, figs[33]) + plot_system_power_profiles_over_time(ax24, m_file, scan, figs[34]) def create_thickness_builds(m_file, scan: int): @@ -13935,7 +13962,7 @@ def main(args=None): # create main plot # Increase range when adding new page - pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(34)] + pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(35)] # run main_plot main_plot( diff --git a/process/core/scan.py b/process/core/scan.py index b2a7035ad..0ac7c495f 100644 --- a/process/core/scan.py +++ b/process/core/scan.py @@ -615,6 +615,12 @@ def post_optimise(self, ifail: int): f"(ineq_con{numerics.icc[i]:03d})", -constraint.normalised_residual, ) + process_output.ovarre( + constants.MFILE, + f"{name} error", + f"(ineq_err{numerics.icc[i]:03d})", + -constraint.residual, + ) process_output.ovarre( constants.MFILE, f"{name} physical value", diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index 949b4439b..25629ed7b 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -1855,6 +1855,103 @@ def constraint_equation_92(constraint_registration): ) +@ConstraintManager.register_constraint(93, "particles/s", "=") +def constraint_equation_93(constraint_registration): + """ + Tritium particle balance + """ + # pscaling: total transport power per volume (MW/m3) + # NEED TO ADD PROPER FUSION RATES FOR EACH REACTION + + # Numerator shall be the tritium particle balance + numerator = ( + ( + data_structure.physics_variables.eta_plasma_fuelling + * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected + * data_structure.physics_variables.f_molflow_plasma_fuelling_tritium + ) + + data_structure.physics_variables.fusrat_plasma_dd_triton + - data_structure.physics_variables.fusrat_dt_total + ) + denominator = ( + data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg + * data_structure.physics_variables.vol_plasma + * data_structure.physics_variables.f_plasma_fuel_tritium + ) / ( + data_structure.physics_variables.t_energy_confinement + / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) + ) + + return eq(numerator, denominator, constraint_registration) + + +@ConstraintManager.register_constraint(94, "particles/s", "=") +def constraint_equation_94(constraint_registration): + """ + Deuterium particle balance + """ + + numerator = ( + ( + data_structure.physics_variables.eta_plasma_fuelling + * data_structure.physics_variables.molflow_plasma_fuelling_vv_injected + * data_structure.physics_variables.f_molflow_plasma_fuelling_deuterium + ) + - data_structure.physics_variables.fusrat_dt_total + - data_structure.physics_variables.fusrat_plasma_dhe3 + - 2.0 * data_structure.physics_variables.fusrat_plasma_dd_total + ) + denominator = ( + data_structure.physics_variables.nd_plasma_fuel_ions_vol_avg + * data_structure.physics_variables.vol_plasma + * data_structure.physics_variables.f_plasma_fuel_deuterium + ) / ( + data_structure.physics_variables.t_energy_confinement + / (1 - data_structure.physics_variables.f_plasma_particles_lcfs_recycled) + ) + + return eq(numerator, denominator, constraint_registration) + + +@ConstraintManager.register_constraint(95, "particles/s", "=") +def constraint_equation_95(constraint_registration): + """ + Alpha particle balance + """ + + # Alpha particle balance + numerator = ( + data_structure.physics_variables.fusrat_dt_total + + data_structure.physics_variables.fusrat_plasma_dhe3 + ) + denominator = ( + data_structure.physics_variables.nd_plasma_alphas_vol_avg + * data_structure.physics_variables.vol_plasma + ) / ( + data_structure.physics_variables.t_energy_confinement + * data_structure.physics_variables.f_alpha_energy_confinement + ) + + return eq(numerator, denominator, constraint_registration) + + +@ConstraintManager.register_constraint(96, "", "=") +def constraint_equation_96(constraint_registration): + """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_structure.physics_variables.f_molflow_plasma_fuelling_deuterium + + data_structure.physics_variables.f_molflow_plasma_fuelling_tritium + + data_structure.physics_variables.f_molflow_plasma_fuelling_helium3, + 1.0, + constraint_registration, + ) + + def constraint_eqns(m: int, ieqn: int): """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 edae1a409..870f5837c 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -291,6 +291,36 @@ class IterationVariable: 176: IterationVariable( "f_st_coil_aspect", data_structure.stellarator_variables, 0.70, 1.30 ), + 177: IterationVariable( + "f_plasma_particles_lcfs_recycled", data_structure.physics_variables, 0.01, 1.0 + ), + 178: IterationVariable( + "eta_plasma_fuelling", data_structure.physics_variables, 0.01, 1.0 + ), + 179: IterationVariable( + "molflow_plasma_fuelling_vv_injected", + data_structure.physics_variables, + 1e18, + 1e22, + ), + 180: IterationVariable( + "f_molflow_plasma_fuelling_deuterium", + data_structure.physics_variables, + 0.0, + 1.0, + ), + 181: IterationVariable( + "f_molflow_plasma_fuelling_tritium", + data_structure.physics_variables, + 0.0, + 1.0, + ), + 182: IterationVariable( + "f_molflow_plasma_fuelling_helium3", + data_structure.physics_variables, + 0.0, + 1.0, + ), } diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 6eafc3aa9..2a388e366 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -1,9 +1,9 @@ import numpy as np -ipnvars: int = 177 +ipnvars: int = 183 """total number of variables available for iteration""" -ipeqns: int = 92 +ipeqns: int = 96 """number of constraint equations available""" ipnfoms: int = 19 @@ -188,6 +188,7 @@ * (90) Lower Limit on number of stress load cycles for CS * (91) Checking if the design point is ECRH ignitable * (92) D/T/He3 ratio in fuel sums to 1 +* (93) Particle balance consistency constraint """ ixc: list[int] = None @@ -314,6 +315,7 @@ * (115) NOT USED * (116) NOT USED * (117) NOT USED +* (118) NOT USED * (119) temp_plasma_separatrix_kev: separatrix temperature calculated by the Kallenbach divertor model * (120) ttarget: Plasma temperature adjacent to divertor sheath [eV] * (121) neratio: ratio of mean SOL density at OMP to separatrix density at OMP @@ -621,6 +623,10 @@ def init_numerics(): "CS stress load cycles ", "ECRH ignitability ", "Fuel composition consistency ", + "Tritium particle balance consistency ", + "Deuterium particle balance consistency ", + "Alpha particle balance consistency ", + "Fuelling composition consistency ", ] ixc = np.array([0] * ipnvars) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index fccd7a829..5a92fe901 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -253,8 +253,15 @@ """Plasma stored magnetic energy (J)""" -burnup: float = None -"""fractional plasma burnup""" +f_plasma_fuel_burnup: float = None +"""Total fuel burnup fraction in plasma""" + + +f_plasma_tritium_burnup: float = None +"""Tritium burnup fraction in plasma""" + +f_plasma_deuterium_burnup: float = None +"""Deuterium burnup fraction in plasma""" burnup_in: float = None @@ -458,6 +465,27 @@ fusrat_total: float = None """fusion reaction rate, from beams and plasma (reactions/sec)""" +fusrat_plasma_dt: float = None +""" D-T fusion reaction rate in plasma, (reactions/sec)""" + +fusrat_dt_total: float = None +""" Total D-T fusion reaction rate from beams and plasma, (reactions/sec)""" + +fusrat_plasma_dd_helion: float = None +"""D-D fusion reaction rate (helium branch) in plasma, (reactions/sec)""" + +fusrat_plasma_dd_triton: float = None +"""D-D fusion reaction rate (tritium branch) in plasma, (reactions/sec)""" + +fusrat_plasma_dd_total: float = None +"""Total D-D fusion reaction rate in plasma, (reactions/sec)""" + +fusrat_plasma_dhe3: float = None +"""D-3He fusion reaction rate in plasma, (reactions/sec)""" + +fusrat_neutron_production_total: float = None +"""Total neutron production rate from plasma and beams (neutrons/sec)""" + fusrat_plasma_dt_profile: list[float] = None """Profile of D-T fusion reaction rate in plasma, (reactions/sec)""" @@ -1062,12 +1090,32 @@ """ -molflow_plasma_fuelling_required: float = None -"""plasma fuelling rate (nucleus-pairs/s)""" +f_plasma_particles_lcfs_recycled: float = None +"""fraction of plasma particles that are recycled at the LCFS. Recycling coefficent (R)""" + +eta_plasma_fuelling: float = None +"""fuelling efficiency (fraction of fuel particles injected that become confined in the plasma)""" + +molflow_plasma_fuelling_vv_injected: float = None +"""plasma fuelling rate into the vacuum vessel (particles/s)""" + +molflow_plasma_fuelling_vv_injected_moles: float = None +"""plasma fuelling rate into the vacuum vessel (moles/s)""" + +molflow_plasma_fuelling_loss: float = None +"""plasma fuelling rate that dosent make it to plasma (particles/s)""" +molflow_plasma_fuelling_loss_moles: float = None +"""plasma fuelling rate that dosent make it to plasma (moles/s)""" -tauratio: float = None -"""tauratio /1.0/ : ratio of He and pellet particle confinement times""" +f_molflow_plasma_fuelling_deuterium: float = None +"""fraction of plasma fuelling that is deuterium""" + +f_molflow_plasma_fuelling_tritium: float = None +"""fraction of plasma fuelling that is tritium""" + +f_molflow_plasma_fuelling_helium3: float = None +"""fraction of plasma fuelling that is helium-3""" q95_min: float = None @@ -1129,10 +1177,6 @@ """n_carbon / n_e""" -rndfuel: float = None -"""fuel burnup rate (reactions/second)""" - - f_nd_plasma_iron_argon_electron: float = None """n_highZ / n_e""" @@ -1448,7 +1492,9 @@ def init_physics_variables(): b_plasma_toroidal_profile, \ b_plasma_total, \ e_plasma_magnetic_stored, \ - burnup, \ + f_plasma_fuel_burnup, \ + f_plasma_tritium_burnup, \ + f_plasma_deuterium_burnup, \ burnup_in, \ b_plasma_vertical_required, \ c_beta, \ @@ -1494,6 +1540,13 @@ def init_physics_variables(): f_plasma_fuel_tritium, \ fusden_total, \ fusrat_total, \ + fusrat_plasma_dt, \ + fusrat_dt_total, \ + fusrat_plasma_dd_helion, \ + fusrat_plasma_dd_triton, \ + fusrat_plasma_dd_total, \ + fusrat_plasma_dhe3, \ + fusrat_neutron_production_total, \ fusrat_plasma_dt_profile, \ fusrat_plasma_dd_triton_profile, \ fusrat_plasma_dd_helion_profile, \ @@ -1616,8 +1669,15 @@ def init_physics_variables(): pden_ion_transport_loss_mw, \ q0, \ q95, \ - molflow_plasma_fuelling_required, \ - tauratio, \ + f_plasma_particles_lcfs_recycled, \ + eta_plasma_fuelling, \ + molflow_plasma_fuelling_vv_injected, \ + molflow_plasma_fuelling_vv_injected_moles, \ + molflow_plasma_fuelling_loss, \ + molflow_plasma_fuelling_loss_moles, \ + f_molflow_plasma_fuelling_deuterium, \ + f_molflow_plasma_fuelling_tritium, \ + f_molflow_plasma_fuelling_helium3, \ q95_min, \ qstar, \ rad_fraction_sol, \ @@ -1634,7 +1694,6 @@ def init_physics_variables(): rminor, \ f_nd_beam_electron, \ f_nd_plasma_carbon_electron, \ - rndfuel, \ f_nd_plasma_iron_argon_electron, \ f_nd_plasma_oxygen_electron, \ f_res_plasma_neo, \ @@ -1737,7 +1796,9 @@ def init_physics_variables(): b_plasma_toroidal_profile = [] b_plasma_total = 0.0 e_plasma_magnetic_stored = 0.0 - burnup = 0.0 + f_plasma_fuel_burnup = 0.0 + f_plasma_tritium_burnup = 0.0 + f_plasma_deuterium_burnup = 0.0 burnup_in = 0.0 b_plasma_vertical_required = 0.0 c_beta = 0.5 @@ -1783,6 +1844,13 @@ def init_physics_variables(): f_plasma_fuel_tritium = 0.5 fusden_total = 0.0 fusrat_total = 0.0 + fusrat_plasma_dt = 0.0 + fusrat_plasma_dd_helion = 0.0 + fusrat_plasma_dd_triton = 0.0 + fusrat_plasma_dd_total = 0.0 + fusrat_plasma_dhe3 = 0.0 + fusrat_neutron_production_total = 0.0 + fusrat_dt_total = 0.0 fusrat_plasma_dt_profile = [] fusrat_plasma_dd_triton_profile = [] fusrat_plasma_dd_helion_profile = [] @@ -1907,8 +1975,15 @@ def init_physics_variables(): pden_ion_transport_loss_mw = 0.0 q0 = 1.0 q95 = 0.0 - molflow_plasma_fuelling_required = 0.0 - tauratio = 1.0 + f_plasma_particles_lcfs_recycled = 0.0 + eta_plasma_fuelling = 0.0 + molflow_plasma_fuelling_vv_injected = 1e20 + molflow_plasma_fuelling_vv_injected_moles = 0.0 + molflow_plasma_fuelling_loss = 0.0 + molflow_plasma_fuelling_loss_moles = 0.0 + f_molflow_plasma_fuelling_deuterium = 0.5 + f_molflow_plasma_fuelling_tritium = 0.5 + f_molflow_plasma_fuelling_helium3 = 0.0 q95_min = 0.0 qstar = 0.0 rad_fraction_sol = 0.8 @@ -1925,7 +2000,6 @@ def init_physics_variables(): rminor = 0.0 f_nd_beam_electron = 0.005 f_nd_plasma_carbon_electron = 0.0 - rndfuel = 0.0 f_nd_plasma_iron_argon_electron = 0.0 f_nd_plasma_oxygen_electron = 0.0 f_res_plasma_neo = 0.0 diff --git a/process/main.py b/process/main.py index c37e24dad..32b1fa6f0 100644 --- a/process/main.py +++ b/process/main.py @@ -97,6 +97,7 @@ ) from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.exhaust import PlasmaExhaust +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.l_h_transition import PlasmaConfinementTransition from process.models.physics.physics import ( @@ -709,6 +710,7 @@ def __init__(self): self.plasma_inductance = PlasmaInductance() self.plasma_density_limit = PlasmaDensityLimit() self.plasma_exhaust = PlasmaExhaust() + self.plasma_fuelling = PlasmaFuelling() self.plasma_bootstrap_current = PlasmaBootstrapCurrent( plasma_profile=self.plasma_profile ) @@ -726,6 +728,7 @@ def __init__(self): plasma_confinement=self.plasma_confinement, plasma_transition=self.plasma_transition, plasma_current=self.plasma_current, + plasma_fuelling=self.plasma_fuelling, ) self.physics_detailed = DetailedPhysics( plasma_profile=self.plasma_profile, diff --git a/process/models/costs/costs.py b/process/models/costs/costs.py index be7f2198e..05e520411 100644 --- a/process/models/costs/costs.py +++ b/process/models/costs/costs.py @@ -2362,7 +2362,7 @@ def acc2272(self): # New calculation: 2 nuclei * reactions/sec * kg/nucleus * g/kg * sec/day physics_variables.wtgpd = ( 2.0e0 - * physics_variables.rndfuel + * physics_variables.fusrat_total * physics_variables.m_fuel_amu * constants.UMASS * 1000.0e0 diff --git a/process/models/physics/fuelling.py b/process/models/physics/fuelling.py new file mode 100644 index 000000000..35e8d4c6e --- /dev/null +++ b/process/models/physics/fuelling.py @@ -0,0 +1,795 @@ +import logging + +import matplotlib.pyplot as plt +import numpy as np + +import process.core.io.mfile as mf +from process.core import constants +from process.core import process_output as po +from process.data_structure import ( + numerics, + physics_variables, + reinke_variables, +) + +logger = logging.getLogger(__name__) + + +class PlasmaFuelling: + """Class to hold plasma fuelling calculations for plasma processing.""" + + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE + + def run(self): + physics_variables.molflow_plasma_fuelling_vv_injected_moles = ( + physics_variables.molflow_plasma_fuelling_vv_injected + / constants.AVOGADRO_NUMBER + ) + + physics_variables.molflow_plasma_fuelling_loss = ( + physics_variables.molflow_plasma_fuelling_vv_injected + * (1 - physics_variables.eta_plasma_fuelling) + ) + physics_variables.molflow_plasma_fuelling_loss_moles = ( + physics_variables.molflow_plasma_fuelling_loss / constants.AVOGADRO_NUMBER + ) + + physics_variables.f_plasma_fuel_burnup = self.calculate_fuel_burnup_fraction( + fusrat_total=physics_variables.fusrat_total, + molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, + ) + physics_variables.f_plasma_tritium_burnup = self.calculate_tritium_burnup_fraction( + fusrat_dt_total=physics_variables.fusrat_dt_total, + molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, + f_molflow_plasma_fuelling_tritium=physics_variables.f_molflow_plasma_fuelling_tritium, + ) + + physics_variables.f_plasma_deuterium_burnup = self.calculate_deuterium_burnup_fraction( + fusrat_plasma_dd_total=physics_variables.fusrat_plasma_dd_total, + molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, + f_molflow_plasma_fuelling_deuterium=physics_variables.f_molflow_plasma_fuelling_deuterium, + fusrat_dt_total=physics_variables.fusrat_dt_total, + fusrat_plasma_dhe3=physics_variables.fusrat_plasma_dhe3, + ) + + @staticmethod + def calculate_fuel_burnup_fraction( + fusrat_total: float, molflow_plasma_fuelling_vv_injected: float + ) -> float: + """Calculate the fuel burnup fraction + + Parameters + ---------- + fusrat_total : float + Total fusion rate (particles/s). + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate into vacuum vessel (particles/s). + + Returns + ------- + float Fuel burnup fraction (dimensionless). + + + Notes + ----- + The fusion rate is multiplied by two to convert from nucleus pairs to particles, as the fuelling rate is in particles/s. + + """ + + return 2 * fusrat_total / molflow_plasma_fuelling_vv_injected + + @staticmethod + def calculate_tritium_burnup_fraction( + fusrat_dt_total: float, + molflow_plasma_fuelling_vv_injected: float, + f_molflow_plasma_fuelling_tritium: float, + ) -> float: + """Calculate the tritium burnup fraction + + Parameters + ---------- + fusrat_dt_total : float + Total DT fusion rate (particles/s). + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate into vacuum vessel (particles/s). + f_molflow_plasma_fuelling_tritium : float + Fraction of tritium in the plasma fuelling. + + Returns + ------- + float Tritium burnup fraction (dimensionless). + + Notes + ----- + The fusion rate is multiplied by two to convert from nucleus pairs to particles, as the fuelling rate is in particles/s. + + """ + + return fusrat_dt_total / ( + molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_tritium + ) + + @staticmethod + def calculate_deuterium_burnup_fraction( + fusrat_dt_total: float, + molflow_plasma_fuelling_vv_injected: float, + f_molflow_plasma_fuelling_deuterium: float, + fusrat_plasma_dd_total: float, + fusrat_plasma_dhe3: float, + ) -> float: + """Calculate the deuterium burnup fraction + + Parameters + ---------- + fusrat_dt_total : float + Total DT fusion rate (particles/s). + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate into vacuum vessel (particles/s). + f_molflow_plasma_fuelling_deuterium : float + Fraction of deuterium in the plasma fuelling. + fusrat_plasma_dd_total : float + Total deuterium consumption rate from DD fusion (particles/s). + fusrat_plasma_dhe3 : float + Deuterium consumption rate from D-He3 fusion (particles/s). + + Returns + ------- + float Deuterium burnup fraction (dimensionless). + + Notes + ----- + The fusion rate is multiplied by two to convert from nucleus pairs to particles, as the fuelling rate is in particles/s. + + """ + + return (fusrat_dt_total + 2 * fusrat_plasma_dd_total + fusrat_plasma_dhe3) / ( + molflow_plasma_fuelling_vv_injected * f_molflow_plasma_fuelling_deuterium + ) + + @staticmethod + def calculate_plasma_tritium_flow_rate( + f_molflow_plasma_fuelling_tritium: float, + eta_plasma_fuelling: float, + molflow_plasma_fuelling_vv_injected: float, + fusrat_dt_total: float, + fusrat_plasma_dd_triton: float, + t_energy_confinement: float, + f_plasma_particles_lcfs_recycled: float, + nd_plasma_fuel_ions_vol_avg: float, + vol_plasma: float, + f_plasma_fuel_tritium: float, + ) -> float: + """Calculate the tritium flow rate in the plasma exhaust. + + Parameters + ---------- + f_molflow_plasma_fuelling_tritium : float + Fraction of tritium in the plasma fuelling. + eta_plasma_fuelling : float + Fuelling rate efficiency. + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate (particles/s). + fusrat_dt_total : float + Total DT fusion rate (particles/s). + fusrat_plasma_dd_triton : float + Tritium production rate from DD fusion (particles/s). + t_energy_confinement : float + Energy confinement time (s). + f_plasma_particles_lcfs_recycled : float + Fraction of plasma particles recycled at the LCFS. + nd_plasma_fuel_ions_vol_avg : float + Volume-averaged density of fuel ions in the plasma (particles/m^3). + vol_plasma : float + Plasma volume (m^3). + f_plasma_fuel_tritium : float + Fraction of tritium in the plasma fuel. + + Returns + ------- + float + Tritium flow rate in the plasma exhaust (particles/s). + + """ + + return ( + ( + f_molflow_plasma_fuelling_tritium + * eta_plasma_fuelling + * molflow_plasma_fuelling_vv_injected + ) + - fusrat_dt_total + + fusrat_plasma_dd_triton + - ( + (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_tritium) + / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) + ) + ) + + @staticmethod + def calculate_plasma_deuterium_flow_rate( + f_molflow_plasma_fuelling_deuterium: float, + eta_plasma_fuelling: float, + molflow_plasma_fuelling_vv_injected: float, + fusrat_dt_total: float, + fusrat_plasma_dhe3: float, + fusrat_plasma_dd_total: float, + t_energy_confinement: float, + f_plasma_particles_lcfs_recycled: float, + nd_plasma_fuel_ions_vol_avg: float, + vol_plasma: float, + f_plasma_fuel_deuterium: float, + ) -> float: + """Calculate the deuterium flow rate in the plasma exhaust. + + Parameters + ---------- + f_molflow_plasma_fuelling_deuterium : float + Fraction of deuterium in the plasma fuelling. + eta_plasma_fuelling : float + Fuelling rate efficiency. + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate (particles/s). + fusrat_dt_total : float + Total DT fusion rate (particles/s). + fusrat_plasma_dhe3 : float + Deuterium consumption rate from D-He3 fusion (particles/s). + fusrat_plasma_dd_total : float + Total deuterium consumption rate from DD fusion (particles/s). + t_energy_confinement : float + Energy confinement time (s). + f_plasma_particles_lcfs_recycled : float + Fraction of plasma particles recycled at the LCFS. + nd_plasma_fuel_ions_vol_avg : float + Volume-averaged density of fuel ions in the plasma (particles/m^3). + vol_plasma : float + Plasma volume (m^3). + f_plasma_fuel_deuterium : float + Fraction of deuterium in the plasma fuel. + + Returns + ------- + float + Deuterium flow rate in the plasma exhaust (particles/s). + + + """ + return ( + ( + f_molflow_plasma_fuelling_deuterium + * eta_plasma_fuelling + * molflow_plasma_fuelling_vv_injected + ) + - fusrat_dt_total + - 2 * fusrat_plasma_dd_total + - fusrat_plasma_dhe3 + - ( + (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_deuterium) + / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) + ) + ) + + @staticmethod + def calculate_plasma_helium3_flow_rate( + f_molflow_plasma_fuelling_helium3: float, + eta_plasma_fuelling: float, + molflow_plasma_fuelling_vv_injected: float, + fusrat_plasma_dhe3: float, + t_energy_confinement: float, + f_plasma_particles_lcfs_recycled: float, + nd_plasma_fuel_ions_vol_avg: float, + vol_plasma: float, + f_plasma_fuel_helium3: float, + ) -> float: + """Calculate the helium-3 flow rate in the plasma exhaust. + + Parameters + ---------- + f_molflow_plasma_fuelling_helium3 : float + Fraction of helium-3 in the plasma fuelling. + eta_plasma_fuelling : float + Fuelling rate efficiency. + molflow_plasma_fuelling_vv_injected : float + Total fuelling rate (particles/s). + fusrat_plasma_dhe3 : float + Deuterium consumption rate from D-He3 fusion (particles/s). + t_energy_confinement : float + Energy confinement time (s). + f_plasma_particles_lcfs_recycled : float + Fraction of plasma particles recycled at the LCFS. + nd_plasma_fuel_ions_vol_avg : float + Volume-averaged density of fuel ions in the plasma (particles/m^3). + vol_plasma : float + Plasma volume (m^3). + f_plasma_fuel_helium3 : float + Fraction of helium-3 in the plasma fuel. + + Returns + ------- + float + Helium-3 flow rate in the plasma exhaust (particles/s). + + """ + + return ( + ( + f_molflow_plasma_fuelling_helium3 + * eta_plasma_fuelling + * molflow_plasma_fuelling_vv_injected + ) + + fusrat_plasma_dhe3 + - ( + (nd_plasma_fuel_ions_vol_avg * vol_plasma * f_plasma_fuel_helium3) + / (t_energy_confinement / (1 - f_plasma_particles_lcfs_recycled)) + ) + ) + + @staticmethod + def calculate_plasma_alphas_flow_rate( + fusrat_dt_total: float, + fusrat_plasma_dhe3: float, + t_energy_confinement: float, + f_t_alpha_energy_confinement: float, + nd_plasma_alphas_vol_avg: float, + vol_plasma: float, + ) -> float: + """Calculate the alpha particle flow rate in the plasma exhaust. + + Parameters + ---------- + fusrat_dt_total : float + Total DT fusion rate (particles/s). + fusrat_plasma_dhe3 : float + Deuterium consumption rate from D-He3 fusion (particles/s). + t_energy_confinement : float + Energy confinement time (s). + f_t_alpha_energy_confinement : float + Ratio of alpha particle confinement time to energy confinement time (dimensionless). + nd_plasma_alphas_vol_avg : float + Volume-averaged density of alpha particles in the plasma (particles/m^3). + vol_plasma : float + Plasma volume (m^3). + + Returns + ------- + float + Alpha particle flow rate in the plasma exhaust (particles/s). + + """ + + # Alpha particle balance + + return ( + fusrat_dt_total + + fusrat_plasma_dhe3 + - (nd_plasma_alphas_vol_avg * vol_plasma) + / (t_energy_confinement * f_t_alpha_energy_confinement) + ) + + def plot_tritium_flow_contour(self, axis: plt.Axes, mfile: mf.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): + tritium_flow[i, j] = self.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 + ), + 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(self, axis: plt.Axes, mfile: mf.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] = self.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 + ), + 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(self, axis: plt.Axes, mfile: mf.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] = self.calculate_plasma_alphas_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_vol_avg=mfile.get( + "nd_plasma_alphas_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_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(self, axis: plt.Axes, mfile: mf.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] = self.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), + 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(self, fig: plt.Figure, mfile: mf.MFile, scan: int): + """Plot fuelling information.""" + msg = ( + f"$\\mathbf{{Plasma \\ Fuelling \\ Information:}}$\n\n" + f"Total fuelling rate:" + f"{mfile.get('molflow_plasma_fuelling_vv_injected', scan=scan):.4e} particles/s\n" + f"Total fuelling rate: " + f"{mfile.get('molflow_plasma_fuelling_vv_injected_moles', scan=scan):.4e} moles/s\n" + f"Total fuelling loss: " + f"{mfile.get('molflow_plasma_fuelling_loss', scan=scan):.4e} particles/s\n" + f"Total fuelling loss: " + 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 output_fuelling_info(self): + """Output fuelling information to mfile.""" + po.oheadr(self.outfile, "Plasma Fuelling") + po.ovarre( + self.outfile, + "Fuelling rate (nucleus-pairs/s)", + "(molflow_plasma_fuelling_vv_injected)", + physics_variables.molflow_plasma_fuelling_vv_injected, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuelling rate (moles/s)", + "(molflow_plasma_fuelling_vv_injected_moles)", + physics_variables.molflow_plasma_fuelling_vv_injected_moles, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuelling loss (nucleus-pairs/s)", + "(molflow_plasma_fuelling_loss)", + physics_variables.molflow_plasma_fuelling_loss, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuelling loss (moles/s)", + "(molflow_plasma_fuelling_loss_moles)", + physics_variables.molflow_plasma_fuelling_loss_moles, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of plasma fuelling that is deuterium", + "(f_molflow_plasma_fuelling_deuterium)", + physics_variables.f_molflow_plasma_fuelling_deuterium, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of plasma fuelling that is tritium", + "(f_molflow_plasma_fuelling_tritium)", + physics_variables.f_molflow_plasma_fuelling_tritium, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of plasma fuelling that is helium-3", + "(f_molflow_plasma_fuelling_helium3)", + physics_variables.f_molflow_plasma_fuelling_helium3, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuelling efficiency", + "(eta_plasma_fuelling)", + physics_variables.eta_plasma_fuelling, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of plasma particles recycled at the LCFS", + "(f_plasma_particles_lcfs_recycled)", + physics_variables.f_plasma_particles_lcfs_recycled, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuel burn-up rate (reactions/s)", + "(fusrat_total)", + physics_variables.fusrat_total, + "OP ", + ) + po.ovarrf( + self.outfile, + "Total fuel burn-up fraction", + "(f_plasma_fuel_burnup)", + physics_variables.f_plasma_fuel_burnup, + "OP ", + ) + po.ovarrf( + self.outfile, + "Tritium burn-up fraction", + "(f_plasma_tritium_burnup)", + physics_variables.f_plasma_tritium_burnup, + "OP ", + ) + po.ovarrf( + self.outfile, + "Deuterium burn-up fraction", + "(f_plasma_deuterium_burnup)", + physics_variables.f_plasma_deuterium_burnup, + "OP ", + ) + + if 78 in numerics.icc: + po.osubhd(self.outfile, "Reinke Criterion :") + po.ovarin( + self.outfile, + "index of impurity to be iterated for divertor detachment", + "(impvardiv)", + reinke_variables.impvardiv, + ) + po.ovarre( + self.outfile, + "Minimum Impurity fraction from Reinke", + "(fzmin)", + reinke_variables.fzmin, + "OP ", + ) + po.ovarre( + self.outfile, + "Actual Impurity fraction", + "(fzactual)", + reinke_variables.fzactual, + ) diff --git a/process/models/physics/fusion_reactions.py b/process/models/physics/fusion_reactions.py index 21c111dc2..6c9be402a 100644 --- a/process/models/physics/fusion_reactions.py +++ b/process/models/physics/fusion_reactions.py @@ -328,6 +328,10 @@ def dhe3_reaction(self): alpha_rate_density = fusion_rate_density proton_rate_density = fusion_rate_density # Proton production rate [m^3/second] + physics_variables.fusrat_plasma_dhe3 = ( + fusion_rate_density * physics_variables.vol_plasma + ) + # Update the cumulative D-3He power density self.dhe3_power_density = fusion_power_density @@ -427,6 +431,10 @@ def dd_helion_reaction(self): alpha_rate_density = 0.0 proton_rate_density = 0.0 + physics_variables.fusrat_plasma_dd_helion = ( + fusion_rate_density * physics_variables.vol_plasma + ) + # Update the cumulative D-D power density self.dd_power_density += fusion_power_density @@ -521,6 +529,10 @@ def dd_triton_reaction(self): alpha_rate_density = 0.0 proton_rate_density = fusion_rate_density # Proton production rate [m^3/second] + physics_variables.fusrat_plasma_dd_triton = ( + fusion_rate_density * physics_variables.vol_plasma + ) + # Update the cumulative D-D power density self.dd_power_density += fusion_power_density diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 34f987e42..3d4f6f044 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -30,6 +30,7 @@ ) from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.exhaust import PlasmaExhaust +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.l_h_transition import PlasmaConfinementTransition from process.models.physics.plasma_current import PlasmaCurrent @@ -214,6 +215,7 @@ def __init__( plasma_confinement: PlasmaConfinementTime, plasma_transition: PlasmaConfinementTransition, plasma_current: PlasmaCurrent, + plasma_fuelling: PlasmaFuelling, ): self.outfile = constants.NOUT self.mfile = constants.MFILE @@ -227,6 +229,7 @@ def __init__( self.confinement = plasma_confinement self.plasma_transition = plasma_transition self.current = plasma_current + self.fuelling = plasma_fuelling def output(self): self.calculate_effective_charge_ionisation_profiles() @@ -655,6 +658,7 @@ def run(self): + (1.0 / (1.0 - constants.DT_NEUTRON_ENERGY_FRACTION)) * physics_variables.p_beam_alpha_mw ) + physics_variables.p_beam_neutron_mw = physics_variables.p_beam_alpha_mw * ( constants.DT_NEUTRON_ENERGY_FRACTION / (1 - constants.DT_NEUTRON_ENERGY_FRACTION) @@ -671,6 +675,21 @@ def run(self): physics_variables.fusrat_total = ( physics_variables.fusden_total * physics_variables.vol_plasma ) + physics_variables.fusrat_plasma_dt = (physics_variables.p_plasma_dt_mw * 1e6) / ( + constants.D_T_ENERGY + ) + physics_variables.fusrat_plasma_dd_total = ( + physics_variables.fusrat_plasma_dd_helion + + physics_variables.fusrat_plasma_dd_triton + ) + + physics_variables.fusrat_neutron_production_total = ( + physics_variables.fusrat_plasma_dd_helion + physics_variables.fusrat_dt_total + ) + + physics_variables.fusrat_dt_total = ( + physics_variables.p_dt_total_mw * 1e6 / (constants.D_T_ENERGY) + ) # Create some derived values and add beam contribution to fusion power ( @@ -923,28 +942,17 @@ def run(self): 0.5e0 * physics_variables.ind_plasma * physics_variables.plasma_current**2 ) - # Calculate auxiliary physics related information - sbar = 1.0e0 ( - physics_variables.burnup, - physics_variables.figmer, - physics_variables.fusrat, - physics_variables.molflow_plasma_fuelling_required, - physics_variables.rndfuel, physics_variables.t_alpha_confinement, physics_variables.f_alpha_energy_confinement, ) = self.phyaux( - physics_variables.aspect, - physics_variables.nd_plasma_fuel_ions_vol_avg, - physics_variables.fusden_total, physics_variables.fusden_alpha_total, - physics_variables.plasma_current, - sbar, physics_variables.nd_plasma_alphas_vol_avg, physics_variables.t_energy_confinement, - physics_variables.vol_plasma, ) + self.fuelling.run() + physics_variables.ntau, physics_variables.nTtau = ( self.confinement.calculate_double_and_triple_product( nd_plasma_electrons_vol_avg=physics_variables.nd_plasma_electrons_vol_avg, @@ -1420,57 +1428,30 @@ def plasma_composition(): @staticmethod @nb.njit(cache=True) def phyaux( - aspect: float, - nd_plasma_fuel_ions_vol_avg: float, - fusden_total: float, fusden_alpha_total: float, - plasma_current: float, - sbar: float, nd_plasma_alphas_vol_avg: float, t_energy_confinement: float, - vol_plasma: float, - ) -> tuple[float, float, float, float, float, float, float, float]: + ) -> tuple[float, float]: """Auxiliary physics quantities Parameters ---------- - aspect : float - Plasma aspect ratio. - nd_plasma_fuel_ions_vol_avg : float - Fuel ion density (/m3). - fusden_total : float - Fusion reaction rate from plasma and beams (/m3/s). fusden_alpha_total : float - Alpha particle production rate (/m3/s). - plasma_current : float - Plasma current (A). - sbar : float - Exponent for aspect ratio (normally 1). + Total alpha particle production rate density (m^-3 s^-1). nd_plasma_alphas_vol_avg : float - Alpha ash density (/m3). + Volume averaged alpha particle density (m^-3). t_energy_confinement : float - Global energy confinement time (s). - vol_plasma : float - Plasma volume (m3). + Energy confinement time (s). Returns ------- tuple A tuple containing: - - burnup (float): Fractional plasma burnup. - - figmer (float): Physics figure of merit. - - fusrat (float): Number of fusion reactions per second. - - molflow_plasma_fuelling_required (float): Fuelling rate for D-T (nucleus-pairs/sec). - - rndfuel (float): Fuel burnup rate (reactions/s). - t_alpha_confinement (float): Alpha particle confinement time (s). - f_alpha_energy_confinement (float): Fraction of alpha energy confinement. This subroutine calculates extra physics related items needed by other parts of the code. """ - figmer = 1e-6 * plasma_current * aspect**sbar - - # Fusion reactions per second - fusrat = fusden_total * vol_plasma # Alpha particle confinement time (s) # Number of alphas / alpha production rate @@ -1479,39 +1460,9 @@ def phyaux( else: # only likely if DD is only active fusion reaction t_alpha_confinement = 0.0 - # Fractional burnup - # (Consider detailed model in: G. L. Jackson, V. S. Chan, R. D. Stambaugh, - # Fusion Science and Technology, vol.64, no.1, July 2013, pp.8-12) - # The ratio of ash to fuel particle confinement times is given by - # tauratio - # Possible logic... - # burnup = fuel ion-pairs burned/m3 / initial fuel ion-pairs/m3; - # fuel ion-pairs burned/m3 = alpha particles/m3 (for both D-T and D-He3 reactions) - # initial fuel ion-pairs/m3 = burnt fuel ion-pairs/m3 + unburnt fuel-ion pairs/m3 - # Remember that unburnt fuel-ion pairs/m3 = 0.5 * unburnt fuel-ions/m3 - if physics_variables.burnup_in <= 1.0e-9: - burnup = ( - nd_plasma_alphas_vol_avg - / (nd_plasma_alphas_vol_avg + 0.5 * nd_plasma_fuel_ions_vol_avg) - / physics_variables.tauratio - ) - else: - burnup = physics_variables.burnup_in - - # Fuel burnup rate (reactions/second) (previously Amps) - rndfuel = fusrat - - # Required fuelling rate (fuel ion pairs/second) (previously Amps) - molflow_plasma_fuelling_required = rndfuel / burnup - f_alpha_energy_confinement = t_alpha_confinement / t_energy_confinement return ( - burnup, - figmer, - fusrat, - molflow_plasma_fuelling_required, - rndfuel, t_alpha_confinement, f_alpha_energy_confinement, ) @@ -2641,6 +2592,56 @@ def outplas(self): physics_variables.fusrat_total, "OP ", ) + po.ovarre( + self.outfile, + "D-T Fusion rate: total (reactions/sec)", + "(fusrat_dt_total)", + physics_variables.fusrat_dt_total, + "OP ", + ) + po.ovarre( + self.outfile, + "D-T Fusion rate: plasma (reactions/sec)", + "(fusrat_plasma_dt)", + physics_variables.fusrat_plasma_dt, + "OP ", + ) + + po.ovarre( + self.outfile, + "D-D -> 3He Fusion rate: plasma (reactions/sec)", + "(fusrat_plasma_dd_helion)", + physics_variables.fusrat_plasma_dd_helion, + "OP ", + ) + po.ovarre( + self.outfile, + "D-D -> T Fusion rate: plasma (reactions/sec)", + "(fusrat_plasma_dd_triton)", + physics_variables.fusrat_plasma_dd_triton, + "OP ", + ) + po.ovarre( + self.outfile, + "D-D Fusion rate: total (reactions/sec)", + "(fusrat_plasma_dd_total)", + physics_variables.fusrat_plasma_dd_total, + "OP ", + ) + po.ovarre( + self.outfile, + "D-3He Fusion rate: total (reactions/sec)", + "(fusrat_plasma_dhe3)", + physics_variables.fusrat_plasma_dhe3, + "OP ", + ) + po.ovarre( + self.outfile, + "Neutron production rate: total (particles/sec)", + "(fusrat_neutron_production_total)", + physics_variables.fusrat_neutron_production_total, + "OP ", + ) po.ovarre( self.outfile, "Fusion rate density: total (reactions/m3/sec)", @@ -3242,56 +3243,7 @@ def outplas(self): if stellarator_variables.istell == 0: self.plasma_bootstrap_current.output() - po.osubhd(self.outfile, "Fuelling :") - po.ovarre( - self.outfile, - "Ratio of He and pellet particle confinement times", - "(tauratio)", - physics_variables.tauratio, - ) - po.ovarre( - self.outfile, - "Fuelling rate (nucleus-pairs/s)", - "(molflow_plasma_fuelling_required)", - physics_variables.molflow_plasma_fuelling_required, - "OP ", - ) - po.ovarre( - self.outfile, - "Fuel burn-up rate (reactions/s)", - "(rndfuel)", - physics_variables.rndfuel, - "OP ", - ) - po.ovarrf( - self.outfile, - "Burn-up fraction", - "(burnup)", - physics_variables.burnup, - "OP ", - ) - - if 78 in numerics.icc: - po.osubhd(self.outfile, "Reinke Criterion :") - po.ovarin( - self.outfile, - "index of impurity to be iterated for divertor detachment", - "(impvardiv)", - reinke_variables.impvardiv, - ) - po.ovarre( - self.outfile, - "Minimum Impurity fraction from Reinke", - "(fzmin)", - reinke_variables.fzmin, - "OP ", - ) - po.ovarre( - self.outfile, - "Actual Impurity fraction", - "(fzactual)", - reinke_variables.fzactual, - ) + self.fuelling.output_fuelling_info() @staticmethod @nb.njit(cache=True) diff --git a/process/models/stellarator/stellarator.py b/process/models/stellarator/stellarator.py index 1df7adc0b..ff59e1d1c 100644 --- a/process/models/stellarator/stellarator.py +++ b/process/models/stellarator/stellarator.py @@ -2355,25 +2355,18 @@ def st_phys(self, output): # Calculate auxiliary physics related information # for the rest of the code - sbar = 1.0e0 ( - physics_variables.burnup, - physics_variables.figmer, - _fusrat, - physics_variables.molflow_plasma_fuelling_required, - physics_variables.rndfuel, physics_variables.t_alpha_confinement, physics_variables.f_alpha_energy_confinement, ) = self.physics.phyaux( - physics_variables.aspect, - physics_variables.nd_plasma_fuel_ions_vol_avg, - physics_variables.fusden_total, physics_variables.fusden_alpha_total, - physics_variables.plasma_current, - sbar, physics_variables.nd_plasma_alphas_vol_avg, physics_variables.t_energy_confinement, - physics_variables.vol_plasma, + ) + + physics_variables.f_plasma_fuel_burnup = self.physics.fuelling.calculate_fuel_burnup_fraction( + fusrat_total=physics_variables.fusrat_total, + molflow_plasma_fuelling_vv_injected=physics_variables.molflow_plasma_fuelling_vv_injected, ) # Calculate the neoclassical sanity check with PROCESS parameters diff --git a/process/models/vacuum.py b/process/models/vacuum.py index 345153a9a..71497d71d 100644 --- a/process/models/vacuum.py +++ b/process/models/vacuum.py @@ -54,7 +54,7 @@ def run(self, output: bool = False): # MDK Check this!! gasld = ( 2.0e0 - * physics_variables.molflow_plasma_fuelling_required + * physics_variables.molflow_plasma_fuelling_vv_injected * physics_variables.m_fuel_amu * constants.UMASS ) @@ -123,7 +123,7 @@ def vacuum_simple(self, output) -> float: # One ITER torus cryopump has a throughput of 50 Pa m3/s = 1.2155e+22 molecules/s # Issue #304 n_iter_vacuum_pumps = ( - physics_variables.molflow_plasma_fuelling_required + physics_variables.molflow_plasma_fuelling_vv_injected / vacuum_variables.molflow_vac_pumps ) @@ -167,8 +167,8 @@ def vacuum_simple(self, output) -> float: process_output.ovarre( self.outfile, "Plasma fuelling rate (nucleus-pairs/s)", - "(molflow_plasma_fuelling_required)", - physics_variables.molflow_plasma_fuelling_required, + "(molflow_plasma_fuelling_vv_injected)", + physics_variables.molflow_plasma_fuelling_vv_injected, "OP ", ) process_output.ocmmnt( diff --git a/tests/unit/test_costs_1990.py b/tests/unit/test_costs_1990.py index 1843269b9..ebeb7c3f8 100644 --- a/tests/unit/test_costs_1990.py +++ b/tests/unit/test_costs_1990.py @@ -164,7 +164,7 @@ def test_acc2272(monkeypatch, costs): :param monkeypatch: Mock fixture :type monkeypatch: object """ - monkeypatch.setattr(physics_variables, "rndfuel", 7.158e20) + monkeypatch.setattr(physics_variables, "fusrat_total", 7.158e20) monkeypatch.setattr(physics_variables, "m_fuel_amu", 2.5) monkeypatch.setattr(cost_variables, "fkind", 1) monkeypatch.setattr(cost_variables, "c2271", 0) @@ -4380,7 +4380,7 @@ class Acc2272Param(NamedTuple): wtgpd: Any = None - rndfuel: Any = None + fusrat_total: Any = None m_fuel_amu: Any = None @@ -4406,7 +4406,7 @@ class Acc2272Param(NamedTuple): gain=0, edrive=5000000, wtgpd=0, - rndfuel=7.0799717510383796e20, + fusrat_total=7.0799717510383796e20, m_fuel_amu=2.5, c227=0, c2272=0, @@ -4422,7 +4422,7 @@ class Acc2272Param(NamedTuple): gain=0, edrive=5000000, wtgpd=507.88376577416528, - rndfuel=7.0777619721108953e20, + fusrat_total=7.0777619721108953e20, m_fuel_amu=2.5, c227=284.96904049038437, c2272=114.02873340990777, @@ -4459,7 +4459,7 @@ def test_acc2272_rut(acc2272param, monkeypatch, costs): monkeypatch.setattr(physics_variables, "wtgpd", acc2272param.wtgpd) - monkeypatch.setattr(physics_variables, "rndfuel", acc2272param.rndfuel) + monkeypatch.setattr(physics_variables, "fusrat_total", acc2272param.fusrat_total) monkeypatch.setattr(physics_variables, "m_fuel_amu", acc2272param.m_fuel_amu) diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 9e4b7157a..b48959fe5 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -23,6 +23,7 @@ ) from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.exhaust import PlasmaExhaust +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.l_h_transition import PlasmaConfinementTransition from process.models.physics.physics import ( @@ -66,6 +67,7 @@ def physics(): PlasmaConfinementTime(), PlasmaConfinementTransition(), PlasmaCurrent(), + plasma_fuelling=PlasmaFuelling(), ) @@ -2147,46 +2149,12 @@ def test_vscalc(voltsecondreqparam): class PhyauxParam(NamedTuple): - tauratio: Any = None - - burnup_in: Any = None - - aspect: Any = None - - nd_plasma_electrons_vol_avg: Any = None - - te: Any = None - - nd_plasma_fuel_ions_vol_avg: Any = None - nd_plasma_alphas_vol_avg: Any = None - fusden_total: Any = None - fusden_alpha_total: Any = None - plasma_current: Any = None - - sbar: Any = None - t_energy_confinement: Any = None - vol_plasma: Any = None - - expected_burnup: Any = None - - expected_ntau: Any = None - - expected_nTtau: Any = None - - expected_figmer: Any = None - - expected_fusrat: Any = None - - expected_molflow_plasma_fuelling_required: Any = None - - expected_rndfuel: Any = None - expected_t_alpha_confinement: Any = None @@ -2194,41 +2162,15 @@ class PhyauxParam(NamedTuple): "phyauxparam", ( PhyauxParam( - tauratio=1, - burnup_in=0, - aspect=3, - nd_plasma_fuel_ions_vol_avg=5.8589175702454272e19, - nd_plasma_alphas_vol_avg=7.5e18, - fusden_total=1.9852091609123786e17, fusden_alpha_total=1.973996644759543e17, - plasma_current=18398455.678867526, - sbar=1, + nd_plasma_alphas_vol_avg=7.5e18, t_energy_confinement=3.401323521525641, - vol_plasma=1888.1711539956691, - expected_burnup=0.20383432558954095, - expected_figmer=55.195367036602576, - expected_fusrat=3.7484146722826997e20, - expected_molflow_plasma_fuelling_required=1.8389516394951276e21, - expected_rndfuel=3.7484146722826997e20, expected_t_alpha_confinement=37.993985551650177, ), PhyauxParam( - tauratio=1, - burnup_in=0, - aspect=3, - nd_plasma_fuel_ions_vol_avg=5.8576156204039725e19, - nd_plasma_alphas_vol_avg=7.5e18, - fusden_total=1.9843269653375773e17, fusden_alpha_total=1.9731194318497056e17, - plasma_current=18398455.678867526, - sbar=1, + nd_plasma_alphas_vol_avg=7.5e18, t_energy_confinement=3.402116961408892, - vol_plasma=1888.1711539956691, - expected_burnup=0.20387039462081086, - expected_figmer=55.195367036602576, - expected_fusrat=3.7467489360461772e20, - expected_molflow_plasma_fuelling_required=1.8378092331723546e21, - expected_rndfuel=3.7467489360461772e20, expected_t_alpha_confinement=38.010876984618747, ), ), @@ -2246,42 +2188,15 @@ def test_phyaux(phyauxparam, monkeypatch, physics): :type monkeypatch: _pytest.monkeypatch.monkeypatch """ - monkeypatch.setattr(physics_variables, "tauratio", phyauxparam.tauratio) - - monkeypatch.setattr(physics_variables, "burnup_in", phyauxparam.burnup_in) - ( - burnup, - figmer, - fusrat, - molflow_plasma_fuelling_required, - rndfuel, t_alpha_confinement, _, ) = physics.phyaux( - aspect=phyauxparam.aspect, - nd_plasma_fuel_ions_vol_avg=phyauxparam.nd_plasma_fuel_ions_vol_avg, nd_plasma_alphas_vol_avg=phyauxparam.nd_plasma_alphas_vol_avg, - fusden_total=phyauxparam.fusden_total, fusden_alpha_total=phyauxparam.fusden_alpha_total, - plasma_current=phyauxparam.plasma_current, - sbar=phyauxparam.sbar, t_energy_confinement=phyauxparam.t_energy_confinement, - vol_plasma=phyauxparam.vol_plasma, - ) - - assert burnup == pytest.approx(phyauxparam.expected_burnup) - - assert figmer == pytest.approx(phyauxparam.expected_figmer) - - assert fusrat == pytest.approx(phyauxparam.expected_fusrat) - - assert molflow_plasma_fuelling_required == pytest.approx( - phyauxparam.expected_molflow_plasma_fuelling_required ) - assert rndfuel == pytest.approx(phyauxparam.expected_rndfuel) - assert t_alpha_confinement == pytest.approx(phyauxparam.expected_t_alpha_confinement) diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 2c5dfe595..11b7a49bf 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -33,6 +33,7 @@ ) from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.exhaust import PlasmaExhaust +from process.models.physics.fuelling import PlasmaFuelling from process.models.physics.l_h_transition import PlasmaConfinementTransition from process.models.physics.physics import ( Physics, @@ -98,6 +99,7 @@ def stellarator(): PlasmaConfinementTime(), PlasmaConfinementTransition(), PlasmaCurrent(), + plasma_fuelling=PlasmaFuelling(), ), Neoclassics(), plasma_beta=PlasmaBeta(), diff --git a/tests/unit/test_vacuum.py b/tests/unit/test_vacuum.py index a20ef5a1e..3af9ff46e 100644 --- a/tests/unit/test_vacuum.py +++ b/tests/unit/test_vacuum.py @@ -42,7 +42,7 @@ def test_simple_model(self, monkeypatch, vacuum): :type tfcoil: tests.unit.test_tfcoil.tfcoil (functional fixture) """ monkeypatch.setattr( - pv, "molflow_plasma_fuelling_required", 7.5745668997694112e22 + pv, "molflow_plasma_fuelling_vv_injected", 7.5745668997694112e22 ) monkeypatch.setattr(pv, "a_plasma_surface", 1500.3146527709359) monkeypatch.setattr(tfv, "n_tf_coils", 18)