diff --git a/documentation/eng-models/central-solenoid.md b/documentation/eng-models/central-solenoid.md index 6fe9dcec63..44f5fe2a4d 100644 --- 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. ‌ \ No newline at end of file diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py index a0136c40ca..fa9d398da5 100644 --- 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 index 26a7291906..0b6f04e1f2 100644 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -9051,7 +9051,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 index 65dc4310bf..b50f421120 100644 --- 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 index c2bfc829c2..7180864642 100644 --- a/tests/integration/test_pfcoil_int.py +++ b/tests/integration/test_pfcoil_int.py @@ -2629,141 +2629,6 @@ def test_superconpf(monkeypatch: pytest.MonkeyPatch): assert pytest.approx(tmarg) == tmarg_exp -def test_axial_stress(cs_coil: CSCoil, monkeypatch: pytest.MonkeyPatch): - """Test axial_stress subroutine. - - axial_stress() requires specific mocked vars in order to work; these were - discovered using gdb to break on the first subroutine call when running the - baseline 2018 IN.DAT. - - :param pfcoil: a PFCoil instance - :type pfcoil: process.pfcoil.PFCoil - :param monkeypatch: mocking fixture - :type monkeypatch: _pytest.monkeypatch.MonkeyPatch - """ - monkeypatch.setattr(pfcoil_variables, "f_a_cs_turn_steel", 0.57874999999999999) - monkeypatch.setattr(pfcoil_variables, "n_cs_pf_coils", 7) - monkeypatch.setattr( - pfcoil_variables, - "r_pf_coil_outer", - np.array([ - 6.8520884119768697, - 6.9480065348448967, - 18.98258241468535, - 18.98258241468535, - 17.22166645654087, - 17.22166645654087, - 2.88462, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - ]), - ) - monkeypatch.setattr( - pfcoil_variables, - "r_pf_coil_inner", - np.array([ - 5.6944236847973242, - 5.5985055619292972, - 17.819978201682968, - 17.819978201682968, - 16.385123084628962, - 16.385123084628962, - 2.3321999999999998, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - ]), - ) - monkeypatch.setattr( - pfcoil_variables, - "c_pf_cs_coils_peak_ma", - np.array([ - 14.742063826112622, - 20.032681634901664, - -8.1098913365453491, - -8.1098913365453491, - -5.5984385047179153, - -5.5984385047179153, - -186.98751599968145, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - ]), - ) - monkeypatch.setattr( - pfcoil_variables, - "z_pf_coil_upper", - np.array([ - 10.184979073267192, - -11.815840508019832, - 3.4490763000495788, - -3.4490763000495788, - 8.4480394278914357, - -8.4480394278914357, - 8.1657810194058289, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - ]), - ) - - s_axial_exp = -7.468967e8 - axial_force_exp = -1.956801e9 - s_axial, axial_force = cs_coil.axial_stress() - - assert pytest.approx(s_axial) == s_axial_exp - assert pytest.approx(axial_force) == axial_force_exp - - def test_induct(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch): """Test induct subroutine. diff --git a/tests/unit/test_pfcoil.py b/tests/unit/test_pfcoil.py index ef46b0eba3..b5c2583ed8 100644 --- a/tests/unit/test_pfcoil.py +++ b/tests/unit/test_pfcoil.py @@ -2341,3 +2341,76 @@ def test_calculate_cs_turn_geometry_eu_demo_output_consistency( ) # The cable space plus conduit thickness should not exceed half the turn width assert radius_cs_turn_cable_space + dz_cs_turn_conduit <= dz_cs_turn / 2 + 1e-8 + + +def test_calculate_cs_self_peak_midplane_axial_stress_basic(cs_coil): + # Typical values for a CS coil + r_cs_outer = 2.0 + r_cs_inner = 1.5 + dz_cs_half = 1.0 + c_cs_peak = 1.2e7 # 12 MA + + s_axial, force_axial = cs_coil.calculate_cs_self_peak_midplane_axial_stress( + r_cs_outer, r_cs_inner, dz_cs_half, c_cs_peak + ) + + # Check actual output numbers + # These expected values should be updated if the implementation changes + s_axial_exp = -40124020.69015699 # Example expected value, update as needed + force_axial_exp = -110296662.5535968 # Example expected value, update as needed + + assert pytest.approx(s_axial, rel=1e-4) == s_axial_exp + assert pytest.approx(force_axial, rel=1e-4) == force_axial_exp + + +def test_calculate_cs_self_peak_midplane_axial_stress_zero_current(cs_coil): + # Zero current should yield zero force and stress + r_cs_outer = 2.0 + r_cs_inner = 1.5 + dz_cs_half = 1.0 + c_cs_peak = 0.0 + + s_axial, force_axial = cs_coil.calculate_cs_self_peak_midplane_axial_stress( + r_cs_outer, r_cs_inner, dz_cs_half, c_cs_peak + ) + + assert pytest.approx(s_axial) == 0.0 + assert pytest.approx(force_axial) == 0.0 + + +def test_calculate_cs_self_peak_midplane_axial_stress_extreme_geometry(cs_coil): + # Very thin coil, large dz_cs_half + r_cs_outer = 1.01 + r_cs_inner = 1.0 + dz_cs_half = 10.0 + c_cs_peak = 1.0e7 + + s_axial, force_axial = cs_coil.calculate_cs_self_peak_midplane_axial_stress( + r_cs_outer, r_cs_inner, dz_cs_half, c_cs_peak + ) + + # Check actual output numbers + s_axial_exp = -15804019.29465635 # Example expected value, update as needed + force_axial_exp = -498980.39867850166 # Example expected value, update as needed + + assert pytest.approx(s_axial, rel=1e-4) == s_axial_exp + assert pytest.approx(force_axial, rel=1e-4) == force_axial_exp + + +def test_calculate_cs_self_peak_midplane_axial_stress_invalid(cs_coil): + # Negative dz_cs_half should still return floats, but may be positive force + r_cs_outer = 2.0 + r_cs_inner = 1.5 + dz_cs_half = 2.0 + c_cs_peak = 1.2e7 + + s_axial, force_axial = cs_coil.calculate_cs_self_peak_midplane_axial_stress( + r_cs_outer, r_cs_inner, dz_cs_half, c_cs_peak + ) + + # Check actual output numbers + s_axial_exp = -16262350.341736654 # Example expected value, update as needed + force_axial_exp = -44703470.31824042 # Example expected value, update as needed + + assert pytest.approx(s_axial, rel=1e-4) == s_axial_exp + assert pytest.approx(force_axial, rel=1e-4) == force_axial_exp