diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index 40a1b1139..b2329f45c 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 4c6557296..136502830 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,8 @@ 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_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 8165ede28..554b9b85b 100644 --- a/process/main.py +++ b/process/main.py @@ -92,6 +92,7 @@ LowerHybrid, NeutralBeam, ) +from process.models.physics.exhaust import PlasmaExhaust from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.physics import ( DetailedPhysics, @@ -698,11 +699,13 @@ def __init__(self): ) self.plasma_beta = PlasmaBeta() self.plasma_inductance = PlasmaInductance() + 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_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 66d06eb33..ac3202531 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -1627,13 +1627,21 @@ def _trapped_particle_fraction_sauter( class Physics: - def __init__(self, plasma_profile, current_drive, plasma_beta, plasma_inductance): + def __init__( + self, + plasma_profile, + current_drive, + plasma_beta, + plasma_inductance, + plasma_exhaust, + ): self.outfile = constants.NOUT self.mfile = constants.MFILE self.plasma_profile = plasma_profile self.current_drive = current_drive self.beta = plasma_beta self.inductance = plasma_inductance + self.exhaust = plasma_exhaust def physics(self): """Routine to calculate tokamak plasma physics information @@ -2403,12 +2411,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.plfux_plasma_surface_neutron_avg_mw = ( @@ -5308,7 +5335,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 ", @@ -5335,25 +5362,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: @@ -5361,25 +5378,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 956912aa7..39af3e01f 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -19,6 +19,7 @@ LowerHybrid, NeutralBeam, ) +from process.models.physics.exhaust import PlasmaExhaust from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.physics import ( DetailedPhysics, @@ -56,6 +57,7 @@ def physics(): ), PlasmaBeta(), PlasmaInductance(), + PlasmaExhaust(), ) diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 268a94b67..0d76a13b0 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -29,7 +29,12 @@ LowerHybrid, NeutralBeam, ) -from process.models.physics.physics import Physics, PlasmaBeta, PlasmaInductance +from process.models.physics.exhaust import PlasmaExhaust +from process.models.physics.physics import ( + Physics, + PlasmaBeta, + PlasmaInductance, +) from process.models.physics.plasma_profiles import PlasmaProfile from process.models.power import Power from process.models.stellarator.build import st_build @@ -82,6 +87,7 @@ def stellarator(): ), PlasmaBeta(), PlasmaInductance(), + PlasmaExhaust(), ), Neoclassics(), plasma_beta=PlasmaBeta(),