♻️ Refactor CS coil

documentation/source/eng-models/central-solenoid.md

@@ -66,21 +66,27 @@ This method calculates the CS geometry parameters. The CS is assumed to be a per
 
 -----------
 
-### Self peak magnetic field | `calculate_cs_self_peak_magnetic_field()`
+### Self peak bore magnetic field | `calculate_cs_bore_magnetic_field()`
 
-The general form for the field at the very centre of the central solenoid bore with uniform current density and rectangular cross-section is given by:
+The self induced peak field at the central bore of a homogenous current density rectangular cross-section solenoid is given by:
 
 $$
-B_0 = J_{\text{CS}}aF(\alpha,\beta)
+B_0 = J_{\text{CS}}r_{\text{CS,inner}}F(\alpha,\beta)
 $$
 
 $$
 F(\alpha,\beta) = \mu_0\beta \ln{\left[\frac{\alpha+\sqrt{\alpha^2+\beta^2}}{1+\sqrt{1+\beta^2}}\right]}
 $$
 
-where $\alpha = \frac{r_{\text{CS,outer}}}{r_{\text{CS,inner}}}$, is the ratio of the outer and inner radii of the solenoid and $\beta = \frac{z_{\text{CS,half}}}{r_{\text{CS,outer}}}$, is the ratio of the solenoid half height to its inboard radius.
+where $\alpha = \frac{r_{\text{CS,outer}}}{r_{\text{CS,inner}}}$, is the ratio of the outer and inner radii of the solenoid and $\beta = \frac{z_{\text{CS,half}}}{r_{\text{CS,inner}}}$, is the ratio of the solenoid half height to its inboard radius.
 
-The peak field at the bore of the central solenoid will not be the same as that felt by the conductors inside the structures. We require to know the peak field on the conductor if we are to design a superconducting central solenoid that has enough margin. Fits to data[^1] for different ranges of $\beta$ have been calculated as follows:
+This is Equation 3.13 from "Case Studies in Superconducting Magnets"[^2].
+
+-------------
+
+### Self peak on coil magnetic field | `calculate_cs_self_peak_magnetic_field()`
+
+The peak field at the bore of the central solenoid will not be the same as that felt by the conductors inside the structures.So wecannot use the bore value directly calculated by [`calculate_cs_bore_magnetic_field()`](#self-peak-bore-magnetic-field--calculate_cs_bore_magnetic_field). We require to know the peak field on the conductor if we are to design a superconducting central solenoid that has enough margin. Fits to data[^1] for different ranges of $\beta$ have been calculated as follows to scale the bore field value by:
 
 - $\beta > 3.0$
 
@@ -123,10 +129,10 @@ The peak field at the bore of the central solenoid will not be the same as that
 
 -----------
 
-### Axial stresses | `calculate_cs_self_peak_midplane_axial_stress()`
+### Peak Axial Stresses | `calculate_cs_self_peak_midplane_axial_stress()`
 
 
-The vertical (axial) force for a "thin-walled" solenoid ($\alpha = 1$) at the midplane is given by[^2]:
+The vertical (axial) force for a "thin-walled" solenoid ($\alpha = 1$) at the midplane where the force is maximum is given by Equation 3.41 in "Case studies in superconducting magnets" [^2]:
 
 $$
 F_{z}(0)=\frac{\mu_0}{2}\left(\frac{N I}{2 \times dz_{\text{half}}}\right) \times \\
@@ -148,41 +154,99 @@ $$
 
 Here $K(k)$ and $E(k)$ are the complete elliptic integrals, respectively of the first and second kinds.
 
+!!! info "Non thin-walled solenoids"
+
+    For solenoids that can be classed as "thick-walled" ($\alpha > 1$), the coil has to be divided radially into several thin-walled solenoids. This is not currently performed though more info can be found in Section 3.5.5 of "Case studies in superconducting magnets" [^2].
 
 The axial compressive force at $z$ in an isolated solenoid increases from 0 at $z  = dz_{\text{half}}$
 to the maximum at the midplane, $F_{z}(0)$.
 
+The axial stress in the steel is given by:
+
+$$
+\sigma_z = \frac{F_z}{f_z A_z}
+$$
+
+where $F_z$ is the axial force, $f_z$ is the fraction of the horizontal cross-section occupied by 
+steel, and $A_z$ is the area of the horizontal cross-section.
+
+The fraction of the horizontal cross-section occupied by steel is calculated assuming that the 
+conductor is square and has a steel jacket with the same thickness on all four sides, giving:
+
+$$
+f_z = 0.5.
+$$
+
 
 
 --------------------------
 
 
-### Hoop stress | `hoop_stress()`
+### Radial stress | `calculate_cs_radial_stress()`
 
 
+The self induced radial stress is calculated using Equation 4.11 from "Superconducting magnets" [^1].
 
+$$
+\sigma_{r} = \frac{K(2+v)}{3(\alpha+1)}\times \left(\alpha^2+\alpha+1-\frac{\alpha^2}{\epsilon^2}-\epsilon(\alpha+1)\right) \\
+- \frac{M(3+v)}{8}\left(\alpha^2+1-\frac{\alpha^2}{\epsilon^2}-\epsilon^2\right)
+$$
 
-The hoop stress is calculated using equations 4.10 and 4.11 from "Superconducting magnets", Martin N. 
-Wilson (1983).  This is divided by the fraction of the area occupied by steel to obtain the hoop 
-stress in the steel, $\sigma_{hoop}$.
+Where:
 
-The axial stress can be calculated using "Case studies in superconducting magnets", Y. Iwasa, p. 
-86, 3.5.2, Special Case 4: Midplane force.  This applies exactly only to a thin-walled solenoid. 
-The axial stress in the steel is given by:
+- $\epsilon = \frac{r}{r_{\text{CS,inner}}}$
+- $\alpha = \frac{r_{\text{CS,outer}}}{r_{\text{CS,inner}}}$
+
+The terms $K$ and $M$ are given by:
 
 $$
-\sigma_z = \frac{F_z}{f_z A_z}
+K = \frac{J r_{\text{CS,inner}}\left(\alpha B_{\text{CS,inner}} - B_{\text{CS,outer}}\right)}{(\alpha-1)}
 $$
 
-where $F_z$ is the axial force, $f_z$ is the fraction of the horizontal cross-section occupied by 
-steel, and $A_z$ is the area of the horizontal cross-section.
+$$
+M = \frac{J r_{\text{CS,inner}}\left(B_{\text{CS,inner}} - B_{\text{CS,outer}}\right)}{(\alpha-1)}
+$$
 
-The fraction of the horizontal cross-section occupied by steel is calculated assuming that the 
-conductor is square and has a steel jacket with the same thickness on all four sides, giving:
+-------------
+
+
+### Hoop stress | `calculate_cs_hoop_stress()`
+
+
+
+
+The hoop stress is calculated using Equation 4.10 from "Superconducting magnets" [^1]. This is divided by the fraction of the area occupied by steel to obtain the hoop 
+stress in the steel, $\sigma_{\theta}$.
+
+$$
+\sigma_{\theta} = \frac{K(2+v)}{3(\alpha+1)}\times \left(\alpha^2+\alpha+1+\frac{\alpha^2}{\epsilon^2}-\epsilon \frac{(1+2v)(\alpha+1)}{(2+v)}\right) \\
+- \frac{M(3+v)}{8}\left(\alpha^2+1+\frac{\alpha^2}{\epsilon^2}-\frac{(1+3v)}{(3+v)}\epsilon^2\right)
+$$
+
+Where:
+
+- $\epsilon = \frac{r}{r_{\text{CS,inner}}}$
+- $\alpha = \frac{r_{\text{CS,outer}}}{r_{\text{CS,inner}}}$
+
+The terms $K$ and $M$ are given by:
 
 $$
-f_z = \frac{f_V}{2}.
+K = \frac{J r_{\text{CS,inner}}\left(\alpha B_{\text{CS,inner}} - B_{\text{CS,outer}}\right)}{(\alpha-1)}
+$$
+
 $$
+M = \frac{J r_{\text{CS,inner}}\left(B_{\text{CS,inner}} - B_{\text{CS,outer}}\right)}{(\alpha-1)}
+$$
+
+!!! warning "Assumption of outer field, $B_{\text{CS,outer}}$" 
+
+    In this case we currently assume that $B_{\text{CS,outer}} = 0$. This is the same as that for an infinite solenoid. Approximation of the outboard field is currently not performed.
+
+--------------------------
+
+
+
+
 
 The radial stress is neglected. The hoop and axial stresses are combined to give the maximum shear 
 stress, as required by the Tresca stress criterion:

process/core/input.py

@@ -972,6 +972,9 @@ def __post_init__(self):
     ),
     "f_dr_tf_outboard_inboard": InputVariable("build", float, range=(0.2, 5.0)),
     "tftmp": InputVariable(data_structure.tfcoil_variables, float, range=(0.01, 293.0)),
+    "temp_cs_superconductor_operating": InputVariable(
+        data_structure.pfcoil_variables, float, range=(0.01, 293.0)
+    ),
     "tgain": InputVariable("ife", float, range=(1.0, 500.0)),
     "th_joint_contact": InputVariable(
         data_structure.tfcoil_variables, float, range=(0.0, 1.0)

process/core/io/plot/summary.py

@@ -24,6 +24,7 @@
 from process.data_structure import impurity_radiation_module
 from process.data_structure.pfcoil_variables import NFIXMX
 from process.models.build import Build
+from process.models.cs_fatigue import CsFatigue
 from process.models.geometry.blanket import (
     blanket_geometry_double_null,
     blanket_geometry_single_null,
@@ -47,6 +48,7 @@
     vacuum_vessel_geometry_double_null,
     vacuum_vessel_geometry_single_null,
 )
+from process.models.pfcoil import CSCoil
 from process.models.physics.bootstrap_current import BootstrapCurrentFractionModel
 from process.models.physics.confinement_time import (
     ConfinementTimeModel,
@@ -9754,6 +9756,7 @@ def plot_cs_coil_structure(
         f"CS full height: {mfile.get('dz_cs_full', scan=scan):.4f} m\n"
         f"CS full width: {mfile.get('dr_cs_full', scan=scan):.4f} m\n"
         f"CS poloidal area: {mfile.get('a_cs_poloidal', scan=scan):.4f} m$^2$\n"
+        f"CS top-down toroidal area: {mfile.get('a_cs_toroidal', scan=scan):.4f} m$^2$\n"
         f"$N_{{\\text{{turns}}}}:$ {mfile.get('n_pf_coil_turns[n_cs_pf_coils-1]', scan=scan):,.2f}\n"
         f"$I_{{\\text{{peak}}}}:$ {mfile.get('c_pf_cs_coils_peak_ma[n_cs_pf_coils-1]', scan=scan):.3f}$ \\ MA$\n"
         f"$B_{{\\text{{peak}}}}:$ {mfile.get('b_pf_coil_peak[n_cs_pf_coils-1]', scan=scan):.3f}$ \\ T$\n"
@@ -9835,14 +9838,14 @@ def plot_cs_stress_time_profile(axis: plt.Axes, mfile: MFile, scan: int) -> None
         stress_z_cs_self_midplane_profile / 1e6,
         "o-",
         linewidth=2,
-        markersize=8,
-        label="Midplane Axial Stress",
+        markersize=4,
+        label="$\\sigma_{z}$,Axial Stress",
     )
     axis.set_xlabel("Pulse Time (s)")
-    axis.set_ylabel("Stress (MPa)")
+    axis.set_ylabel("Midplane Axial Stress (MPa)")
     axis.minorticks_on()
     axis.legend(loc="best")
-    axis.set_title("Central Solenoid Stress")
+    axis.set_title("CS Midplane Axial Stress Time Profile")
     axis.grid(True, alpha=0.3)
 
 
@@ -14496,14 +14499,32 @@ def main_plot(
 
     plot_current_profiles_over_time(figs[29].add_subplot(111), m_file, scan)
 
-    plot_cs_stress_time_profile(axis=figs[30].add_subplot(311), mfile=m_file, scan=scan)
+    plot_cs_stress_time_profile(axis=figs[30].add_subplot(431), mfile=m_file, scan=scan)
+
+    cs_coil = CSCoil(cs_fatigue=CsFatigue())
+    cs_coil.plot_cs_radial_hoop_stress_profile(
+        axis=figs[30].add_subplot(432),
+        mfile=m_file,
+        scan=scan,
+        j_cs=m_file.get("j_cs_pulse_start", scan=scan),
+        b_cs_inner=m_file.get("b_cs_peak_pulse_start", scan=scan),
+    )
+
+    cs_coil.plot_cs_radial_stress_profile(
+        axis=figs[30].add_subplot(433),
+        mfile=m_file,
+        scan=scan,
+        j_cs=m_file.get("j_cs_pulse_start", scan=scan),
+        b_cs_inner=m_file.get("b_cs_peak_pulse_start", scan=scan),
+    )
 
     plot_cs_coil_structure(
         figs[30].add_subplot(223, aspect="equal"), figs[30], m_file, scan
     )
     plot_cs_turn_structure(
         figs[30].add_subplot(326, aspect="equal"), figs[30], m_file, scan
     )
+    figs[30].subplots_adjust(wspace=0.3)
 
     plot_first_wall_top_down_cross_section(
         figs[31].add_subplot(221, aspect="equal"), m_file, scan

process/data_structure/pfcoil_variables.py

@@ -42,7 +42,13 @@ class PFCoilData:
     )
     """Axial stress (z) in central solenoid at midplane due to its own field at each time point (Pa)"""
 
-    sig_hoop: float = 0.0
+    stress_hoop_cs_inner: float = 0.0
+    """Hoop stress in central solenoid at inboard edge due to its own field (when at peak current) (Pa)"""
+
+    stress_hoop_cs_inner_profile: list[float] = field(
+        default_factory=lambda: np.zeros(6)
+    )
+    """Hoop stress in central solenoid at inboard edge due to its own field at each time point (Pa)"""
 
     forc_z_cs_self_peak_midplane: float = 0.0
     """Axial force (z) on central solenoid at midplane due to its own field (when at peak current) (N)"""
@@ -100,21 +106,33 @@ class PFCoilData:
     - =1 Hoop + Axial stress
     """
 
+    a_cs_toroidal: float = 0.0
+    """Central solenoid top-down toroidal cross-sectional area (m2)"""
+
     a_cs_poloidal: float = 0.0
     """Central solenoid vertical cross-sectional area (m2)"""
 
+    a_cs_steel_poloidal: float = 0.0
+    """Central solenoid vertical cross-sectional area of steel (m2)"""
+
     a_cs_turn: float = 0.0
     """Central solenoid (OH) trun cross-sectional area (m2)"""
 
-    awpoh: float = 0.0
+    a_cs_cable_space: float = 0.0
     """central solenoid conductor+void area with area of steel subtracted (m2)"""
 
+    b_cs_self_outer_midplane: float = 0.0
+    """magnetic field at outer edge of central solenoid at midplane due to its own current (T)"""
+
     b_cs_peak_flat_top_end: float = 0.0
     """maximum field in central solenoid at end of flat-top (EoF) (T)"""
 
     b_cs_peak_pulse_start: float = 0.0
     """maximum field in central solenoid at beginning of pulse (T)"""
 
+    b_cs_bore_peak: float = 0.0
+    """peak field in central solenoid bore (T)"""
+
     b_pf_coil_peak: list[float] = field(default_factory=lambda: np.zeros(NGC2))
     """peak field at coil i (T)"""
 
@@ -347,6 +365,12 @@ class PFCoilData:
     r_cs_middle: float = 0.0
     """radius to the centre of the central solenoid (m)"""
 
+    r_cs_inner: float = 0.0
+    """inner radius of the central solenoid (m)"""
+
+    r_cs_outer: float = 0.0
+    """outer radius of the central solenoid (m)"""
+
     dz_cs_full: float = 0.0
     """Full height of the central solenoid (m)"""
 
@@ -395,6 +419,9 @@ class PFCoilData:
     temp_cs_superconductor_margin: float = 0.0
     """Central solenoid temperature margin (K)"""
 
+    temp_cs_superconductor_operating: float = 4.75
+    """Central solenoid operating temperature (K)"""
+
     n_pf_coil_turns: list[float] = field(default_factory=lambda: np.zeros(NGC2))
     """number of turns in PF coil i"""
 

process/models/costs/costs.py

@@ -1716,7 +1716,7 @@ def acc2222(self):
                 if self.data.pf_coil.i_pf_conductor == 0:
                     costpfsc = (
                         self.data.costs.ucsc[self.data.pf_coil.i_cs_superconductor - 1]
-                        * self.data.pf_coil.awpoh
+                        * self.data.pf_coil.a_cs_cable_space
                         * (1 - self.data.pf_coil.f_a_cs_void)
                         * (1 - self.data.pf_coil.fcuohsu)
                         / self.data.pf_coil.n_pf_coil_turns[
@@ -1746,7 +1746,7 @@ def acc2222(self):
             if self.data.pf_coil.i_pf_conductor == 0:
                 costpfcu = (
                     self.data.costs.uccu
-                    * self.data.pf_coil.awpoh
+                    * self.data.pf_coil.a_cs_cable_space
                     * (1 - self.data.pf_coil.f_a_cs_void)
                     * self.data.pf_coil.fcuohsu
                     / self.data.pf_coil.n_pf_coil_turns[
@@ -1758,7 +1758,7 @@ def acc2222(self):
                 # MDK I don't know if this is ccorrect as we never use the resistive model
                 costpfcu = (
                     self.data.costs.uccu
-                    * self.data.pf_coil.awpoh
+                    * self.data.pf_coil.a_cs_cable_space
                     * (1 - self.data.pf_coil.f_a_cs_void)
                     / self.data.pf_coil.n_pf_coil_turns[
                         self.data.pf_coil.n_cs_pf_coils - 1