From dda3ad9347d051a3dc6b94aa43ab4c088b01f1b6 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Wed, 29 Oct 2025 21:51:03 +0000 Subject: [PATCH 1/7] Refactor axial stress calculation in CSCoil class to improve parameter handling and documentation --- process/pfcoil.py | 79 +++++++++++++++++++++++++++++++---------------- 1 file changed, 53 insertions(+), 26 deletions(-) diff --git a/process/pfcoil.py b/process/pfcoil.py index 65dc4310bf..919215f8f1 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -3344,7 +3344,19 @@ def ohcalc(self): # New calculation from Y. Iwasa for axial stress pfcoil_variables.sig_axial, pfcoil_variables.axial_force = ( - self.axial_stress() + self.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 @@ -3679,60 +3691,75 @@ def output_cs_structure(self): "OP ", ) - def axial_stress(self): - """Calculation of axial stress of central solenoid. + def 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. - 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 + :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). - :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] + :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). - # Half height of central Solenoid [m] - hl = pfcoil_variables.z_pf_coil_upper[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. - # Central Solenoid current [A] - ni = ( - pfcoil_variables.c_pf_cs_coils_peak_ma[pfcoil_variables.n_cs_pf_coils - 1] - * 1.0e6 - ) + :references: + - Case Studies in Superconducting Magnets. Boston, MA: Springer US, 2009. + doi: https://doi.org/10.1007/b112047. + + """ # 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) # 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) From 70c801486cba08e67bb9f2787ab988811f5e84c7 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Wed, 29 Oct 2025 21:52:22 +0000 Subject: [PATCH 2/7] =?UTF-8?q?=F0=9F=94=84=20-=20Rename=20axial=5Fstress?= =?UTF-8?q?=20method=20to=20`calculate=5Fcs=5Fself=5Fpeak=5Fmidplane=5Faxi?= =?UTF-8?q?al=5Fstress`=20and=20update=20references=20in=20tests?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/pfcoil.py | 4 ++-- tests/integration/test_pfcoil_int.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/process/pfcoil.py b/process/pfcoil.py index 919215f8f1..4fa6a7a9f9 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -3344,7 +3344,7 @@ def ohcalc(self): # New calculation from Y. Iwasa for axial stress pfcoil_variables.sig_axial, pfcoil_variables.axial_force = ( - self.axial_stress( + self.calculate_cs_self_peak_midplane_axial_stress( r_cs_outer=pfcoil_variables.r_pf_coil_outer[ pfcoil_variables.n_cs_pf_coils - 1 ], @@ -3691,7 +3691,7 @@ def output_cs_structure(self): "OP ", ) - def axial_stress( + def calculate_cs_self_peak_midplane_axial_stress( self, r_cs_outer: float, r_cs_inner: float, diff --git a/tests/integration/test_pfcoil_int.py b/tests/integration/test_pfcoil_int.py index c2bfc829c2..a9d1b48426 100644 --- a/tests/integration/test_pfcoil_int.py +++ b/tests/integration/test_pfcoil_int.py @@ -2632,7 +2632,7 @@ def test_superconpf(monkeypatch: pytest.MonkeyPatch): 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 + calculate_cs_self_peak_midplane_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. @@ -2758,7 +2758,7 @@ def test_axial_stress(cs_coil: CSCoil, monkeypatch: pytest.MonkeyPatch): s_axial_exp = -7.468967e8 axial_force_exp = -1.956801e9 - s_axial, axial_force = cs_coil.axial_stress() + s_axial, axial_force = cs_coil.calculate_cs_self_peak_midplane_axial_stress() assert pytest.approx(s_axial) == s_axial_exp assert pytest.approx(axial_force) == axial_force_exp From 5524ce074f39e9ffe6a561480a082dd6427583e9 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Wed, 29 Oct 2025 21:55:40 +0000 Subject: [PATCH 3/7] =?UTF-8?q?=F0=9F=94=84=20-=20Rename=20variable=20`sig?= =?UTF-8?q?=5Faxial`=20to=20`stress=5Fz=5Fcs=5Fself=5Fpeak=5Fmidplane`=20a?= =?UTF-8?q?nd=20update=20references=20in=20PFCoil=20and=20CSCoil=20classes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/data_structure/pfcoil_variables.py | 7 ++-- process/pfcoil.py | 40 ++++++++++++---------- 2 files changed, 26 insertions(+), 21 deletions(-) diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py index a0136c40ca..c43a965c33 100644 --- a/process/data_structure/pfcoil_variables.py +++ b/process/data_structure/pfcoil_variables.py @@ -29,7 +29,8 @@ 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 @@ -567,7 +568,7 @@ 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 r_pf_cs_current_filaments @@ -587,7 +588,7 @@ 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 diff --git a/process/pfcoil.py b/process/pfcoil.py index 4fa6a7a9f9..384b79f0bc 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( @@ -3343,20 +3343,21 @@ def ohcalc(self): ) # New calculation from Y. Iwasa for axial stress - pfcoil_variables.sig_axial, pfcoil_variables.axial_force = ( - 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, - ) + ( + pfcoil_variables.stress_z_cs_self_peak_midplane, + pfcoil_variables.axial_force, + ) = 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 @@ -3384,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: From 1a356bc73fa8ba2a9c7b97655b671980a10b954c Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Wed, 29 Oct 2025 21:58:58 +0000 Subject: [PATCH 4/7] =?UTF-8?q?=F0=9F=94=84=20-=20Rename=20axial=5Fforce?= =?UTF-8?q?=20variable=20to=20`forc=5Fz=5Fcs=5Fself=5Fpeak=5Fmidplane`=20a?= =?UTF-8?q?nd=20update=20references=20in=20PFCoil=20and=20CSCoil=20classes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/data_structure/pfcoil_variables.py | 7 ++++--- process/pfcoil.py | 14 ++++++++------ tests/integration/test_pfcoil_int.py | 6 ++++-- 3 files changed, 16 insertions(+), 11 deletions(-) diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py index c43a965c33..fa9d398da5 100644 --- a/process/data_structure/pfcoil_variables.py +++ b/process/data_structure/pfcoil_variables.py @@ -34,7 +34,8 @@ 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""" @@ -570,7 +571,7 @@ def init_pfcoil_module(): global ssq0 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 @@ -590,7 +591,7 @@ def init_pfcoil_module(): ssq0 = 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/pfcoil.py b/process/pfcoil.py index 384b79f0bc..e6b36870cd 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -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( @@ -3345,7 +3345,7 @@ def ohcalc(self): # New calculation from Y. Iwasa for axial stress ( pfcoil_variables.stress_z_cs_self_peak_midplane, - pfcoil_variables.axial_force, + 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 @@ -3760,15 +3760,17 @@ def calculate_cs_self_peak_midplane_axial_stress( ) # 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 * (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) + s_axial = forc_z_cs_self_peak_midplane / ( + pfcoil_variables.f_a_cs_turn_steel * 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 a9d1b48426..09d1a816f3 100644 --- a/tests/integration/test_pfcoil_int.py +++ b/tests/integration/test_pfcoil_int.py @@ -2758,10 +2758,12 @@ def test_axial_stress(cs_coil: CSCoil, monkeypatch: pytest.MonkeyPatch): s_axial_exp = -7.468967e8 axial_force_exp = -1.956801e9 - s_axial, axial_force = cs_coil.calculate_cs_self_peak_midplane_axial_stress() + s_axial, forc_z_cs_self_peak_midplane = ( + cs_coil.calculate_cs_self_peak_midplane_axial_stress() + ) assert pytest.approx(s_axial) == s_axial_exp - assert pytest.approx(axial_force) == axial_force_exp + assert pytest.approx(forc_z_cs_self_peak_midplane) == axial_force_exp def test_induct(pfcoil: PFCoil, monkeypatch: pytest.MonkeyPatch): From 1188e971633054cfe391897011cdd9cdf80badb8 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Wed, 29 Oct 2025 22:30:28 +0000 Subject: [PATCH 5/7] :memo: Update axial stress calculation in CSCoil class to use average axial stress at coil interface --- documentation/eng-models/central-solenoid.md | 31 +++++++++++++++++++- process/pfcoil.py | 5 ++-- 2 files changed, 32 insertions(+), 4 deletions(-) diff --git a/documentation/eng-models/central-solenoid.md b/documentation/eng-models/central-solenoid.md index 6fe9dcec63..6aed63f436 100644 --- a/documentation/eng-models/central-solenoid.md +++ b/documentation/eng-models/central-solenoid.md @@ -123,7 +123,35 @@ 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}} +$$ + +The axial compressive force at $z$ in an isolated solenoid increases from 0 at $z = dz_{\text{half}}^2}$ +to the maximum at the midplane, $F_{z}(0)$. + + + +-------------------------- ### Hoop stress | `hoop_stress()` @@ -290,4 +318,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/pfcoil.py b/process/pfcoil.py index e6b36870cd..6d16ad9e3c 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -3766,9 +3766,8 @@ def calculate_cs_self_peak_midplane_axial_stress( area_ax = constants.PI * (r_cs_outer**2 - r_cs_inner**2) # calculate unsmeared axial stress [MPa] - s_axial = forc_z_cs_self_peak_midplane / ( - pfcoil_variables.f_a_cs_turn_steel * 0.5 * area_ax - ) + # 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, forc_z_cs_self_peak_midplane From 4930e2adc1355c4be51a03aa65165bf4a5e67728 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 30 Oct 2025 08:44:39 +0000 Subject: [PATCH 6/7] Update tests --- documentation/eng-models/central-solenoid.md | 5 +- process/pfcoil.py | 2 +- tests/integration/test_pfcoil_int.py | 136 ------------------- tests/unit/test_pfcoil.py | 69 ++++++++++ 4 files changed, 74 insertions(+), 138 deletions(-) diff --git a/documentation/eng-models/central-solenoid.md b/documentation/eng-models/central-solenoid.md index 6aed63f436..44f5fe2a4d 100644 --- a/documentation/eng-models/central-solenoid.md +++ b/documentation/eng-models/central-solenoid.md @@ -146,7 +146,10 @@ $$ k_{2b} = \sqrt{\frac{4r_{\text{CS,outer}}^2}{4r_{\text{CS,outer}}^2+dz_{\text{half}}^2}} $$ -The axial compressive force at $z$ in an isolated solenoid increases from 0 at $z = 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)$. diff --git a/process/pfcoil.py b/process/pfcoil.py index 6d16ad9e3c..b50f421120 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -3765,7 +3765,7 @@ def calculate_cs_self_peak_midplane_axial_stress( # axial area [m2] area_ax = constants.PI * (r_cs_outer**2 - r_cs_inner**2) - # calculate unsmeared axial stress [MPa] + # 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) diff --git a/tests/integration/test_pfcoil_int.py b/tests/integration/test_pfcoil_int.py index 09d1a816f3..0ecb1306c8 100644 --- a/tests/integration/test_pfcoil_int.py +++ b/tests/integration/test_pfcoil_int.py @@ -2629,142 +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. - - calculate_cs_self_peak_midplane_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, forc_z_cs_self_peak_midplane = ( - cs_coil.calculate_cs_self_peak_midplane_axial_stress() - ) - - assert pytest.approx(s_axial) == s_axial_exp - assert pytest.approx(forc_z_cs_self_peak_midplane) == 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..92f00562a0 100644 --- a/tests/unit/test_pfcoil.py +++ b/tests/unit/test_pfcoil.py @@ -2341,3 +2341,72 @@ 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 \ No newline at end of file From d69ca429d3869ff65c8f55c3279a8032cabbbabf Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 30 Oct 2025 08:56:02 +0000 Subject: [PATCH 7/7] =?UTF-8?q?=F0=9F=94=84=20-=20Update=20plot=5Fcs=5Fcoi?= =?UTF-8?q?l=5Fstructure=20to=20include=20peak=20magnetic=20field,=20axial?= =?UTF-8?q?=20force,=20and=20stress=20in=20the=20plot=20output=20?= =?UTF-8?q?=F0=9F=94=84=20-=20Refactor=20test=20cases=20for=20calculate=5F?= =?UTF-8?q?cs=5Fself=5Fpeak=5Fmidplane=5Faxial=5Fstress=20to=20ensure=20pr?= =?UTF-8?q?oper=20assertions=20and=20formatting?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/io/plot_proc.py | 4 +++- tests/integration/test_pfcoil_int.py | 1 - tests/unit/test_pfcoil.py | 8 ++++++-- 3 files changed, 9 insertions(+), 4 deletions(-) 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/tests/integration/test_pfcoil_int.py b/tests/integration/test_pfcoil_int.py index 0ecb1306c8..7180864642 100644 --- a/tests/integration/test_pfcoil_int.py +++ b/tests/integration/test_pfcoil_int.py @@ -2629,7 +2629,6 @@ def test_superconpf(monkeypatch: pytest.MonkeyPatch): assert pytest.approx(tmarg) == tmarg_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 92f00562a0..b5c2583ed8 100644 --- a/tests/unit/test_pfcoil.py +++ b/tests/unit/test_pfcoil.py @@ -2342,6 +2342,7 @@ 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 @@ -2361,6 +2362,7 @@ def test_calculate_cs_self_peak_midplane_axial_stress_basic(cs_coil): 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 @@ -2375,6 +2377,7 @@ def test_calculate_cs_self_peak_midplane_axial_stress_zero_current(cs_coil): 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 @@ -2393,6 +2396,7 @@ def test_calculate_cs_self_peak_midplane_axial_stress_extreme_geometry(cs_coil): 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 @@ -2406,7 +2410,7 @@ def test_calculate_cs_self_peak_midplane_axial_stress_invalid(cs_coil): # 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 + 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 \ No newline at end of file + assert pytest.approx(force_axial, rel=1e-4) == force_axial_exp