diff --git a/documentation/source/physics-models/plasma_density.md b/documentation/source/physics-models/plasma_density.md index 5c4fa4ab3..5b07350b9 100644 --- a/documentation/source/physics-models/plasma_density.md +++ b/documentation/source/physics-models/plasma_density.md @@ -1,7 +1,7 @@ -# Density Limit +# Density Limit | `PlasmaDensityLimit` Several density limit models are available in PROCESS. These are -calculated in routine `calculate_density_limit()`, which is called by `physics`. +calculated in the run routine `density_limit.run()`, which is called by `physics()`. This constraint can be activated by stating `icc = 5` in the input file. @@ -13,7 +13,7 @@ For the models below $P_{\perp}$ is the mean heat flux density across the separa ----------------- -## ASDEX model +## ASDEX model | `calculate_asdex_density_limit()` Switch value: `i_density_limit = 1`[^1][^2] @@ -23,7 +23,7 @@ $$ ----------------- -## Borrass model for ITER, I +## Borrass model for ITER, I | `calculate_borrass_iter_i_density_limit()` Switch value: `i_density_limit = 2` [^1] @@ -35,7 +35,7 @@ $C \approx$ 1.8 for ITER-like conditions. ----------------- -## Borrass model for ITER, II +## Borrass model for ITER, II | `calculate_borrass_iter_ii_density_limit()` Switch value: `i_density_limit = 3` [^1] @@ -45,7 +45,7 @@ $$ ----------------- -## JET edge radiation model +## JET edge radiation model | `calculate_jet_edge_radiation_density_limit()` Switch value: `i_density_limit = 4` [^1] @@ -55,7 +55,7 @@ $$ ----------------- -## JET simplified model +## JET simplified model | `calculate_jet_simple_density_limit()` Switch value: `i_density_limit = 5` [^1] @@ -73,7 +73,7 @@ where $\kappa \approx 1.5, \Delta \approx 0.1a$ has been taken from JET. ----------------- -## Hugill-Murakami model +## Hugill-Murakami model | `calculate_hugill_murakami_density_limit()` Switch value: `i_density_limit = 6` [^2] @@ -84,7 +84,7 @@ $$ ----------------- -## Greenwald model +## Greenwald model | `calculate_greenwald_density_limit()` Switch value: `i_density_limit = 7` [^3][^4] @@ -96,7 +96,7 @@ For the Greenwald model the limit applies to the line-averaged electron density, --------------------- -## ASDEX New model +## ASDEX New model | `calculate_asdex_new_density_limit()` Switch value: `i_density_limit = 8` [^5][^6] diff --git a/process/main.py b/process/main.py index 8165ede28..7388e750c 100644 --- a/process/main.py +++ b/process/main.py @@ -92,6 +92,7 @@ LowerHybrid, NeutralBeam, ) +from process.models.physics.density_limit import PlasmaDensityLimit 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_density_limit = PlasmaDensityLimit() 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, ) self.physics_detailed = DetailedPhysics( plasma_profile=self.plasma_profile, diff --git a/process/models/physics/density_limit.py b/process/models/physics/density_limit.py new file mode 100644 index 000000000..0819ff76f --- /dev/null +++ b/process/models/physics/density_limit.py @@ -0,0 +1,645 @@ +import logging +from enum import IntEnum + +import numpy as np + +from process.core import constants +from process.core import process_output as po +from process.core.exceptions import ProcessValueError +from process.data_structure import ( + current_drive_variables, + divertor_variables, + physics_variables, +) + +logger = logging.getLogger(__name__) + + +class DensityLimitModel(IntEnum): + """Electron density model types""" + + ASDEX = 1 + BORRASS_ITER_I = 2 + BORRASS_ITER_II = 3 + JET_EDGE_RADIATION = 4 + JET_SIMPLE = 5 + HUGILL_MURAKAMI = 6 + GREENWALD = 7 + ASDEX_NEW = 8 + + +class PlasmaDensityLimit: + """Class to hold plasma density limit calculations for plasma processing.""" + + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE + + def run(self): + physics_variables.nd_plasma_electron_max_array, _ = self.calculate_density_limit( + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + i_density_limit=physics_variables.i_density_limit, + p_plasma_separatrix_mw=physics_variables.p_plasma_separatrix_mw, + p_hcd_injected_total_mw=current_drive_variables.p_hcd_injected_total_mw, + plasma_current=physics_variables.plasma_current, + prn1=divertor_variables.prn1, + qcyl=physics_variables.qstar, + q95=physics_variables.q95, + rmajor=physics_variables.rmajor, + rminor=physics_variables.rminor, + a_plasma_surface=physics_variables.a_plasma_surface, + zeff=physics_variables.n_charge_plasma_effective_vol_avg, + ) + + # Calculate beta_norm_max based on i_beta_norm_max + try: + model = DensityLimitModel(int(physics_variables.i_density_limit)) + physics_variables.nd_plasma_electrons_max = self.get_density_limit_value( + model + ) + except ValueError: + raise ProcessValueError( + "Illegal value of i_density_limit", + i_density_limit=physics_variables.i_density_limit, + ) from None + + def get_density_limit_value(self, model: DensityLimitModel) -> float: + """ + Get the density limit value (n_e_max) for the specified model. + + Parameters + ---------- + model : DensityLimitModel + The density limit model type. + + Returns + ------- + float + The density limit value (m⁻³). + """ + model_map = { + DensityLimitModel.ASDEX: physics_variables.nd_plasma_electron_max_array[0], + DensityLimitModel.BORRASS_ITER_I: physics_variables.nd_plasma_electron_max_array[ + 1 + ], + DensityLimitModel.BORRASS_ITER_II: physics_variables.nd_plasma_electron_max_array[ + 2 + ], + DensityLimitModel.JET_EDGE_RADIATION: physics_variables.nd_plasma_electron_max_array[ + 3 + ], + DensityLimitModel.JET_SIMPLE: physics_variables.nd_plasma_electron_max_array[ + 4 + ], + DensityLimitModel.HUGILL_MURAKAMI: physics_variables.nd_plasma_electron_max_array[ + 5 + ], + DensityLimitModel.GREENWALD: physics_variables.nd_plasma_electron_max_array[ + 6 + ], + DensityLimitModel.ASDEX_NEW: physics_variables.nd_plasma_electron_max_array[ + 7 + ], + } + return model_map[model] + + @staticmethod + def calculate_asdex_density_limit( + p_perp: float, + b_plasma_toroidal_on_axis: float, + q95: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the ASDEX density limit. + + Parameters + ---------- + p_perp : float + Perpendicular power density (MW/m²). + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + q95 : float + Safety factor at 95% of the plasma poloidal flux. + rmajor : float + Plasma major radius (m). + prn1 : float + Edge density / average plasma density. + + Returns + ------- + float + The ASDEX density limit (m⁻³). + + References + ---------- + T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + return ( + 1.54e20 + * p_perp**0.43 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.45 + ) / prn1 + + @staticmethod + def calculate_borrass_iter_i_density_limit( + p_perp: float, + b_plasma_toroidal_on_axis: float, + q95: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the Borrass ITER I density limit. + + Parameters + ---------- + p_perp : float + Perpendicular power density (MW/m²). + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + q95 : float + Safety factor at 95% of the plasma poloidal flux. + rmajor : float + Plasma major radius (m). + prn1 : float + Edge density / average plasma density. + + Returns + ------- + float + The Borrass ITER I density limit (m⁻³). + + References + ---------- + T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + return ( + 1.8e20 + * p_perp**0.53 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.22 + ) / prn1 + + @staticmethod + def calculate_borrass_iter_ii_density_limit( + p_perp: float, + b_plasma_toroidal_on_axis: float, + q95: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the Borrass ITER II density limit. + + Parameters + ---------- + p_perp : float + Perpendicular power density (MW/m²). + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + q95 : float + Safety factor at 95% of the plasma poloidal flux. + rmajor : float + Plasma major radius (m). + prn1 : float + Edge density / average plasma density. + + Returns + ------- + float + The Borrass ITER II density limit (m⁻³). + + References + ---------- + T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + return ( + 0.5e20 + * p_perp**0.57 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.09 + ) / prn1 + + @staticmethod + def calculate_jet_edge_radiation_density_limit( + zeff: float, p_hcd_injected_total_mw: float, prn1: float, qcyl: float + ) -> float: + """ + Calculate the JET edge radiation density limit. + + Parameters + ---------- + zeff : float + Effective charge (Z_eff). + p_hcd_injected_total_mw : float + Power injected into the plasma (MW). + prn1 : float + Edge density / average plasma density. + qcyl : float + Equivalent cylindrical safety factor (qstar). + + Returns + ------- + float + The JET edge radiation density limit (m⁻³). + + References + ---------- + T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + + denom = (zeff - 1.0) * (1.0 - 4.0 / (3.0 * qcyl)) + if denom <= 0.0: + return 0.0 + return (1.0e20 * np.sqrt(p_hcd_injected_total_mw / denom)) / prn1 + + @staticmethod + def calculate_jet_simple_density_limit( + b_plasma_toroidal_on_axis: float, + p_plasma_separatrix_mw: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the JET simple density limit. + + Parameters + ---------- + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + p_plasma_separatrix_mw : float + Power crossing the separatrix (MW). + rmajor : float + Plasma major radius (m). + prn1 : float + Edge density / average plasma density. + + Returns + ------- + float + The JET simple density limit (m⁻³). + + References + ---------- + T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + return ( + 0.237e20 + * b_plasma_toroidal_on_axis + * np.sqrt(p_plasma_separatrix_mw) + / rmajor + ) / prn1 + + @staticmethod + def calculate_hugill_murakami_density_limit( + b_plasma_toroidal_on_axis: float, rmajor: float, qcyl: float + ) -> float: + """ + Calculate the Hugill-Murakami density limit. + + Parameters + ---------- + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + rmajor : float + Plasma major radius (m). + qcyl : float + Equivalent cylindrical safety factor (qstar). + + Returns + ------- + float + The Hugill-Murakami density limit (m⁻³). + + References + ---------- + N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989' + """ + + return 3.0e20 * b_plasma_toroidal_on_axis / (rmajor * qcyl) + + @staticmethod + def calculate_greenwald_density_limit(c_plasma: float, rminor: float) -> float: + """ + Calculate the Greenwald density limit (n_GW). + + Parameters + ---------- + c_plasma : float + Plasma current (A). + rminor : float + Plasma minor radius (m). + + Returns + ------- + float + The Greenwald density limit (m⁻³). + + Notes + ----- + The Greenwald limit is typically applied to the line averaged electron density. + + References + ---------- + M. Greenwald et al., "A new look at density limits in tokamaks," + Nuclear Fusion, vol. 28, no. 12, pp. 2199-2207, Dec. 1988, + doi: https://doi.org/10.1088/0029-5515/28/12/009. + + M. Greenwald, "Density limits in toroidal plasmas," + Plasma Physics and Controlled Fusion, vol. 44, no. 8, pp. R27-R53, Jul. 2002, + doi: https://doi.org/10.1088/0741-3335/44/8/201. + """ + + return 1.0e14 * c_plasma / (np.pi * rminor**2) + + @staticmethod + def calculate_asdex_new_density_limit( + p_hcd_injected_total_mw: float, c_plasma: float, q95: float, prn1: float + ) -> float: + """ + Calculate the ASDEX Upgrade new density limit. + + Parameters + ---------- + p_hcd_injected_total_mw : float + Power injected into the plasma (MW). + c_plasma : float + Plasma current (A). + q95 : float + Safety factor at 95% surface. + prn1 : float + Edge density / average plasma density. + + Returns + ------- + float + The ASDEX Upgrade new density limit (m⁻³). + + Notes + ----- + This limit is for the separatrix density so we scale by `prn1` to get it as a volume average. + + References + ---------- + J. W. Berkery et al., "Density limits as disruption forecasters for spherical tokamaks," + Plasma Physics and Controlled Fusion, vol. 65, no. 9, pp. 095003-095003, Jul. 2023, + doi: https://doi.org/10.1088/1361-6587/ace476. + + M. Bernert et al., "The H-mode density limit in the full tungsten ASDEX Upgrade tokamak," vol. 57, no. 1, pp. 014038-014038, Nov. 2014, + doi: https://doi.org/10.1088/0741-3335/57/1/014038. + """ + return ( + 1.0e20 + * 0.506 + * (p_hcd_injected_total_mw**0.396 * (c_plasma / 1.0e6) ** 0.265) + / (q95**0.323) + ) / prn1 + + def calculate_density_limit( + self, + b_plasma_toroidal_on_axis: float, + i_density_limit: int, + p_plasma_separatrix_mw: float, + p_hcd_injected_total_mw: float, + plasma_current: float, + prn1: float, + qcyl: float, + q95: float, + rmajor: float, + rminor: float, + a_plasma_surface: float, + zeff: float, + ) -> tuple[np.ndarray, float]: + """ + Calculate the density limit using various models. + + Parameters + ---------- + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + i_density_limit : int + Switch denoting which formula to enforce (1-7). + p_plasma_separatrix_mw : float + Power flowing to the edge plasma via charged particles (MW). + p_hcd_injected_total_mw : float + Power injected into the plasma (MW). + plasma_current : float + Plasma current (A). + prn1 : float + Edge density / average plasma density. + qcyl : float + Equivalent cylindrical safety factor (qstar). + q95 : float + Safety factor at 95% surface. + rmajor : float + Plasma major radius (m). + rminor : float + Plasma minor radius (m). + a_plasma_surface : float + Plasma surface area (m²). + zeff : float + Plasma effective charge. + + Returns + ------- + tuple[np.ndarray, float] + A tuple containing: + - nd_plasma_electron_max_array : Average plasma density limit using eight different models (m⁻³). + - nd_plasma_electrons_max : Enforced average plasma density limit (m⁻³). + + Raises + ------ + ValueError + If i_density_limit is not between 1 and 7. + + Notes + ----- + This routine calculates several different formulae for the density limit and enforces the one chosen by the user. + For i_density_limit = 1-5, 8, we scale the separatrix density limit output by the ratio of the separatrix to volume averaged density. + + References + ---------- + AEA FUS 172: Physics Assessment for the European Reactor Study + + N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989' + + M. Bernert et al., "The H-mode density limit in the full tungsten ASDEX Upgrade tokamak," + vol. 57, no. 1, pp. 014038-014038, Nov. 2014, doi: https://doi.org/10.1088/0741-3335/57/1/014038. + """ + + if i_density_limit < 1 or i_density_limit > 7: + raise ProcessValueError( + "Illegal value for i_density_limit", i_density_limit=i_density_limit + ) + + nd_plasma_electron_max_array = np.empty((8,)) + + # Power per unit area crossing the plasma edge + # (excludes radiation and neutrons) + + p_perp = p_plasma_separatrix_mw / a_plasma_surface + + # Old ASDEX density limit formula + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + + nd_plasma_electron_max_array[0] = self.calculate_asdex_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + q95=q95, + rmajor=rmajor, + prn1=prn1, + ) + + # Borrass density limit model for ITER (I) + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # Borrass et al, ITER-TN-PH-9-6 (1989) + + nd_plasma_electron_max_array[1] = self.calculate_borrass_iter_i_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + q95=q95, + rmajor=rmajor, + prn1=prn1, + ) + + # Borrass density limit model for ITER (II) + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # This formula is (almost) identical to that in the original routine + # denlim (now deleted). + + nd_plasma_electron_max_array[2] = self.calculate_borrass_iter_ii_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + q95=q95, + rmajor=rmajor, + prn1=prn1, + ) + + # JET edge radiation density limit model + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # qcyl=qstar here, but literature is not clear. + + nd_plasma_electron_max_array[3] = ( + self.calculate_jet_edge_radiation_density_limit( + zeff=zeff, + p_hcd_injected_total_mw=p_hcd_injected_total_mw, + prn1=prn1, + qcyl=qcyl, + ) + ) + + # JET simplified density limit model + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + + nd_plasma_electron_max_array[4] = self.calculate_jet_simple_density_limit( + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + p_plasma_separatrix_mw=p_plasma_separatrix_mw, + rmajor=rmajor, + prn1=prn1, + ) + + # Hugill-Murakami M.q limit + # qcyl=qstar here, which is okay according to the literature + + nd_plasma_electron_max_array[5] = self.calculate_hugill_murakami_density_limit( + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, rmajor=rmajor, qcyl=qcyl + ) + + # Greenwald limit + + nd_plasma_electron_max_array[6] = self.calculate_greenwald_density_limit( + c_plasma=plasma_current, rminor=rminor + ) + + nd_plasma_electron_max_array[7] = self.calculate_asdex_new_density_limit( + p_hcd_injected_total_mw=p_hcd_injected_total_mw, + c_plasma=plasma_current, + q95=q95, + prn1=prn1, + ) + + # Enforce the chosen density limit + + return nd_plasma_electron_max_array, nd_plasma_electron_max_array[ + i_density_limit - 1 + ] + + def output_density_limit_information(self): + """Output density limit information to file.""" + + po.osubhd(self.outfile, "Density Limit using different models :") + po.ovarre( + self.outfile, + "Old ASDEX model", + "(nd_plasma_electron_max_array(1))", + physics_variables.nd_plasma_electron_max_array[0], + "OP ", + ) + po.ovarre( + self.outfile, + "Borrass ITER model I", + "(nd_plasma_electron_max_array(2))", + physics_variables.nd_plasma_electron_max_array[1], + "OP ", + ) + po.ovarre( + self.outfile, + "Borrass ITER model II", + "(nd_plasma_electron_max_array(3))", + physics_variables.nd_plasma_electron_max_array[2], + "OP ", + ) + po.ovarre( + self.outfile, + "JET edge radiation model", + "(nd_plasma_electron_max_array(4))", + physics_variables.nd_plasma_electron_max_array[3], + "OP ", + ) + po.ovarre( + self.outfile, + "JET simplified model", + "(nd_plasma_electron_max_array(5))", + physics_variables.nd_plasma_electron_max_array[4], + "OP ", + ) + po.ovarre( + self.outfile, + "Hugill-Murakami Mq model", + "(nd_plasma_electron_max_array(6))", + physics_variables.nd_plasma_electron_max_array[5], + "OP ", + ) + po.ovarre( + self.outfile, + "Greenwald model", + "(nd_plasma_electron_max_array(7))", + physics_variables.nd_plasma_electron_max_array[6], + "OP ", + ) + po.ovarre( + self.outfile, + "ASDEX New", + "(nd_plasma_electron_max_array(8))", + physics_variables.nd_plasma_electron_max_array[7], + "OP ", + ) + po.ovarre( + self.outfile, + "Density limit from scaling (/m3)", + "(nd_plasma_electrons_max)", + physics_variables.nd_plasma_electrons_max, + "OP ", + ) + + po.oblnkl(self.outfile) + po.ostars(self.outfile, 110) + po.oblnkl(self.outfile) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index d2405c263..02e925af3 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -1625,13 +1625,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_density_limit, + ): 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.density_limit = plasma_density_limit def physics(self): """Routine to calculate tokamak plasma physics information @@ -2412,24 +2420,7 @@ def physics(self): physics_variables.kappa95, ) - # Density limit - ( - physics_variables.nd_plasma_electron_max_array, - physics_variables.nd_plasma_electrons_max, - ) = self.calculate_density_limit( - physics_variables.b_plasma_toroidal_on_axis, - physics_variables.i_density_limit, - physics_variables.p_plasma_separatrix_mw, - current_drive_variables.p_hcd_injected_total_mw, - physics_variables.plasma_current, - divertor_variables.prn1, - physics_variables.qstar, - physics_variables.q95, - physics_variables.rmajor, - physics_variables.rminor, - physics_variables.a_plasma_surface, - physics_variables.n_charge_plasma_effective_vol_avg, - ) + self.density_limit.run() # Calculate transport losses and energy confinement time using the # chosen scaling law @@ -2697,192 +2688,6 @@ def calculate_current_profile_index_wesson(qstar: float, q0: float) -> float: """ return qstar / q0 - 1.0 - @staticmethod - def calculate_density_limit( - b_plasma_toroidal_on_axis: float, - i_density_limit: int, - p_plasma_separatrix_mw: float, - p_hcd_injected_total_mw: float, - plasma_current: float, - prn1: float, - qcyl: float, - q95: float, - rmajor: float, - rminor: float, - a_plasma_surface: float, - zeff: float, - ) -> tuple[np.ndarray, float]: - """Calculate the density limit using various models. - - Parameters - ---------- - b_plasma_toroidal_on_axis : float - Toroidal field on axis (T). - i_density_limit : int - Switch denoting which formula to enforce. - p_plasma_separatrix_mw : float - Power flowing to the edge plasma via charged particles (MW). - p_hcd_injected_total_mw : float - Power injected into the plasma (MW). - plasma_current : float - Plasma current (A). - prn1 : float - Edge density / average plasma density. - qcyl : float - Equivalent cylindrical safety factor (qstar). - q95 : float - Safety factor at 95% surface. - rmajor : float - Plasma major radius (m). - rminor : float - Plasma minor radius (m). - a_plasma_surface : float - Plasma surface area (m^2). - zeff : float - Plasma effective charge. - - Returns - ------- - Tuple[np.ndarray, float] - A tuple containing: - - nd_plasma_electron_max_array (np.ndarray): Average plasma density limit using seven different models (m^-3). - - nd_plasma_electrons_max (float): Enforced average plasma density limit (m^-3). - - Raises - ------ - ValueError - If i_density_limit is not between 1 and 7. - Notes - ----- - ValueError - If i_density_limit is not between 1 and 7. - Notes - ----- - This routine calculates several different formulae for the density limit and enforces the one chosen by the user. - ValueError - If i_density_limit is not between 1 and 7. - Notes - ----- - This routine calculates several different formulae for the density limit and enforces the one chosen by the user. - For i_density_limit = 1-5, 8, we scale the sepatrix density limit output by the ratio of the separatrix to volume averaged density - - References - ---------- - - AEA FUS 172 - Physics Assessment for the European Reactor Study - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines - 1989 - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines - 1989 - - M. Bernert et al., “The H-mode density limit in the full tungsten ASDEX Upgrade tokamak,” - vol. 57, no. 1, pp. 014038-014038, Nov. 2014, doi: https://doi.org/10.1088/0741-3335/57/1/014038. ‌ - - """ - - if i_density_limit < 1 or i_density_limit > 7: - raise ProcessValueError( - "Illegal value for i_density_limit", i_density_limit=i_density_limit - ) - - nd_plasma_electron_max_array = np.empty((8,)) - - # Power per unit area crossing the plasma edge - # (excludes radiation and neutrons) - - p_perp = p_plasma_separatrix_mw / a_plasma_surface - - # Old ASDEX density limit formula - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - - nd_plasma_electron_max_array[0] = ( - 1.54e20 - * p_perp**0.43 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.45 - ) / prn1 - - # Borrass density limit model for ITER (I) - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # Borrass et al, ITER-TN-PH-9-6 (1989) - - nd_plasma_electron_max_array[1] = ( - 1.8e20 - * p_perp**0.53 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.22 - ) / prn1 - - # Borrass density limit model for ITER (II) - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # This formula is (almost) identical to that in the original routine - # denlim (now deleted). - - nd_plasma_electron_max_array[2] = ( - 0.5e20 - * p_perp**0.57 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.09 - ) / prn1 - - # JET edge radiation density limit model - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # qcyl=qstar here, but literature is not clear. - - denom = (zeff - 1.0) * (1.0 - 4.0 / (3.0 * qcyl)) - if denom <= 0.0: - if i_density_limit == 4: - logger.error( - f"qcyl < 4/3; nd_plasma_electron_max_array(4) set to zero; model 5 will be enforced instead. {denom=} {qcyl=}" - ) - i_density_limit = 5 - - nd_plasma_electron_max_array[3] = 0.0 - else: - nd_plasma_electron_max_array[3] = ( - 1.0e20 * np.sqrt(p_hcd_injected_total_mw / denom) - ) / prn1 - - # JET simplified density limit model - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - - nd_plasma_electron_max_array[4] = ( - 0.237e20 - * b_plasma_toroidal_on_axis - * np.sqrt(p_plasma_separatrix_mw) - / rmajor - ) / prn1 - - # Hugill-Murakami M.q limit - # qcyl=qstar here, which is okay according to the literature - - nd_plasma_electron_max_array[5] = ( - 3.0e20 * b_plasma_toroidal_on_axis / (rmajor * qcyl) - ) - - # Greenwald limit - - nd_plasma_electron_max_array[6] = ( - 1.0e14 * plasma_current / (np.pi * rminor * rminor) - ) - - nd_plasma_electron_max_array[7] = ( - 1.0e20 - * 0.506 - * (p_hcd_injected_total_mw**0.396 * (plasma_current / 1.0e6) ** 0.265) - / (q95**0.323) - ) / prn1 - - # Enforce the chosen density limit - - return nd_plasma_electron_max_array, nd_plasma_electron_max_array[ - i_density_limit - 1 - ] - @staticmethod def plasma_composition(): """Calculates various plasma component fractional makeups. @@ -4586,74 +4391,7 @@ def outplas(self): po.oblnkl(self.outfile) if stellarator_variables.istell == 0: - po.osubhd(self.outfile, "Density Limit using different models :") - po.ovarre( - self.outfile, - "Old ASDEX model", - "(nd_plasma_electron_max_array(1))", - physics_variables.nd_plasma_electron_max_array[0], - "OP ", - ) - po.ovarre( - self.outfile, - "Borrass ITER model I", - "(nd_plasma_electron_max_array(2))", - physics_variables.nd_plasma_electron_max_array[1], - "OP ", - ) - po.ovarre( - self.outfile, - "Borrass ITER model II", - "(nd_plasma_electron_max_array(3))", - physics_variables.nd_plasma_electron_max_array[2], - "OP ", - ) - po.ovarre( - self.outfile, - "JET edge radiation model", - "(nd_plasma_electron_max_array(4))", - physics_variables.nd_plasma_electron_max_array[3], - "OP ", - ) - po.ovarre( - self.outfile, - "JET simplified model", - "(nd_plasma_electron_max_array(5))", - physics_variables.nd_plasma_electron_max_array[4], - "OP ", - ) - po.ovarre( - self.outfile, - "Hugill-Murakami Mq model", - "(nd_plasma_electron_max_array(6))", - physics_variables.nd_plasma_electron_max_array[5], - "OP ", - ) - po.ovarre( - self.outfile, - "Greenwald model", - "(nd_plasma_electron_max_array(7))", - physics_variables.nd_plasma_electron_max_array[6], - "OP ", - ) - po.ovarre( - self.outfile, - "ASDEX New", - "(nd_plasma_electron_max_array(8))", - physics_variables.nd_plasma_electron_max_array[7], - "OP ", - ) - po.ovarre( - self.outfile, - "Density limit from scaling (/m3)", - "(nd_plasma_electrons_max)", - physics_variables.nd_plasma_electrons_max, - "OP ", - ) - - po.oblnkl(self.outfile) - po.ostars(self.outfile, 110) - po.oblnkl(self.outfile) + self.density_limit.output_density_limit_information() po.oheadr(self.outfile, "Plasma Reactions :") diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 956912aa7..72be4b201 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -19,6 +19,7 @@ LowerHybrid, NeutralBeam, ) +from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.physics import ( DetailedPhysics, @@ -56,6 +57,7 @@ def physics(): ), PlasmaBeta(), PlasmaInductance(), + PlasmaDensityLimit(), ) @@ -2473,7 +2475,7 @@ def test_calculate_density_limit(calculatedensitylimitparam, physics): """ nd_plasma_electron_max_array, nd_plasma_electrons_max = ( - physics.calculate_density_limit( + physics.density_limit.calculate_density_limit( i_density_limit=calculatedensitylimitparam.i_density_limit, b_plasma_toroidal_on_axis=calculatedensitylimitparam.b_plasma_toroidal_on_axis, p_plasma_separatrix_mw=calculatedensitylimitparam.p_plasma_separatrix_mw, diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 268a94b67..3809b0fa8 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.density_limit import PlasmaDensityLimit +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(), + PlasmaDensityLimit(), ), Neoclassics(), plasma_beta=PlasmaBeta(),