feat: Describing fluxes & inductances of the plasma

cfspopcon/algorithms/__init__.py

@@ -7,9 +7,11 @@
 from .composite_algorithm import predictive_popcon
 from .core_radiated_power import calc_core_radiated_power
 from .extrinsic_core_radiator import calc_extrinsic_core_radiator
+from .fluxes import calc_fluxes
 from .fusion_gain import calc_fusion_gain
 from .geometry import calc_geometry
 from .heat_exhaust import calc_heat_exhaust
+from .inductances import calc_inductances
 from .ohmic_power import calc_ohmic_power
 from .peaked_profiles import calc_peaked_profiles
 from .plasma_current_from_q_star import calc_plasma_current_from_q_star
@@ -24,6 +26,8 @@
 
 ALGORITHMS: dict[Algorithms, Union[Algorithm, CompositeAlgorithm]] = {
     Algorithms["calc_beta"]: calc_beta,
+    Algorithms["calc_inductances"]: calc_inductances,
+    Algorithms["calc_fluxes"]: calc_fluxes,
     Algorithms["calc_core_radiated_power"]: calc_core_radiated_power,
     Algorithms["calc_extrinsic_core_radiator"]: calc_extrinsic_core_radiator,
     Algorithms["calc_fusion_gain"]: calc_fusion_gain,
@@ -58,6 +62,8 @@ def get_algorithm(algorithm: Union[Algorithms, str]) -> Union[Algorithm, Composi
     "calc_extrinsic_core_radiator",
     "calc_fusion_gain",
     "calc_geometry",
+    "calc_fluxes",
+    "calc_inductances",
     "calc_heat_exhaust",
     "calc_ohmic_power",
     "calc_peaked_profiles",

cfspopcon/algorithms/fluxes.py

@@ -0,0 +1,63 @@
+"""Calculate PF flux contribution and resistive, internal, and external flux consumed over the ramp-up."""
+from .. import formulas
+from ..unit_handling import Unitfull, convert_to_default_units
+from .algorithm_class import Algorithm
+
+RETURN_KEYS = [
+    "internal_flux",
+    "external_flux",
+    "resistive_flux",
+    "poloidal_field_flux",
+    "max_flux_for_flattop",
+    "max_flattop_duration",
+    "breakdown_flux_consumption",
+    "flux_needed_from_CS_over_rampup",
+]
+
+
+def run_calc_fluxes(
+    plasma_current: Unitfull,
+    major_radius: Unitfull,
+    internal_inductance: Unitfull,
+    external_inductance: Unitfull,
+    ejima_coefficient: Unitfull,
+    vertical_field_mutual_inductance: Unitfull,
+    vertical_magnetic_field: Unitfull,
+    loop_voltage: Unitfull,
+    total_flux_available_from_CS: Unitfull,
+) -> dict[str, Unitfull]:
+    """Calculate PF flux contribution and resistive, internal, and external flux consumed over the ramp-up.
+
+    Args:
+        plasma_current: :term:`glossary link<plasma_current>`
+        major_radius: :term:`glossary link<major_radius>`
+        internal_inductance: :term:`glossary link<internal_inductance>`
+        external_inductance: :term:`glossary link<external_inductance>`
+        vertical_field_mutual_inductance: :term:`glossary link<vertical_field_mutual_inductance>`
+        vertical_magnetic_field: :term:`glossary link<vertical_magnetic_field>`
+        ejima_coefficient: :term:`glossary link<ejima_coefficient>`
+        loop_voltage: :term:`glossary link<loop_voltage>`
+        total_flux_available_from_CS: :term:`glossary link<total_flux_available_from_CS>`
+
+    Returns:
+        :term:`resistive_flux`, :term:`internal_flux`, :term:`external_flux`, :term:`max_flattop_duration`, :term:`max_flux_for_flattop`, :term:`breakdown_flux_consumption`, :term:`glossary link<flux_needed_from_CS_over_rampup>`
+    """
+    internal_flux = formulas.calc_flux_internal(plasma_current, internal_inductance)
+    external_flux = formulas.calc_flux_external(plasma_current, external_inductance)
+    resistive_flux = formulas.calc_flux_res(plasma_current, major_radius, ejima_coefficient)
+    poloidal_field_flux = formulas.calc_flux_PF(vertical_field_mutual_inductance, vertical_magnetic_field, major_radius)
+
+    flux_needed_from_CS_over_rampup = internal_flux + external_flux + resistive_flux - poloidal_field_flux
+    max_flux_for_flattop = total_flux_available_from_CS - internal_flux - external_flux - resistive_flux + poloidal_field_flux
+    max_flattop_duration = max_flux_for_flattop / loop_voltage
+
+    breakdown_flux_consumption = formulas.calc_breakdown_flux_consumption(major_radius)
+
+    local_vars = locals()
+    return {key: convert_to_default_units(local_vars[key], key) for key in RETURN_KEYS}
+
+
+calc_fluxes = Algorithm(
+    function=run_calc_fluxes,
+    return_keys=RETURN_KEYS,
+)

cfspopcon/algorithms/geometry.py

@@ -10,6 +10,7 @@
     "vertical_minor_radius",
     "plasma_volume",
     "surface_area",
+    "poloidal_circumference",
 ]
 
 
@@ -32,7 +33,7 @@ def run_calc_geometry(
         triangularity_ratio_sep_to_psi95: :term:`glossary link<triangularity_ratio_sep_to_psi95>`
 
     Returns:
-        :term:`separatrix_elongation`, :term:`separatrix_triangularity`, :term:`minor_radius`, :term:`vertical_minor_radius`, :term:`plasma_volume`, :term:`surface_area`
+        :term:`separatrix_elongation`, :term:`separatrix_triangularity`, :term:`minor_radius`, :term:`vertical_minor_radius`, :term:`plasma_volume`, :term:`surface_area`, :term:`poloidal_circumference`
     """
     separatrix_elongation = areal_elongation * elongation_ratio_sep_to_areal
 
@@ -43,6 +44,7 @@ def run_calc_geometry(
 
     plasma_volume = formulas.calc_plasma_volume(major_radius, inverse_aspect_ratio, areal_elongation)
     surface_area = formulas.calc_plasma_surface_area(major_radius, inverse_aspect_ratio, areal_elongation)
+    poloidal_circumference = formulas.calc_plasma_poloidal_circumference(minor_radius, areal_elongation)
 
     local_vars = locals()
     return {key: convert_to_default_units(local_vars[key], key) for key in RETURN_KEYS}

cfspopcon/algorithms/inductances.py

@@ -0,0 +1,99 @@
+"""Calculate the vertical magnetic field, as well as the plasma surface's mutual inductance with the vertical field, internal inductivity, external inductance and internal inductance."""
+from .. import formulas, named_options
+from ..unit_handling import Unitfull, convert_to_default_units
+from .algorithm_class import Algorithm
+
+RETURN_KEYS = [
+    "internal_inductivity",
+    "internal_inductance",
+    "external_inductance",
+    "vertical_field_mutual_inductance",
+    "vertical_magnetic_field",
+]
+
+
+def run_calc_inductances(
+    major_radius: Unitfull,
+    plasma_volume: Unitfull,
+    poloidal_circumference: Unitfull,
+    internal_inductance_geometry: named_options.InternalInductanceGeometry,
+    plasma_current: Unitfull,
+    magnetic_field_on_axis: Unitfull,
+    minor_radius: Unitfull,
+    safety_factor_on_axis: Unitfull,
+    inverse_aspect_ratio: Unitfull,
+    areal_elongation: Unitfull,
+    beta_poloidal: Unitfull,
+    vertical_magnetic_field_equation: named_options.VertMagneticFieldEq,
+    surface_inductance_coefficients: named_options.SurfaceInductanceCoeffs,
+    internal_inductivity: Unitfull = None,
+) -> dict[str, Unitfull]:
+    """Calculate the vertical magnetic field, as well as the plasma surface's mutual inductance with the vertical field, internal inductivity, external inductance and internal inductance.
+
+    Args:
+        major_radius: :term:`glossary link<major_radius>`
+        plasma_volume: [m**3] :term:`glossary<plasma_volume>`
+        poloidal_circumference: [m] :term:`glossary<poloidal_circumference>`
+        internal_inductance_geometry: [~] :term:`glossary<internal_inductance_geometry>`
+        plasma_current: :term:`glossary link<plasma_current>`
+        magnetic_field_on_axis: [T] :term:`glossary<magnetic_field_on_axis>`
+        minor_radius: :term:`glossary link<minor_radius>`
+        safety_factor_on_axis: [~] :term:`glossary<safety_factor_on_axis>`
+        inverse_aspect_ratio: [~] :term:`glossary link<inverse_aspect_ratio>`
+        areal_elongation: [~] :term:`glossary<areal_elongation>`
+        beta_poloidal: [~] :term:`glossary link<beta_poloidal>`
+        vertical_magnetic_field_equation: [~] :term:`glossary link<vertical_magnetic_field_equation>`
+        surface_inductance_coefficients: [~] :term:`glossary link<surface_inductance_coefficients>`
+        internal_inductivity: [~] :term:`glossary<internal_inductivity>`
+
+    Returns:
+        :term:`internal_inductivity`,
+        :term:`internal_inductance`,
+        :term:`external_inductance`,
+        :term:`vertical_field_mutual_inductance`,
+        :term:`vertical_magnetic_field`
+    """
+    if internal_inductivity is None:
+        internal_inductivity = formulas.calc_internal_inductivity(
+            plasma_current, major_radius, magnetic_field_on_axis, minor_radius, safety_factor_on_axis
+        )
+
+    internal_inductance = formulas.calc_internal_inductance(
+        major_radius, internal_inductivity, plasma_volume, poloidal_circumference, internal_inductance_geometry
+    )
+    external_inductance = formulas.calc_external_inductance(
+        inverse_aspect_ratio, areal_elongation, beta_poloidal, major_radius, internal_inductivity, surface_inductance_coefficients
+    )
+    vertical_field_mutual_inductance = formulas.calc_vertical_field_mutual_inductance(
+        inverse_aspect_ratio, areal_elongation, surface_inductance_coefficients
+    )
+    invmu_0_dLedR = formulas.inductances.calc_invmu_0_dLedR(
+        inverse_aspect_ratio,
+        areal_elongation,
+        beta_poloidal,
+        internal_inductivity,
+        external_inductance,
+        major_radius,
+        surface_inductance_coefficients,
+    )
+    vertical_magnetic_field = formulas.calc_vertical_magnetic_field(
+        inverse_aspect_ratio,
+        areal_elongation,
+        beta_poloidal,
+        internal_inductivity,
+        external_inductance,
+        major_radius,
+        plasma_current,
+        invmu_0_dLedR,
+        vertical_magnetic_field_equation,
+        surface_inductance_coefficients,
+    )
+
+    local_vars = locals()
+    return {key: convert_to_default_units(local_vars[key], key) for key in RETURN_KEYS}
+
+
+calc_inductances = Algorithm(
+    function=run_calc_inductances,
+    return_keys=RETURN_KEYS,
+)

cfspopcon/formulas/__init__.py

@@ -24,9 +24,17 @@
 from .divertor_metrics import calc_B_pol_omp, calc_B_tor_omp
 from .energy_confinement_time_scalings import calc_tau_e_and_P_in_from_scaling
 from .figures_of_merit import calc_normalised_collisionality, calc_peak_pressure, calc_rho_star, calc_triple_product
+from .fluxes import calc_breakdown_flux_consumption, calc_flux_external, calc_flux_internal, calc_flux_PF, calc_flux_res
 from .fusion_rates import calc_fusion_power, calc_neutron_flux_to_walls
-from .geometry import calc_plasma_surface_area, calc_plasma_volume
+from .geometry import calc_plasma_poloidal_circumference, calc_plasma_surface_area, calc_plasma_volume
 from .impurity_effects import calc_change_in_dilution, calc_change_in_zeff, calc_impurity_charge_state
+from .inductances import (
+    calc_external_inductance,
+    calc_internal_inductance,
+    calc_internal_inductivity,
+    calc_vertical_field_mutual_inductance,
+    calc_vertical_magnetic_field,
+)
 from .operational_limits import calc_greenwald_density_limit, calc_greenwald_fraction, calc_troyon_limit
 from .plasma_profiles import calc_1D_plasma_profiles
 from .Q_thermal_gain_factor import thermal_calc_gain_factor
@@ -58,8 +66,19 @@
     "calc_bootstrap_fraction",
     "calc_current_relaxation_time",
     "calc_f_shaping",
+    "calc_plasma_poloidal_circumference",
+    "calc_flux_external",
+    "calc_flux_internal",
+    "calc_flux_res",
+    "calc_flux_PF",
+    "calc_breakdown_flux_consumption",
     "calc_fusion_power",
     "calc_greenwald_fraction",
+    "calc_internal_inductivity",
+    "calc_internal_inductance",
+    "calc_external_inductance",
+    "calc_vertical_field_mutual_inductance",
+    "calc_vertical_magnetic_field",
     "calc_LH_transition_threshold_power",
     "calc_LI_transition_threshold_power",
     "calc_loop_voltage",

cfspopcon/formulas/fluxes.py

@@ -0,0 +1,109 @@
+"""Resistive flux consumption, inductive flux consumption (internally and externally on the plasma surface) during purely ohmic ramp-up."""
+import numpy as np
+from scipy import constants  # type: ignore[import]
+
+from ..unit_handling import Unitfull, ureg, wraps_ufunc
+
+
+@wraps_ufunc(
+    input_units=dict(plasma_current=ureg.A, major_radius=ureg.m, ejima_coefficient=ureg.dimensionless),
+    return_units=dict(flux_res=ureg.weber),
+)
+def calc_flux_res(plasma_current: float, major_radius: float, ejima_coefficient: float = 0.4) -> float:
+    """Calculate the resistive flux.
+
+    Chapter 8: Plasma operation and control: Physics cite:`Gribov_2007`
+    NOTE: CE, the Ejima coefficient, is chosen by default to be 0.4, See for example...
+    Chapter 8: Plasma operation and control: cite:`Gribov_2007`
+    Ohmic flux consumption during initial operation of the NSTX spherical torus :cite:`Menard_2001`
+
+    Args:
+        plasma_current: [A] :term:`glossary link<plasma_current>`
+        major_radius: [m] :term:`glossary link<major_radius>`
+        ejima_coefficient: [~] :term:`glossary link<ejima_coefficient>`
+
+
+    Returns:
+        [weber] :term:`resistive_flux`
+    """
+    return float(ejima_coefficient * constants.mu_0 * plasma_current * major_radius)
+
+
+### COMPONENTS OF INDUCTIVE FLUX: INTERNAL AND EXTERNAL ###
+
+
+@wraps_ufunc(input_units=dict(plasma_current=ureg.A, internal_inductance=ureg.henry), return_units=dict(internal_flux=ureg.weber))
+def calc_flux_internal(plasma_current: float, internal_inductance: float) -> Unitfull:
+    """Calculate the flux due to the plasma current and internal inductance of the plasma (assuming a circular cross-section).
+
+    From: A power-balance model for local helicity injection startup in a spherical tokamak :cite:`Barr_2018`
+    NOTE: This is (plasma current times) equation 25 from Barr but applied to a plasma with a
+    circular cross-section and a non-cirular cross-section.
+
+    Args:
+        plasma_current: [A] :term:`glossary link<plasma_current>`
+        internal_inductance: [henry] :term:`glossary link<internal_inductance>`
+
+    Returns:
+        [weber] :term:`internal_flux`
+    """
+    internal_flux = plasma_current * internal_inductance
+
+    return float(internal_flux)
+
+
+@wraps_ufunc(
+    input_units=dict(
+        vertical_field_mutual_inductance=ureg.dimensionless,
+        vertical_magnetic_field=ureg.T,
+        major_radius=ureg.m,
+    ),
+    return_units=dict(flux_PF=ureg.weber),
+)
+def calc_flux_PF(vertical_field_mutual_inductance: float, vertical_magnetic_field: float, major_radius: float) -> Unitfull:
+    """Calculate the surface flux contribution from the vertical magnetic field required for radial force balance (which arises from the poloidal field coils).
+
+    From: A power-balance model for local helicity injection startup in a spherical tokamak :cite:`Barr_2018`
+
+    Args:
+        vertical_field_mutual_inductance: [~] :term:`glossary link<vertical_field_mutual_inductance>`
+        vertical_magnetic_field: [henry] :term:`glossary link<vertical_magnetic_field>`
+        major_radius: [m] :term:`glossary link<major_radius>`
+
+    Returns:
+        [weber] :term:`poloidal_field_flux`
+    """
+    return float(np.pi * major_radius**2 * vertical_field_mutual_inductance * vertical_magnetic_field)
+
+
+@wraps_ufunc(input_units=dict(plasma_current=ureg.A, external_inductance=ureg.henry), return_units=dict(external_flux=ureg.weber))
+def calc_flux_external(plasma_current: float, external_inductance: float) -> Unitfull:
+    """Calculate the surface flux generated by the plasma current.
+
+    From: A power-balance model for local helicity injection startup in a spherical tokamak :cite:`Barr_2018`
+
+    Args:
+        plasma_current: [A] :term:`glossary link<plasma_current>`
+        external_inductance: [henry] :term:`glossary link<external_inductance>`
+
+    Returns:
+        [weber] :term:`external_flux`
+    """
+    return float(plasma_current * external_inductance)
+
+
+@wraps_ufunc(input_units=dict(major_radius=ureg.m), return_units=dict(breakdown_flux_consumption=ureg.weber))
+def calc_breakdown_flux_consumption(major_radius: float) -> Unitfull:
+    """Calculate the resistive flux required for breakdown.
+
+    Plasma Design Considerations of Near Term Tokamak Fusion Experimental Reactor :cite:`Sugihara`
+    NOTE: given the way the ejima_coefficient is emprically derived in :cite:`Gribov_2007` (i.e., with an implicit assumption that the ramp is defined from Ip=0)
+    there is reason to believe a seperate calculation for flux consumed over breakdown is not necessary, but it is included here anyways.
+
+    Args:
+        major_radius: [m] :term:`glossary link<major_radius>`
+
+    Returns:
+        [weber] :term:`breakdown_flux_consumption`
+    """
+    return float(0.073 * major_radius - 0.00665)