Refactor CS coil stress routines and variables (#3962)

chris-ashe · web-flow · commit 5eeaeff8145c · 2025-11-04T15:33:58.000Z
* Refactor axial stress calculation in CSCoil class to improve parameter handling and documentation

* 🔄 - Rename axial_stress method to `calculate_cs_self_peak_midplane_axial_stress` and update references in tests

* 🔄 - Rename variable `sig_axial` to `stress_z_cs_self_peak_midplane` and update references in PFCoil and CSCoil classes

* 🔄 - Rename axial_force variable to `forc_z_cs_self_peak_midplane` and update references in PFCoil and CSCoil classes

* 📝 Update axial stress calculation in CSCoil class to use average axial stress at coil interface

* Update tests

* 🔄 - Update plot_cs_coil_structure to include peak magnetic field, axial force, and stress in the plot output
🔄 - Refactor test cases for calculate_cs_self_peak_midplane_axial_stress to ensure proper assertions and formatting
diff --git a/documentation/eng-models/central-solenoid.md b/documentation/eng-models/central-solenoid.md
@@ -123,7 +123,38 @@ The peak field at the bore of the central solenoid will not be the same as that
 
 -----------
 
-### Axial stresses | `axial_stress()`
+### 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]:
+
+$$
+F_{z}(0)=\frac{\mu_0}{2}\left(\frac{N I}{2 \times dz_{\text{half}}}\right) \times \\
+ \left(2dz_{\text{half}} \sqrt{4r_{\text{CS,outer}}^2 + dz_{\text{half}}^2} \left[K(k_b) - E(k_b)\right]\\
+  - 2dz_{\text{half}} \sqrt{4r_{\text{CS,outer}}^2 + 4dz_{\text{half}}^2} \left[K(k_{2b}) - E(k_{2b})\right]  \right)
+$$
+
+where the modulus $k_{2b}$ is given by
+
+$$
+k_{2b} = \sqrt{\frac{4r_{\text{CS,outer}}^2}{4r_{\text{CS,outer}}^2+4dz_{\text{half}}^2}}
+$$
+
+and $k_b$ is given by:
+
+$$
+k_{2b} = \sqrt{\frac{4r_{\text{CS,outer}}^2}{4r_{\text{CS,outer}}^2+dz_{\text{half}}^2}}
+$$
+
+Here $K(k)$ and $E(k)$ are the complete elliptic integrals, respectively of the first and second kinds.
+
+
+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)$.
+
+
+
+--------------------------
 
 
 ### Hoop stress | `hoop_stress()`
@@ -290,4 +321,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.
 ‌
diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py
@@ -29,11 +29,13 @@
 
 ssq0: float = None
 
-sig_axial: float = None
+stress_z_cs_self_peak_midplane: float = None
+"""Peak axial stress (z) in central solenoid at midplane due to its own field (when at peak current) (Pa)"""
 
 sig_hoop: float = None
 
-axial_force: float = None
+forc_z_cs_self_peak_midplane: float = None
+"""Axial force (z) on central solenoid at midplane due to its own field (when at peak current) (N)"""
 
 r_pf_cs_current_filaments: list[float] = None
 """array of radial positions of current filaments in central solenoid"""
@@ -567,9 +569,9 @@ def init_pfcoil_module():
     global nfxf
     global ricpf
     global ssq0
-    global sig_axial
+    global stress_z_cs_self_peak_midplane
     global sig_hoop
-    global axial_force
+    global forc_z_cs_self_peak_midplane
     global r_pf_cs_current_filaments
     global z_pf_cs_current_filaments
     global c_pf_cs_current_filaments
@@ -587,9 +589,9 @@ def init_pfcoil_module():
     nfxf = 0
     ricpf = 0.0
     ssq0 = 0.0
-    sig_axial = 0.0
+    stress_z_cs_self_peak_midplane = 0.0
     sig_hoop = 0.0
-    axial_force = 0.0
+    forc_z_cs_self_peak_midplane = 0.0
 
     r_pf_cs_current_filaments = np.zeros(NFIXMX)
     z_pf_cs_current_filaments = np.zeros(NFIXMX)
diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py
@@ -9131,7 +9131,9 @@ def plot_cs_coil_structure(axis, fig, mfile_data, scan, colour_scheme=1):
         f"CS poloidal area: {mfile_data.data['a_cs_poloidal'].get_scan(scan):.4f} m$^2$\n"
         f"$N_{{\\text{{turns}}}}:$ {mfile_data.data['n_pf_coil_turns[n_cs_pf_coils-1]'].get_scan(scan):,.2f}\n"
         f"$I_{{\\text{{peak}}}}:$ {mfile_data.data['c_pf_cs_coils_peak_ma[n_cs_pf_coils-1]'].get_scan(scan):.3f}$ \\ MA$\n"
-        f"$B_{{\\text{{peak}}}}:$ {mfile_data.data['b_pf_coil_peak[n_cs_pf_coils-1]'].get_scan(scan):.3f}$ \\ T$\n\n"
+        f"$B_{{\\text{{peak}}}}:$ {mfile_data.data['b_pf_coil_peak[n_cs_pf_coils-1]'].get_scan(scan):.3f}$ \\ T$\n"
+        f"$F_{{\\text{{z,self,peak}}}}:$ {mfile_data.data['forc_z_cs_self_peak_midplane'].get_scan(scan) / 1e6:.3f}$ \\ MN$\n"
+        f"$\\sigma_{{\\text{{z,self,peak}}}}:$ {mfile_data.data['stress_z_cs_self_peak_midplane'].get_scan(scan) / 1e6:.3f}$ \\ MPa$ "
     )
 
     axis.text(
diff --git a/process/pfcoil.py b/process/pfcoil.py
@@ -2208,8 +2208,8 @@ def outpf(self):
                 op.ovarre(
                     self.outfile,
                     "Axial stress in CS steel (Pa)",
-                    "(sig_axial)",
-                    pfcoil_variables.sig_axial,
+                    "(stress_z_cs_self_peak_midplane)",
+                    pfcoil_variables.stress_z_cs_self_peak_midplane,
                     "OP ",
                 )
                 op.ovarre(
@@ -2222,8 +2222,8 @@ def outpf(self):
                 op.ovarre(
                     self.outfile,
                     "Axial force in CS (N)",
-                    "(axial_force)",
-                    pfcoil_variables.axial_force,
+                    "(forc_z_cs_self_peak_midplane)",
+                    pfcoil_variables.forc_z_cs_self_peak_midplane,
                     "OP ",
                 )
                 op.ovarre(
@@ -3343,8 +3343,21 @@ def ohcalc(self):
             )
 
             # New calculation from Y. Iwasa for axial stress
-            pfcoil_variables.sig_axial, pfcoil_variables.axial_force = (
-                self.axial_stress()
+            (
+                pfcoil_variables.stress_z_cs_self_peak_midplane,
+                pfcoil_variables.forc_z_cs_self_peak_midplane,
+            ) = self.calculate_cs_self_peak_midplane_axial_stress(
+                r_cs_outer=pfcoil_variables.r_pf_coil_outer[
+                    pfcoil_variables.n_cs_pf_coils - 1
+                ],
+                r_cs_inner=pfcoil_variables.r_pf_coil_inner[
+                    pfcoil_variables.n_cs_pf_coils - 1
+                ],
+                dz_cs_half=pfcoil_variables.dz_cs_full / 2.0,
+                c_cs_peak=pfcoil_variables.c_pf_cs_coils_peak_ma[
+                    pfcoil_variables.n_cs_pf_coils - 1
+                ]
+                * 1.0e6,
             )
 
             # Allowable (hoop) stress (Pa) alstroh
@@ -3372,8 +3385,11 @@ def ohcalc(self):
 
             if pfcoil_variables.i_cs_stress == 1:
                 pfcoil_variables.s_shear_cs_peak = max(
-                    abs(pfcoil_variables.sig_hoop - pfcoil_variables.sig_axial),
-                    abs(pfcoil_variables.sig_axial - 0.0e0),
+                    abs(
+                        pfcoil_variables.sig_hoop
+                        - pfcoil_variables.stress_z_cs_self_peak_midplane
+                    ),
+                    abs(pfcoil_variables.stress_z_cs_self_peak_midplane - 0.0e0),
                     abs(0.0e0 - pfcoil_variables.sig_hoop),
                 )
             else:
@@ -3679,65 +3695,81 @@ def output_cs_structure(self):
             "OP ",
         )
 
-    def axial_stress(self):
-        """Calculation of axial stress of central solenoid.
+    def calculate_cs_self_peak_midplane_axial_stress(
+        self,
+        r_cs_outer: float,
+        r_cs_inner: float,
+        dz_cs_half: float,
+        c_cs_peak: float,
+    ) -> tuple[float, float]:
+        """
+        Calculate axial stress and axial force for the central solenoid.
+
+        :param float r_cs_outer: Outer radius of the central solenoid (m).
+        :param float r_cs_inner: Inner radius of the central solenoid (m).
+        :param float dz_cs_half: Half-height of the central solenoid (m).
+        :param float c_cs_peak: Peak CS coil current (A).
 
-        author: J Morris, CCFE, Culham Science Centre
-        This routine calculates the axial stress of the central solenoid
-        from "Case studies in superconducting magnets", Y. Iwasa, Springer
+        :returns: A tuple containing the unsmeared axial stress and the axial force.
+        :rtype: tuple(float, float)
+            The first element is the unsmeared axial stress in MPa.
+            The second element is the axial force in newtons (N).
 
-        :return: unsmeared axial stress [MPa], axial force [N]
-        :rtype: tuple[float, float]
-        """
-        b = pfcoil_variables.r_pf_coil_outer[pfcoil_variables.n_cs_pf_coils - 1]
+        :note:
+           The axial force is computed using elliptic-integral based terms and the
+           unsmeared axial stress is obtained by dividing the axial force by
+           the effective steel area associated with the CS turns.
 
-        # Half height of central Solenoid [m]
-        hl = pfcoil_variables.z_pf_coil_upper[pfcoil_variables.n_cs_pf_coils - 1]
+        :references:
+            - Case Studies in Superconducting Magnets. Boston, MA: Springer US, 2009.
+              doi: https://doi.org/10.1007/b112047.
 
-        # Central Solenoid current [A]
-        ni = (
-            pfcoil_variables.c_pf_cs_coils_peak_ma[pfcoil_variables.n_cs_pf_coils - 1]
-            * 1.0e6
-        )
+        """
 
         # kb term for elliptical integrals
         # kb2 = SQRT((4.0e0*b**2)/(4.0e0*b**2 + hl**2))
-        kb2 = (4.0e0 * b**2) / (4.0e0 * b**2 + hl**2)
+        kb2 = (4.0e0 * r_cs_outer**2) / (4.0e0 * r_cs_outer**2 + dz_cs_half**2)
 
         # k2b term for elliptical integrals
         # k2b2 = SQRT((4.0e0*b**2)/(4.0e0*b**2 + 4.0e0*hl**2))
-        k2b2 = (4.0e0 * b**2) / (4.0e0 * b**2 + 4.0e0 * hl**2)
+        k2b2 = (4.0e0 * r_cs_outer**2) / (4.0e0 * r_cs_outer**2 + 4.0e0 * dz_cs_half**2)
 
         # term 1
-        axial_term_1 = -(constants.RMU0 / 2.0e0) * (ni / (2.0e0 * hl)) ** 2
+        axial_term_1 = (
+            -(constants.RMU0 / 2.0e0) * (c_cs_peak / (2.0e0 * dz_cs_half)) ** 2
+        )
 
         # term 2
         ekb2_1 = ellipk(kb2)
         ekb2_2 = ellipe(kb2)
         axial_term_2 = (
-            2.0e0 * hl * (math.sqrt(4.0e0 * b**2 + hl**2)) * (ekb2_1 - ekb2_2)
+            2.0e0
+            * dz_cs_half
+            * (math.sqrt(4.0e0 * r_cs_outer**2 + dz_cs_half**2))
+            * (ekb2_1 - ekb2_2)
         )
 
         # term 3
         ek2b2_1 = ellipk(k2b2)
         ek2b2_2 = ellipe(k2b2)
         axial_term_3 = (
-            2.0e0 * hl * (math.sqrt(4.0e0 * b**2 + 4.0e0 * hl**2)) * (ek2b2_1 - ek2b2_2)
+            2.0e0
+            * dz_cs_half
+            * (math.sqrt(4.0e0 * r_cs_outer**2 + 4.0e0 * dz_cs_half**2))
+            * (ek2b2_1 - ek2b2_2)
         )
 
         # calculate axial force [N]
-        axial_force = axial_term_1 * (axial_term_2 - axial_term_3)
+        forc_z_cs_self_peak_midplane = axial_term_1 * (axial_term_2 - axial_term_3)
 
         # axial area [m2]
-        area_ax = constants.PI * (
-            pfcoil_variables.r_pf_coil_outer[pfcoil_variables.n_cs_pf_coils - 1] ** 2
-            - pfcoil_variables.r_pf_coil_inner[pfcoil_variables.n_cs_pf_coils - 1] ** 2
-        )
+        area_ax = constants.PI * (r_cs_outer**2 - r_cs_inner**2)
 
-        # calculate unsmeared axial stress [MPa]
-        s_axial = axial_force / (pfcoil_variables.f_a_cs_turn_steel * 0.5 * area_ax)
+        # Calculate unsmeared axial stress
+        # Average axial stress at the interface of each half of the coil
+        s_axial = forc_z_cs_self_peak_midplane / (0.5 * area_ax)
 
-        return s_axial, axial_force
+        return s_axial, forc_z_cs_self_peak_midplane
 
     def hoop_stress(self, r):
         """Calculation of hoop stress of central solenoid.
diff --git a/tests/integration/test_pfcoil_int.py b/tests/integration/test_pfcoil_int.py
diff --git a/tests/unit/test_pfcoil.py b/tests/unit/test_pfcoil.py
