Add cs stress profiles

process/core/io/plot_proc.py

@@ -58,6 +58,7 @@
     vacuum_vessel_geometry_double_null,
     vacuum_vessel_geometry_single_null,
 )
+from process.models.pfcoil import CSCoil
 from process.models.physics.confinement_time import (
     ConfinementTimeModel,
     PlasmaConfinementTime,
@@ -9802,8 +9803,8 @@ def plot_cs_coil_structure(
     )
 
     axis.text(
-        0.6,
-        0.75,
+        0.45,
+        0.375,
         textstr_cs,
         fontsize=9,
         verticalalignment="bottom",
@@ -9917,8 +9918,8 @@ def plot_cs_turn_structure(axis: plt.Axes, fig, mfile: mf.MFile, scan: int):
     )
 
     axis.text(
-        0.6,
-        0.5,
+        0.7,
+        0.375,
         textstr_turn,
         fontsize=9,
         verticalalignment="bottom",
@@ -13878,11 +13879,15 @@ def main_plot(
 
     plot_current_profiles_over_time(figs[28].add_subplot(111), m_file, scan)
 
+    CSCoil.plot_stress_time_profile(
+        axis=figs[29].add_subplot(311), mfile=m_file, scan=scan
+    )
+
     plot_cs_coil_structure(
-        figs[29].add_subplot(121, aspect="equal"), figs[29], m_file, scan
+        figs[29].add_subplot(223, aspect="equal"), figs[29], m_file, scan
     )
     plot_cs_turn_structure(
-        figs[29].add_subplot(224, aspect="equal"), figs[29], m_file, scan
+        figs[29].add_subplot(326, aspect="equal"), figs[29], m_file, scan
     )
 
     plot_first_wall_top_down_cross_section(

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 = np.zeros(6)
     sig_hoop = 0.0
     forc_z_cs_self_peak_midplane = 0.0
 

process/models/pfcoil.py

@@ -1,12 +1,14 @@
 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
@@ -2541,6 +2543,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
@@ -3538,6 +3547,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,
@@ -3762,6 +3772,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.
 
@@ -3834,6 +3867,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