♻️ Refactor CS coil

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

@@ -33,7 +33,7 @@ This method calculates the CS geometry parameters. The CS is assumed to be a per
 2. The half height of the CS is set relative to that of the inside height of the TF and can be scaled by changing the input value of `f_z_cs_tf_internal`:
 
     $$
-    \overbrace{z_{\text{CS,half}}}^{\texttt{z_cs_inside_half}} = \overbrace{z_{\text{TF,inside-half}}}^{\texttt{z_tf_inside_half}} \times \texttt{f_z_cs_tf_internal}
+    \overbrace{z_{\text{CS,half}}}^{\texttt{z_cs_half}} = \overbrace{z_{\text{TF,inside-half}}}^{\texttt{z_tf_inside_half}} \times \texttt{f_z_cs_tf_internal}
     $$
 
 3. The full height of the CS is thus simply given by
@@ -42,17 +42,91 @@ This method calculates the CS geometry parameters. The CS is assumed to be a per
     \overbrace{dz_{\text{CS}}}^{\texttt{dz_cs_full}} = z_{\text{CS,half}} \times 2
     $$
 
+    Along with the upper and lower dimensions:
+
+    $$
+    z_{\text{CS,upper}} = z_{\text{CS,half}}
+    $$
+
+    $$
+    z_{\text{CS,lower}} = -z_{\text{CS,upper}}
+    $$
+
+
 4. The outboard edge of the CS is given by:
 
 
     $$
     r_{\text{CS,outer}} = r_{\text{CS,middle}} + \frac{dr_{\text{CS}}}{2}
     $$
 
-5. The full poloidal cross-sectional area is given by:
+5. The inboard edge of the CS is given by:
+
+    $$
+    r_{\text{CS,inner}} = r_{\text{CS,outer}} - dr_{\text{CS}}
+    $$
+
+6. The full radial width of the CS is given by:
+
+    $$
+    dr_{\text{CS,full}}  = 2 \times r_{\text{CS,outer}}
+    $$
+
+
+7. The full poloidal cross-sectional area is given by:
 
     $$
-    \overbrace{A_{\text{CS,poloidal}}}^{\texttt{a_cs_poloidal}} = 2 \times dr_{\text{CS}} \times dz_{\text{CS}}
+    \overbrace{A_{\text{CS,poloidal}}}^{\texttt{a_cs_poloidal}} = dr_{\text{CS}} \times dz_{\text{CS}}
+    $$
+
+8. The full top-down toroidal cross-sectional area is given by:
+
+    $$
+    \overbrace{A_{\text{CS,toroidal}}}^{\texttt{a_cs_toroidal}} = \pi \left(r_{\text{CS,outer}}^2 - r_{\text{CS,inner}}^2 \right)
+    $$
+
+------------
+
+### EU-DEMO Turn Geometry | `calculate_cs_turn_geometry_eu_demo()`
+
+This superconducting turn structure for the CS assumes a rectangular turn shape with a "stadium"-shaped cable area[^eu_demo_turn].
+
+
+![CS turn layout](../eng-models/images/cs_eu_demo_turn.PNG "CS EU-DEMO like turn")
+
+The turn geometry is calculated as follows:
+
+1. The vertical height of the turn is given by:
+
+    $$
+    dz_{\text{CS,turn}} = \sqrt{\frac{A_{\text{CS,turn}}}{\texttt{f_dr_dz_cs_turn}}}
+    $$
+
+    $\texttt{f_dr_dz_cs_turn}$ is the intended length to height ratio of the turn
+
+2. The radial width or length of the turn is now given by:
+
+    $$
+    dr_{\text{CS,turn}} = \texttt{f_dr_dz_cs_turn} \times dz_{\text{CS,turn}}
+    $$
+
+3. The radius of the corners of the cable space is given by:
+
+    $$
+    r_{\text{CS,cable space corner}} = - \frac{dr_{\text{CS,turn}}-dz_{\text{CS,turn}}}{\pi} \\
+    + \sqrt{\left(\frac{dr_{\text{CS,turn}}-dz_{\text{CS,turn}}}{\pi}\right)^2+ \frac{dr_{\text{CS,turn}}dz_{\text{CS,turn}}(4-\pi)r_{\text{CS,turn corners}}^2 - (A_{\text{CS,turn}}\times \texttt{f_a_cs_turn_steel})}{\pi}}
+    $$
+
+4. The thickness of the conduit around the cable space is given by:
+
+    $$
+    dz_{\text{CS,turn conduit}} = \frac{dz_{\text{CS,turn}}}{2} - r_{\text{CS,cable space corner}}
+    $$
+
+    In this model the vertical and radial conduit thicknesses are equal so:
+
+    $$
+    dr_{\text{CS,turn conduit}}  = dz_{\text{CS,turn conduit}}
     $$
 
 
@@ -66,21 +140,32 @@ 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.
+
+This is Equation 3.13 from "Case Studies in Superconducting Magnets"[^2].
+
+-------------
 
-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:
+### 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 we cannot 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 +208,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 +233,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)}
+$$
+
+-------------
+
+
+### Hoop stress | `calculate_cs_hoop_stress()`
+
+
 
-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:
 
+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)
 $$
-f_z = \frac{f_V}{2}.
+
+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:
+
 $$
+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:
@@ -320,4 +463,5 @@ constraints (26 and 27) are activated.
 
 [^1]: M. N. Wilson, Superconducting Magnets. Oxford University Press, USA, 1983, ISBN 13: 9780198548102
 [^2]: Case Studies in Superconducting Magnets. Boston, MA: Springer US, 2009. doi: https://doi.org/10.1007/b112047.
+[^eu_demo_turn]: R. Wesche et al., “Central solenoid winding pack design for DEMO,” Fusion Engineering and Design, vol. 124, pp. 82-85, Apr. 2017, doi: https://doi.org/10.1016/j.fusengdes.2017.04.052.
 ‌

process/core/input.py

@@ -798,6 +798,9 @@ def __post_init__(self):
     "dx_tf_wp_insertion_gap": InputVariable("tfcoil", float, range=(1e-10, 0.1)),
     "f_dr_tf_outboard_inboard": InputVariable("build", float, range=(0.2, 5.0)),
     "tftmp": InputVariable("tfcoil", 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("tfcoil", float, range=(0.0, 1.0)),
     "theta1_coil": InputVariable("tfcoil", float, range=(0.1, 60.0)),

process/core/io/plot/summary.py

@@ -24,6 +24,7 @@
 from process.data_structure.numerics import FiguresOfMerit, PROCESSRunMode
 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,
@@ -10050,6 +10052,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"
@@ -10131,14 +10134,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)
 
 
@@ -15513,14 +15516,32 @@ def main_plot(
 
     plot_current_profiles_over_time(figs[30].add_subplot(111), m_file, scan)
 
-    plot_cs_stress_time_profile(axis=figs[31].add_subplot(311), mfile=m_file, scan=scan)
+    plot_cs_stress_time_profile(axis=figs[31].add_subplot(431), mfile=m_file, scan=scan)
+
+    cs_coil = CSCoil(cs_fatigue=CsFatigue())
+    cs_coil.plot_cs_radial_hoop_stress_profile(
+        axis=figs[31].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[31].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[31].add_subplot(223, aspect="equal"), figs[31], m_file, scan
     )
     plot_cs_turn_structure(
         figs[31].add_subplot(326, aspect="equal"), figs[31], m_file, scan
     )
+    figs[31].subplots_adjust(wspace=0.3)
 
     plot_first_wall_top_down_cross_section(
         figs[32].add_subplot(221, aspect="equal"), m_file, scan

process/data_structure/build_variables.py

@@ -138,10 +138,10 @@ class BuildData:
     """
 
     dr_cs: float = 0.811
-    """Central solenoid thickness (m) (`iteration variable 16`)"""
+    """Central solenoid radial thickness (m) (`iteration variable 16`)"""
 
     dr_cs_precomp: float = 0.0
-    """CS coil precompression structure thickness (m)"""
+    """CS coil precompression structure radial thickness (m)"""
 
     rbld: float = 0.0
     """sum of thicknesses to the major radius (m)"""