From dac93d94c946e20e5ca59906ee2f0bbc818b6c6c Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 1 Apr 2026 10:43:52 +0100 Subject: [PATCH 1/2] Add axial stress profile variable for central solenoid midplane --- process/data_structure/pfcoil_variables.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/process/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py index 0223332a1c..fa262b2e1b 100644 --- a/process/data_structure/pfcoil_variables.py +++ b/process/data_structure/pfcoil_variables.py @@ -32,6 +32,9 @@ 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)""" +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 = None forc_z_cs_self_peak_midplane: float = None @@ -563,6 +566,7 @@ def init_pfcoil_module(): ricpf, \ ssq0, \ stress_z_cs_self_peak_midplane, \ + stress_z_cs_self_midplane_profile, \ sig_hoop, \ forc_z_cs_self_peak_midplane, \ r_pf_cs_current_filaments, \ @@ -583,6 +587,7 @@ def init_pfcoil_module(): ricpf = 0.0 ssq0 = 0.0 stress_z_cs_self_peak_midplane = 0.0 + stress_z_cs_self_midplane_profile = [] sig_hoop = 0.0 forc_z_cs_self_peak_midplane = 0.0 From f8b0363c3e8ac920d026bccd961b2d7b87fc9417 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 1 Apr 2026 11:26:54 +0100 Subject: [PATCH 2/2] Add central solenoid stress profile calculations and plotting --- process/core/io/plot/summary.py | 17 +++-- process/data_structure/pfcoil_variables.py | 2 +- process/models/pfcoil.py | 84 ++++++++++++++++++++++ 3 files changed, 96 insertions(+), 7 deletions(-) diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index b340d0ab20..c87437b6c1 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -47,6 +47,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, @@ -9763,8 +9764,8 @@ def plot_cs_coil_structure( ) axis.text( - 0.6, - 0.75, + 0.45, + 0.375, textstr_cs, fontsize=9, verticalalignment="bottom", @@ -9878,8 +9879,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", @@ -14448,11 +14449,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/data_structure/pfcoil_variables.py b/process/data_structure/pfcoil_variables.py index fa262b2e1b..3427b39083 100644 --- a/process/data_structure/pfcoil_variables.py +++ b/process/data_structure/pfcoil_variables.py @@ -587,7 +587,7 @@ def init_pfcoil_module(): ricpf = 0.0 ssq0 = 0.0 stress_z_cs_self_peak_midplane = 0.0 - stress_z_cs_self_midplane_profile = [] + stress_z_cs_self_midplane_profile = np.zeros(6) sig_hoop = 0.0 forc_z_cs_self_peak_midplane = 0.0 diff --git a/process/models/pfcoil.py b/process/models/pfcoil.py index 8448146963..e4bd9a2473 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 @@ -2538,6 +2541,13 @@ def outpf(self): f"(b_pf_coil_peak[{k}])", pfcoil_variables.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 @@ -3543,6 +3553,7 @@ def ohcalc(self): pfcoil_variables.p_pf_coil_resistive_total_flat_top += ( pfcoil_variables.p_cs_resistive_flat_top ) + self.calculate_cs_self_midplane_axial_stress_time_profile() def calculate_cs_self_peak_magnetic_field( self, @@ -3766,6 +3777,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. @@ -3838,6 +3872,56 @@ def hoop_stress(self, r): return s_hoop_nom / pfcoil_variables.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