diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index af8b0dd94..cc0e3db50 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -2776,8 +2776,8 @@ def plot_main_plasma_information( # Add divertor information textstr_div = ( f"\n$P_{{\\text{{sep}}}}$: {mfile.get('p_plasma_separatrix_mw', scan=scan):.2f} MW \n" - f"$\\frac{{P_{{\\text{{sep}}}}}}{{R}}$: {mfile.get('p_plasma_separatrix_mw/rmajor', scan=scan):.2f} MW/m \n" - f"$\\frac{{P_{{\\text{{sep}}}}}}{{B_T q_a R}}$: {mfile.get('pdivtbt_over_qar', scan=scan):.2f} MW T/m " + f"$\\frac{{P_{{\\text{{sep}}}}}}{{R}}$: {mfile.get('p_plasma_separatrix_rmajor_mw', scan=scan):.2f} MW/m \n" + f"$\\frac{{P_{{\\text{{sep}}}}}}{{B_T q_a R}}$: {mfile.get('p_div_bt_q_aspect_rmajor_mw', scan=scan):.2f} MW T/m " ) axis.text( diff --git a/process/core/io/variable_metadata.py b/process/core/io/variable_metadata.py index 3e9dfd5d2..2e4cdbdf4 100644 --- a/process/core/io/variable_metadata.py +++ b/process/core/io/variable_metadata.py @@ -247,7 +247,7 @@ class VariableMetadata: "p_plasma_rad_mw": VariableMetadata( latex=r"$P_{\mathrm{rad}}$ [$MW$]", description="Radiation power", units="MW" ), - "pdivtbt_over_qar": VariableMetadata( + "p_div_bt_q_aspect_rmajor_mw": VariableMetadata( latex=r"$\frac{P_{\mathrm{sep}}B_T}{q_{95}AR_{\mathrm{maj}}}$ [$MWTm^{-1}$]", description="", units="MWTm^{-1}", diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 6ff494442..92c9d9bb9 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -911,6 +911,12 @@ p_plasma_separatrix_mw: float = None """power to conducted to the divertor region (MW)""" +p_plasma_separatrix_rmajor_mw: float = None +"""Power to conducted to the divertor region per major radius (MW/m)""" + +p_div_bt_q_aspect_rmajor_mw: float = None +"""EU DEMO divertor protection parameter (MW/T/m)""" + p_div_lower_separatrix_mw: float = None """Separatrix power conducted to the lower divertor region (calculated if `i_single_null = 0`) (MW)""" @@ -1567,6 +1573,8 @@ def init_physics_variables(): p_dd_total_mw, \ p_dhe3_total_mw, \ p_plasma_separatrix_mw, \ + p_plasma_separatrix_rmajor_mw, \ + p_div_bt_q_aspect_rmajor_mw, \ p_div_lower_separatrix_mw, \ p_div_upper_separatrix_mw, \ p_div_separatrix_max_mw, \ @@ -1838,6 +1846,9 @@ def init_physics_variables(): p_dd_total_mw = 0.0 p_dhe3_total_mw = 0.0 p_plasma_separatrix_mw = 0.0 + p_plasma_separatrix_rmajor_mw = 0.0 + p_div_bt_q_aspect_rmajor_mw = 0.0 + p_div_lower_separatrix_mw = 0.0 p_div_lower_separatrix_mw = 0.0 p_div_upper_separatrix_mw = 0.0 p_div_separatrix_max_mw = 0.0 diff --git a/process/main.py b/process/main.py index 7388e750c..97e2db0eb 100644 --- a/process/main.py +++ b/process/main.py @@ -93,6 +93,7 @@ NeutralBeam, ) from process.models.physics.density_limit import PlasmaDensityLimit +from process.models.physics.exhaust import PlasmaExhaust from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.physics import ( DetailedPhysics, @@ -700,12 +701,14 @@ def __init__(self): self.plasma_beta = PlasmaBeta() self.plasma_inductance = PlasmaInductance() self.plasma_density_limit = PlasmaDensityLimit() + self.plasma_exhaust = PlasmaExhaust() self.physics = Physics( plasma_profile=self.plasma_profile, current_drive=self.current_drive, plasma_beta=self.plasma_beta, plasma_inductance=self.plasma_inductance, plasma_density_limit=self.plasma_density_limit, + plasma_exhaust=self.plasma_exhaust, ) self.physics_detailed = DetailedPhysics( plasma_profile=self.plasma_profile, diff --git a/process/models/physics/exhaust.py b/process/models/physics/exhaust.py new file mode 100644 index 000000000..eb0021fb8 --- /dev/null +++ b/process/models/physics/exhaust.py @@ -0,0 +1,120 @@ +import logging + +from process.core import constants + +logger = logging.getLogger(__name__) + + +class PlasmaExhaust: + """Class to hold plasma exhaust calculations for plasma processing.""" + + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE + + @staticmethod + def calculate_separatrix_power( + f_p_alpha_plasma_deposited: float, + p_alpha_total_mw: float, + p_non_alpha_charged_mw: float, + p_hcd_injected_total_mw: float, + p_plasma_ohmic_mw: float, + p_plasma_rad_mw: float, + ) -> float: + """ + Calculate the power crossing the separatrix (P_sep). + + Parameters + ---------- + f_p_alpha_plasma_deposited : float + Fraction of alpha power deposited in plasma. + p_alpha_total_mw : float + Total alpha power produced (MW). + p_non_alpha_charged_mw : float + Power from non-alpha charged particles (MW). + p_hcd_injected_total_mw : float + Total power injected by heating and current drive (MW). + p_plasma_ohmic_mw : float + Ohmic heating power (MW). + p_plasma_rad_mw : float + Radiated power from plasma (MW). + + Returns + ------- + float + Power crossing the separatrix (MW). + """ + + return ( + f_p_alpha_plasma_deposited * p_alpha_total_mw + + p_non_alpha_charged_mw + + p_hcd_injected_total_mw + + p_plasma_ohmic_mw + - p_plasma_rad_mw + ) + + @staticmethod + def calculate_psep_over_r_metric( + p_plasma_separatrix_mw: float, rmajor: float + ) -> float: + """ + Calculate the power crossing the separatrix per unit major radius (P_sep/R). + + Parameters + ---------- + p_plasma_separatrix_mw : float + Power crossing the separatrix (MW). + rmajor : float + Plasma major radius (m). + + Returns + ------- + float + Power crossing the separatrix per unit major radius (MW/m). + """ + return p_plasma_separatrix_mw / rmajor + + @staticmethod + def calculate_eu_demo_re_attachment_metric( + p_plasma_separatrix_mw: float, + b_plasma_toroidal_on_axis: float, + q95: float, + aspect: float, + rmajor: float, + ) -> float: + """Calculate the EU DEMO divertor protection re-attachment metric for plasma exhaust. + + Parameters + ---------- + p_plasma_separatrix_mw : float + Power crossing the separatrix (MW). + b_plasma_toroidal_on_axis : float + Toroidal magnetic field on axis (T). + q95 : float + Safety factor at 95% flux surface. + aspect : float + Aspect ratio of the plasma. + rmajor : float + Plasma major radius (m). + + Returns + ------- + float + EU DEMO re-attachment metric (MW T /m). + + References + ---------- + - M. Siccinio, G. Federici, R. Kembleton, H. Lux, F. Maviglia, and J. Morris, + "Figure of merit for divertor protection in the preliminary design of the EU-DEMO reactor," + Nuclear Fusion, vol. 59, no. 10, pp. 106026-106026, Jul. 2019, + doi: https://doi.org/10.1088/1741-4326/ab3153. + + - H. Zohm et al., + "A stepladder approach to a tokamak fusion power plant," + Nuclear Fusion, vol. 57, no. 8, pp. 086002-086002, May 2017, + doi: https://doi.org/10.1088/1741-4326/aa739e. + """ + + return (p_plasma_separatrix_mw * b_plasma_toroidal_on_axis) / ( + q95 * aspect * rmajor + ) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 02e925af3..5316dd361 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -1632,6 +1632,7 @@ def __init__( plasma_beta, plasma_inductance, plasma_density_limit, + plasma_exhaust, ): self.outfile = constants.NOUT self.mfile = constants.MFILE @@ -1640,6 +1641,7 @@ def __init__( self.beta = plasma_beta self.inductance = plasma_inductance self.density_limit = plasma_density_limit + self.exhaust = plasma_exhaust def physics(self): """Routine to calculate tokamak plasma physics information @@ -2377,12 +2379,31 @@ def physics(self): ) physics_variables.p_plasma_separatrix_mw = ( - physics_variables.f_p_alpha_plasma_deposited - * physics_variables.p_alpha_total_mw - + physics_variables.p_non_alpha_charged_mw - + pinj - + physics_variables.p_plasma_ohmic_mw - - physics_variables.p_plasma_rad_mw + self.exhaust.calculate_separatrix_power( + f_p_alpha_plasma_deposited=physics_variables.f_p_alpha_plasma_deposited, + p_alpha_total_mw=physics_variables.p_alpha_total_mw, + p_non_alpha_charged_mw=physics_variables.p_non_alpha_charged_mw, + p_hcd_injected_total_mw=pinj, + p_plasma_ohmic_mw=physics_variables.p_plasma_ohmic_mw, + p_plasma_rad_mw=physics_variables.p_plasma_rad_mw, + ) + ) + + physics_variables.p_plasma_separatrix_rmajor_mw = ( + self.exhaust.calculate_psep_over_r_metric( + p_plasma_separatrix_mw=physics_variables.p_plasma_separatrix_mw, + rmajor=physics_variables.rmajor, + ) + ) + + physics_variables.p_div_bt_q_aspect_rmajor_mw = ( + self.exhaust.calculate_eu_demo_re_attachment_metric( + p_plasma_separatrix_mw=physics_variables.p_plasma_separatrix_mw, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + q95=physics_variables.q95, + aspect=physics_variables.aspect, + rmajor=physics_variables.rmajor, + ) ) physics_variables.pflux_plasma_surface_neutron_avg_mw = ( @@ -4912,7 +4933,7 @@ def outplas(self): po.oblnkl(self.outfile) po.ovarre( self.outfile, - "Power into divertor zone via charged particles (MW)", + "Plasma separatrix power (MW)", "(p_plasma_separatrix_mw)", physics_variables.p_plasma_separatrix_mw, "OP ", @@ -4939,25 +4960,15 @@ def outplas(self): po.ovarre( self.outfile, "Pdivt / R ratio (MW/m) (On peak divertor)", - "(p_div_separatrix_max_mw/physics_variables.rmajor)", - physics_variables.p_div_separatrix_max_mw / physics_variables.rmajor, + "(p_plasma_separatrix_rmajor_mw)", + physics_variables.p_plasma_separatrix_rmajor_mw, "OP ", ) po.ovarre( self.outfile, "Pdivt Bt / qAR ratio (MWT/m) (On peak divertor)", - "(pdivmaxbt/qar)", - ( - ( - physics_variables.p_div_separatrix_max_mw - * physics_variables.b_plasma_toroidal_on_axis - ) - / ( - physics_variables.q95 - * physics_variables.aspect - * physics_variables.rmajor - ) - ), + "(p_div_bt_q_aspect_rmajor_mw)", + physics_variables.p_div_bt_q_aspect_rmajor_mw, "OP ", ) else: @@ -4965,25 +4976,15 @@ def outplas(self): po.ovarre( self.outfile, "Psep / R ratio (MW/m)", - "(p_plasma_separatrix_mw/rmajor)", - physics_variables.p_plasma_separatrix_mw / physics_variables.rmajor, + "(p_plasma_separatrix_rmajor_mw)", + physics_variables.p_plasma_separatrix_rmajor_mw, "OP ", ) po.ovarre( self.outfile, "Psep Bt / qAR ratio (MWT/m)", - "(pdivtbt_over_qar)", - ( - ( - physics_variables.p_plasma_separatrix_mw - * physics_variables.b_plasma_toroidal_on_axis - ) - / ( - physics_variables.q95 - * physics_variables.aspect - * physics_variables.rmajor - ) - ), + "(p_div_bt_q_aspect_rmajor_mw)", + physics_variables.p_div_bt_q_aspect_rmajor_mw, "OP ", ) diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 72be4b201..7722b2f47 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -20,6 +20,7 @@ NeutralBeam, ) from process.models.physics.density_limit import PlasmaDensityLimit +from process.models.physics.exhaust import PlasmaExhaust from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.physics import ( DetailedPhysics, @@ -58,6 +59,7 @@ def physics(): PlasmaBeta(), PlasmaInductance(), PlasmaDensityLimit(), + PlasmaExhaust(), ) diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 3809b0fa8..0fa5e0444 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -30,6 +30,7 @@ NeutralBeam, ) from process.models.physics.density_limit import PlasmaDensityLimit +from process.models.physics.exhaust import PlasmaExhaust from process.models.physics.physics import ( Physics, PlasmaBeta, @@ -88,6 +89,7 @@ def stellarator(): PlasmaBeta(), PlasmaInductance(), PlasmaDensityLimit(), + PlasmaExhaust(), ), Neoclassics(), plasma_beta=PlasmaBeta(),