Create plasma exhaust class

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(

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}",

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

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,

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
+        )

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,51 +4960,31 @@ 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:
             # Single null divertor configuration
             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 ",
             )
 

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(),
     )
 
 

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(),