From 5567a32127cfe357a677339ecddc1897ad95d139 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 1 Apr 2026 10:43:52 +0100 Subject: [PATCH 1/4] Add axial stress profile variable for central solenoid midplane --- process/data_structure/pfcoil_variables.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py index a20efcfa4..593c62703 100644 --- a/process/data_structure/pfcoil_variables.py +++ b/process/data_structure/pfcoil_variables.py @@ -37,6 +37,9 @@ class PFCoilData: stress_z_cs_self_peak_midplane: float = 0.0 """Peak axial stress (z) in central solenoid at midplane due to its own field (when at peak current) (Pa)""" +stress_z_cs_self_midplane_profile: list[float] = None +"""Axial stress (z) in central solenoid at midplane due to its own field at each time point (Pa)""" + sig_hoop: float = 0.0 forc_z_cs_self_peak_midplane: float = 0.0 From 90f2e3a7a1d767e6fe2e3164e7d1d7d47d1192e3 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 1 Apr 2026 11:26:54 +0100 Subject: [PATCH 2/4] Add central solenoid stress profile calculations and plotting --- process/core/io/plot/summary.py | 17 ++++--- process/models/pfcoil.py | 84 +++++++++++++++++++++++++++++++++ 2 files changed, 95 insertions(+), 6 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index bce8c6173..5a4a76194 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -48,6 +48,7 @@ vacuum_vessel_geometry_single_null, ) from process.models.physics.bootstrap_current import BootstrapCurrentFractionModel +from process.models.pfcoil import CSCoil from process.models.physics.confinement_time import ( ConfinementTimeModel, PlasmaConfinementTime, @@ -9762,8 +9763,8 @@ def plot_cs_coil_structure( ) axis.text( - 0.6, - 0.75, + 0.45, + 0.375, textstr_cs, fontsize=9, verticalalignment="bottom", @@ -9877,8 +9878,8 @@ def plot_cs_turn_structure(axis: plt.Axes, fig, mfile: MFile, scan: int): ) axis.text( - 0.6, - 0.5, + 0.7, + 0.375, textstr_turn, fontsize=9, verticalalignment="bottom", @@ -14447,11 +14448,15 @@ def main_plot( plot_current_profiles_over_time(figs[29].add_subplot(111), m_file, scan) + CSCoil.plot_stress_time_profile( + axis=figs[30].add_subplot(311), mfile=m_file, scan=scan + ) + plot_cs_coil_structure( - figs[30].add_subplot(121, aspect="equal"), figs[30], m_file, scan + figs[30].add_subplot(223, aspect="equal"), figs[30], m_file, scan ) plot_cs_turn_structure( - figs[30].add_subplot(224, aspect="equal"), figs[30], m_file, scan + figs[30].add_subplot(326, aspect="equal"), figs[30], m_file, scan ) plot_first_wall_top_down_cross_section( diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py index 3d7731a2b..cd1171b62 100644 --- a/process/models/pfcoil.py +++ b/process/models/pfcoil.py @@ -1,12 +1,15 @@ import logging import math +import matplotlib.pyplot as plt import numba import numpy as np from scipy import optimize from scipy.linalg import svd from scipy.special import ellipe, ellipk +import process.core.io.mfile as mf +import process.models.superconductors as superconductors from process.core import constants from process.core import process_output as op from process.core.exceptions import ProcessValueError @@ -2543,6 +2546,13 @@ def outpf(self): f"(b_pf_coil_peak[{k}])", self.data.pf_coil.b_pf_coil_peak[k], ) + for time in range(6): + op.ovarre( + self.mfile, + f"CS coil midplane axial stress at time point {time} (MPa)", + f"(stress_z_cs_self_midplane_profile[{time}])", + pfcoil_variables.stress_z_cs_self_midplane_profile[time], + ) self.tf_pf_collision_detector() # Central Solenoid, if present @@ -3558,6 +3568,7 @@ def ohcalc(self): self.data.pf_coil.p_pf_coil_resistive_total_flat_top += ( self.data.pf_coil.p_cs_resistive_flat_top ) + self.calculate_cs_self_midplane_axial_stress_time_profile() def calculate_cs_self_peak_magnetic_field( self, @@ -3781,6 +3792,29 @@ def calculate_cs_self_peak_midplane_axial_stress( return s_axial, forc_z_cs_self_peak_midplane + def calculate_cs_self_midplane_axial_stress_time_profile( + self, + ) -> None: + for time in range(6): + stress_value, _ = 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_coil_turn[ + pfcoil_variables.n_cs_pf_coils - 1, time + ] + * pfcoil_variables.n_pf_coil_turns[ + pfcoil_variables.n_cs_pf_coils - 1 + ] + ), + ) + pfcoil_variables.stress_z_cs_self_midplane_profile[time] = stress_value + def hoop_stress(self, r): """Calculation of hoop stress of central solenoid. @@ -3853,6 +3887,56 @@ def hoop_stress(self, r): return s_hoop_nom / self.data.pf_coil.f_a_cs_turn_steel + @staticmethod + def plot_stress_time_profile(axis: plt.Axes, mfile: mf.MFile, scan: int): + t_plant_pulse_coil_precharge = mfile.get( + "t_plant_pulse_coil_precharge", scan=scan + ) + t_plant_pulse_plasma_current_ramp_up = mfile.get( + "t_plant_pulse_plasma_current_ramp_up", scan=scan + ) + t_plant_pulse_fusion_ramp = mfile.get("t_plant_pulse_fusion_ramp", scan=scan) + t_plant_pulse_burn = mfile.get("t_plant_pulse_burn", scan=scan) + t_plant_pulse_plasma_current_ramp_down = mfile.get( + "t_plant_pulse_plasma_current_ramp_down", scan=scan + ) + + # Define a cumulative sum list for each point in the pulse + t_steps = np.cumsum([ + 0, + t_plant_pulse_coil_precharge, + t_plant_pulse_plasma_current_ramp_up, + t_plant_pulse_fusion_ramp, + t_plant_pulse_burn, + t_plant_pulse_plasma_current_ramp_down, + ]) + + stress_times = t_steps[ + :6 + ] # Get the first 6 time points corresponding to the stress profile + + stress_z_cs_self_midplane_profile = np.zeros(6) + for i in range(6): + stress_z_cs_self_midplane_profile[i] = mfile.get( + f"stress_z_cs_self_midplane_profile[{i}]", scan=scan + ) + + # Plot stress vs time + axis.plot( + stress_times, + stress_z_cs_self_midplane_profile / 1e6, + "o-", + linewidth=2, + markersize=8, + label="Midplane Axial Stress", + ) + axis.set_xlabel("Pulse Time (s)") + axis.set_ylabel("Stress (MPa)") + axis.minorticks_on() + axis.legend(loc="best") + axis.set_title("Central Solenoid Stress") + axis.grid(True, alpha=0.3) + def peak_b_field_at_pf_coil( n_coil: int, n_coil_group: int, t_b_field_peak: int, data: DataStructure From b7c9e4ace21a6cfaed3432da977ac9cbebc7a254 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 15 May 2026 13:17:02 +0100 Subject: [PATCH 3/4] Refactor import statements in summary.py and pfcoil.py for clarity --- process/core/io/plot/summary.py | 2 +- process/models/pfcoil.py | 1 - 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 5a4a76194..4c44ee325 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -47,8 +47,8 @@ vacuum_vessel_geometry_double_null, vacuum_vessel_geometry_single_null, ) -from process.models.physics.bootstrap_current import BootstrapCurrentFractionModel from process.models.pfcoil import CSCoil +from process.models.physics.bootstrap_current import BootstrapCurrentFractionModel from process.models.physics.confinement_time import ( ConfinementTimeModel, PlasmaConfinementTime, diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py index cd1171b62..de9924c1a 100644 --- a/process/models/pfcoil.py +++ b/process/models/pfcoil.py @@ -9,7 +9,6 @@ from scipy.special import ellipe, ellipk import process.core.io.mfile as mf -import process.models.superconductors as superconductors from process.core import constants from process.core import process_output as op from process.core.exceptions import ProcessValueError From bf3f817a77643268d6b49668c628076ceaf0790d Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 18 May 2026 13:32:16 +0100 Subject: [PATCH 4/4] Post rebase changes --- process/core/io/plot/summary.py | 54 +++++++++++++- process/data_structure/pfcoil_variables.py | 6 +- process/models/pfcoil.py | 87 +++++----------------- 3 files changed, 72 insertions(+), 75 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index 4c44ee325..ac1e0951d 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -47,7 +47,6 @@ 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, @@ -9798,6 +9797,55 @@ def plot_cs_coil_structure( axis.legend() +def plot_cs_stress_time_profile(axis: plt.Axes, mfile: MFile, scan: int) -> None: + """Function to plot the time profile of the CS stress during the pulse.""" + t_plant_pulse_coil_precharge = mfile.get("t_plant_pulse_coil_precharge", scan=scan) + t_plant_pulse_plasma_current_ramp_up = mfile.get( + "t_plant_pulse_plasma_current_ramp_up", scan=scan + ) + t_plant_pulse_fusion_ramp = mfile.get("t_plant_pulse_fusion_ramp", scan=scan) + t_plant_pulse_burn = mfile.get("t_plant_pulse_burn", scan=scan) + t_plant_pulse_plasma_current_ramp_down = mfile.get( + "t_plant_pulse_plasma_current_ramp_down", scan=scan + ) + + # Define a cumulative sum list for each point in the pulse + t_steps = np.cumsum([ + 0, + t_plant_pulse_coil_precharge, + t_plant_pulse_plasma_current_ramp_up, + t_plant_pulse_fusion_ramp, + t_plant_pulse_burn, + t_plant_pulse_plasma_current_ramp_down, + ]) + + stress_times = t_steps[ + :6 + ] # Get the first 6 time points corresponding to the stress profile + + stress_z_cs_self_midplane_profile = np.zeros(6) + for i in range(6): + stress_z_cs_self_midplane_profile[i] = mfile.get( + f"stress_z_cs_self_midplane_profile[{i}]", scan=scan + ) + + # Plot stress vs time + axis.plot( + stress_times, + stress_z_cs_self_midplane_profile / 1e6, + "o-", + linewidth=2, + markersize=8, + label="Midplane Axial Stress", + ) + axis.set_xlabel("Pulse Time (s)") + axis.set_ylabel("Stress (MPa)") + axis.minorticks_on() + axis.legend(loc="best") + axis.set_title("Central Solenoid Stress") + axis.grid(True, alpha=0.3) + + def plot_cs_turn_structure(axis: plt.Axes, fig, mfile: MFile, scan: int): a_cs_turn = mfile.get("a_cs_turn", scan=scan) dz_cs_turn = mfile.get("dz_cs_turn", scan=scan) @@ -14448,9 +14496,7 @@ def main_plot( plot_current_profiles_over_time(figs[29].add_subplot(111), m_file, scan) - CSCoil.plot_stress_time_profile( - axis=figs[30].add_subplot(311), mfile=m_file, scan=scan - ) + plot_cs_stress_time_profile(axis=figs[30].add_subplot(311), mfile=m_file, scan=scan) plot_cs_coil_structure( figs[30].add_subplot(223, aspect="equal"), figs[30], m_file, scan diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py index 593c62703..720e813ae 100644 --- a/process/data_structure/pfcoil_variables.py +++ b/process/data_structure/pfcoil_variables.py @@ -37,8 +37,10 @@ class PFCoilData: stress_z_cs_self_peak_midplane: float = 0.0 """Peak axial stress (z) in central solenoid at midplane due to its own field (when at peak current) (Pa)""" -stress_z_cs_self_midplane_profile: list[float] = None -"""Axial stress (z) in central solenoid at midplane due to its own field at each time point (Pa)""" + stress_z_cs_self_midplane_profile: list[float] = field( + default_factory=lambda: np.zeros(6) + ) + """Axial stress (z) in central solenoid at midplane due to its own field at each time point (Pa)""" sig_hoop: float = 0.0 diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py index de9924c1a..71f73afe9 100644 --- a/process/models/pfcoil.py +++ b/process/models/pfcoil.py @@ -1,14 +1,12 @@ import logging import math -import matplotlib.pyplot as plt import numba import numpy as np from scipy import optimize from scipy.linalg import svd from scipy.special import ellipe, ellipk -import process.core.io.mfile as mf from process.core import constants from process.core import process_output as op from process.core.exceptions import ProcessValueError @@ -2550,7 +2548,7 @@ def outpf(self): self.mfile, f"CS coil midplane axial stress at time point {time} (MPa)", f"(stress_z_cs_self_midplane_profile[{time}])", - pfcoil_variables.stress_z_cs_self_midplane_profile[time], + self.data.pf_coil.stress_z_cs_self_midplane_profile[time], ) self.tf_pf_collision_detector() @@ -3736,15 +3734,16 @@ def calculate_cs_self_peak_midplane_axial_stress( The first element is the unsmeared axial stress in MPa. The second element is the axial force in newtons (N). - :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. + Notes + ----- + 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. References ---------- - - Case Studies in Superconducting Magnets. Boston, MA: Springer US, 2009. - doi: https://doi.org/10.1007/b112047. + [1] 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)) @@ -3796,23 +3795,23 @@ def calculate_cs_self_midplane_axial_stress_time_profile( ) -> None: for time in range(6): stress_value, _ = 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_outer=self.data.pf_coil.r_pf_coil_outer[ + self.data.pf_coil.n_cs_pf_coils - 1 ], - r_cs_inner=pfcoil_variables.r_pf_coil_inner[ - pfcoil_variables.n_cs_pf_coils - 1 + r_cs_inner=self.data.pf_coil.r_pf_coil_inner[ + self.data.pf_coil.n_cs_pf_coils - 1 ], - dz_cs_half=pfcoil_variables.dz_cs_full / 2.0, + dz_cs_half=self.data.pf_coil.dz_cs_full / 2.0, c_cs_peak=( - pfcoil_variables.c_pf_coil_turn[ - pfcoil_variables.n_cs_pf_coils - 1, time + self.data.pf_coil.c_pf_coil_turn[ + self.data.pf_coil.n_cs_pf_coils - 1, time ] - * pfcoil_variables.n_pf_coil_turns[ - pfcoil_variables.n_cs_pf_coils - 1 + * self.data.pf_coil.n_pf_coil_turns[ + self.data.pf_coil.n_cs_pf_coils - 1 ] ), ) - pfcoil_variables.stress_z_cs_self_midplane_profile[time] = stress_value + self.data.pf_coil.stress_z_cs_self_midplane_profile[time] = stress_value def hoop_stress(self, r): """Calculation of hoop stress of central solenoid. @@ -3886,56 +3885,6 @@ def hoop_stress(self, r): return s_hoop_nom / self.data.pf_coil.f_a_cs_turn_steel - @staticmethod - def plot_stress_time_profile(axis: plt.Axes, mfile: mf.MFile, scan: int): - t_plant_pulse_coil_precharge = mfile.get( - "t_plant_pulse_coil_precharge", scan=scan - ) - t_plant_pulse_plasma_current_ramp_up = mfile.get( - "t_plant_pulse_plasma_current_ramp_up", scan=scan - ) - t_plant_pulse_fusion_ramp = mfile.get("t_plant_pulse_fusion_ramp", scan=scan) - t_plant_pulse_burn = mfile.get("t_plant_pulse_burn", scan=scan) - t_plant_pulse_plasma_current_ramp_down = mfile.get( - "t_plant_pulse_plasma_current_ramp_down", scan=scan - ) - - # Define a cumulative sum list for each point in the pulse - t_steps = np.cumsum([ - 0, - t_plant_pulse_coil_precharge, - t_plant_pulse_plasma_current_ramp_up, - t_plant_pulse_fusion_ramp, - t_plant_pulse_burn, - t_plant_pulse_plasma_current_ramp_down, - ]) - - stress_times = t_steps[ - :6 - ] # Get the first 6 time points corresponding to the stress profile - - stress_z_cs_self_midplane_profile = np.zeros(6) - for i in range(6): - stress_z_cs_self_midplane_profile[i] = mfile.get( - f"stress_z_cs_self_midplane_profile[{i}]", scan=scan - ) - - # Plot stress vs time - axis.plot( - stress_times, - stress_z_cs_self_midplane_profile / 1e6, - "o-", - linewidth=2, - markersize=8, - label="Midplane Axial Stress", - ) - axis.set_xlabel("Pulse Time (s)") - axis.set_ylabel("Stress (MPa)") - axis.minorticks_on() - axis.legend(loc="best") - axis.set_title("Central Solenoid Stress") - axis.grid(True, alpha=0.3) - def peak_b_field_at_pf_coil( n_coil: int, n_coil_group: int, t_b_field_peak: int, data: DataStructure