diff --git a/process/core/io/plot/summary.py b/process/core/io/plot/summary.py index bce8c6173..f0afe4ae2 100644 --- a/process/core/io/plot/summary.py +++ b/process/core/io/plot/summary.py @@ -75,6 +75,7 @@ ) from process.models.superconductors import SuperconductorModel from process.models.tfcoil.base import TFCoilShapeModel, TFPlasmaCaseType +from process.models.tfcoil.quench import plot_quench_time_evolution @dataclass @@ -7474,7 +7475,8 @@ def _pack_strands_rectangular_with_obstacles( textstr_superconductor = ( f"$\\mathbf{{Superconductor:}}$\n \n" - f"Superconductor used: {SuperconductorModel(mfile.get('i_tf_sc_mat', scan=scan)).full_name}\n" + f"Superconductor used: \n" + f"{SuperconductorModel(mfile.get('i_tf_sc_mat', scan=scan)).full_name}\n" f"Critical field at zero \ntemperature and strain: {mfile.get('b_tf_superconductor_critical_zero_temp_strain', scan=scan):.4f} T\n" f"Critical temperature at \nzero field and strain: {mfile.get('temp_tf_superconductor_critical_zero_field_strain', scan=scan):.4f} K\n" f"Temperature at conductor: {mfile.get('tftmp', scan=scan):.4f} K\n" @@ -7483,8 +7485,8 @@ def _pack_strands_rectangular_with_obstacles( f"Critcal current ratio: {mfile.get('f_c_tf_turn_operating_critical', scan=scan):,.4f}\n" f"Superconductor temperature \nmargin: {mfile.get('temp_tf_superconductor_margin', scan=scan):,.4f} K\n" f"\n$\\mathbf{{Quench:}}$\n \n" - f"Quench dump time: {mfile.get('t_tf_superconductor_quench', scan=scan):.4e} s\n" - f"Quench detection time: {mfile.get('t_tf_quench_detection', scan=scan):.4e} s\n" + f"Quench dump time: {mfile.get('t_tf_superconductor_quench', scan=scan):.4f} s\n" + f"Quench detection time: {mfile.get('t_tf_quench_detection', scan=scan):.4f} s\n" f"User input max temperature \nduring quench: {mfile.get('temp_tf_conductor_quench_max', scan=scan):.2f} K\n" f"Required maxium WP current \ndensity for heat protection:\n{mfile.get('j_tf_wp_quench_heat_max', scan=scan):.2e} A/m$^2$\n" ) @@ -14432,46 +14434,60 @@ def main_plot( plot_205 = figs[25].add_subplot(223, aspect="equal") plot_205.set_position([0.075, 0.1, 0.3, 0.3]) plot_cable_in_conduit_cable(plot_205, figs[25], m_file, scan) + plot_quench_time_evolution( + tau_discharge=m_file.get("t_tf_superconductor_quench", scan=scan), + b_peak=m_file.get("b_tf_inboard_peak_with_ripple", scan=scan), + f_cu=m_file.get("f_a_tf_turn_cable_copper", scan=scan), + f_he=m_file.get("f_a_tf_turn_cable_space_cooling", scan=scan), + t_he_peak=m_file.get("tftmp", scan=scan), + temp_quench_max=m_file.get("temp_tf_conductor_quench_max", scan=scan), + cu_rrr=m_file.get("rrr_tf_cu", scan=scan), + t_quench_detection=m_file.get("t_tf_quench_detection", scan=scan), + fluence=m_file.get("nflutfmax", scan=scan), + j_operating=m_file.get("j_tf_wp", scan=scan), + axes_1=figs[26].add_subplot(211), + axes_2=figs[26].add_subplot(212), + ) else: ax19 = figs[24].add_subplot(211, aspect="equal") ax19.set_position([0.06, 0.55, 0.675, 0.4]) plot_resistive_tf_wp(ax19, m_file, scan, figs[24]) plot_resistive_tf_info(ax19, m_file, scan, figs[24]) plot_tf_coil_structure( - figs[26].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme + figs[27].add_subplot(111, aspect="equal"), m_file, scan, colour_scheme ) - plot_plasma_outboard_toroidal_ripple_map(figs[27], m_file, scan) + plot_plasma_outboard_toroidal_ripple_map(figs[28], m_file, scan) - plot_tf_stress(figs[28].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) + plot_tf_stress(figs[29].subplots(nrows=3, ncols=1, sharex=True).flatten(), m_file) - plot_current_profiles_over_time(figs[29].add_subplot(111), m_file, scan) + plot_current_profiles_over_time(figs[30].add_subplot(111), m_file, scan) plot_cs_coil_structure( - figs[30].add_subplot(121, aspect="equal"), figs[30], m_file, scan + figs[31].add_subplot(121, aspect="equal"), figs[31], m_file, scan ) plot_cs_turn_structure( - figs[30].add_subplot(224, aspect="equal"), figs[30], m_file, scan + figs[31].add_subplot(224, aspect="equal"), figs[31], m_file, scan ) plot_first_wall_top_down_cross_section( - figs[31].add_subplot(221, aspect="equal"), m_file, scan + figs[32].add_subplot(221, aspect="equal"), m_file, scan ) - plot_first_wall_poloidal_cross_section(figs[31].add_subplot(122), m_file, scan) - plot_fw_90_deg_pipe_bend(figs[31].add_subplot(337), m_file, scan) + plot_first_wall_poloidal_cross_section(figs[32].add_subplot(122), m_file, scan) + plot_fw_90_deg_pipe_bend(figs[32].add_subplot(337), m_file, scan) - plot_blkt_pipe_bends(figs[32], m_file, scan) - ax_blanket = figs[32].add_subplot(122, aspect="equal") - plot_blkt_structure(ax_blanket, figs[32], m_file, scan, radial_build, colour_scheme) + plot_blkt_pipe_bends(figs[33], m_file, scan) + ax_blanket = figs[33].add_subplot(122, aspect="equal") + plot_blkt_structure(ax_blanket, figs[33], m_file, scan, radial_build, colour_scheme) plot_main_power_flow( - figs[33].add_subplot(111, aspect="equal"), m_file, scan, figs[33] + figs[34].add_subplot(111, aspect="equal"), m_file, scan, figs[34] ) - ax24 = figs[34].add_subplot(111) + ax24 = figs[35].add_subplot(111) # set_position([left, bottom, width, height]) -> height ~ 0.66 => ~2/3 of page height ax24.set_position([0.08, 0.35, 0.84, 0.57]) - plot_system_power_profiles_over_time(ax24, m_file, scan, figs[34]) + plot_system_power_profiles_over_time(ax24, m_file, scan, figs[35]) def create_thickness_builds(m_file, scan: int): @@ -14553,7 +14569,7 @@ def plot_summary( ): # create main plot # Increase range when adding new page - pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(35)] + pages = [plt.figure(figsize=(12, 9), dpi=80) for i in range(36)] # run main_plot main_plot( diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 79f1feb67..4ce0c79db 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -127,7 +127,7 @@ * (32) TF coil conduit stress upper limit (SCTF) (itv 49,56,57,58,59,60,24) * (33) I_op / I_critical (TF coil) (SCTF) (itv 50,56,57,58,59,60,24) * (34) Dump voltage upper limit (SCTF) (itv 51,52,56,57,58,59,60,24) -* (35) J_winding pack/J_protection upper limit (SCTF) (itv 53,56,57,58,59,60,24) +* (35) TF Quench Hotspot J limit (SCTF) (itv 53,56,57,58,59,60,24) * (36) TF coil temperature margin lower limit (SCTF) (itv 54,55,56,57,58,59,60,24) * (37) Current drive gamma upper limit (itv 40,47) * (38) First wall coolant temperature rise upper limit (itv 62) @@ -560,7 +560,7 @@ def init_numerics(): "TF coil conduit stress upper lim ", "I_op / I_critical (TF coil) ", "Dump voltage upper limit ", - "J_winding pack/J_protection limit", + "TF Quench Hotspot J limit ", "TF coil temp. margin lower limit ", "Current drive gamma limit ", "1st wall coolant temp rise limit ", diff --git a/process/models/tfcoil/quench.py b/process/models/tfcoil/quench.py index e1a31930f..6ce0e1ae5 100644 --- a/process/models/tfcoil/quench.py +++ b/process/models/tfcoil/quench.py @@ -3,6 +3,7 @@ from typing import Final from warnings import warn +import matplotlib.pyplot as plt import numpy as np # TODO: Use of CoolProp prevents nb.jit at present... @@ -12,8 +13,8 @@ # Material property parameterisations -COPPER_DENSITY = 8960.0 # [kg/m^3] -NB3SN_DENSITY = 8040.0 # [kg/m^3] +COPPER_DENSITY = 8960.0 # [kg/m³] +NB3SN_DENSITY = 8040.0 # [kg/m³] def _copper_specific_heat_capacity(temperature: float) -> float: @@ -36,11 +37,11 @@ def _copper_specific_heat_capacity(temperature: float) -> float: References ---------- - - J. Simon, E. S. Drexler, and R. P. Reed, *NIST Monograph 177*, "Properties of - Copper and Copper Alloys at Cryogenic Temperatures", - U.S. Government Printing Office, February 1992. - https://nvlpubs.nist.gov/nistpubs/Legacy/MONO/nistmonograph177.pdf - Equation 7-1 + [1] J. Simon, E. S. Drexler, and R. P. Reed, *NIST Monograph 177*, "Properties of + Copper and Copper Alloys at Cryogenic Temperatures", + U.S. Government Printing Office, February 1992. + https://nvlpubs.nist.gov/nistpubs/Legacy/MONO/nistmonograph177.pdf + Equation 7-1 """ poly_coeffs: Final[list[float]] = [1.131, -9.454, 12.99, -5.501, 0.7637] logt = np.log10(temperature) @@ -66,16 +67,14 @@ def _copper_rrr_resistivity(temperature: float, rrr: float) -> float: References ---------- - - J. Simon, E. S. Drexler, and R. P. Reed, *NIST Monograph 177*, "Properties of - Copper and Copper Alloys at Cryogenic Temperatures", - U.S. Government Printing Office, February 1992. - https://nvlpubs.nist.gov/nistpubs/Legacy/MONO/nistmonograph177.pdf - Equation 8-1 - - - J. G. Hust, A. B. Lankford, NBSIR 84-3007 "THERMAL CONDUCTIVITY OF - ALUMINUM, COPPER, IRON, AND TUNGSTEN FOR TEMPERATURES FROM 1 K TO THE - MELTING POINT", 1984 - https://nvlpubs.nist.gov/nistpubs/Legacy/IR/nbsir84-3007.pdf + [1] J. Simon, E. S. Drexler, and R. P. Reed, *NIST Monograph 177*, "Properties of + Copper and Copper Alloys at Cryogenic Temperatures", U.S. Government Printing Office, + February 1992. https://nvlpubs.nist.gov/nistpubs/Legacy/MONO/nistmonograph177.pdf + Equation 8-1 + + [2] J. G. Hust, A. B. Lankford, NBSIR 84-3007 "THERMAL CONDUCTIVITY OF ALUMINUM, + COPPER, IRON, AND TUNGSTEN FOR TEMPERATURES FROM 1 K TO THE MELTING POINT", 1984 + https://nvlpubs.nist.gov/nistpubs/Legacy/IR/nbsir84-3007.pdf """ p1: Final[float] = 1.171e-17 p2: Final[float] = 4.49 @@ -131,11 +130,12 @@ def _copper_irradiation_resistivity(fluence: float) -> float: References ---------- - - M. Kovari, 09/11/2012, internal notes (Excel / Mathcad), Technology Program, - WP12, PEX, Super-X Divertor for DEMO. - - M. Nakagawa et al., "High-dose neutron-irradiation effects in fcc metals at - 4.6 K", *Phys. Rev. B*, 16, 5285 (1977). - https://doi.org/10.1103/PhysRevB.16.5285 Figure 6 + [1] M. Kovari, 09/11/2012, internal notes (Excel / Mathcad), Technology Program, + WP12, PEX, Super-X Divertor for DEMO. + + [2] M. Nakagawa et al., "High-dose neutron-irradiation effects in fcc metals at + 4.6 K", *Phys. Rev. B*, 16, 5285 (1977). https://doi.org/10.1103/PhysRevB.16.5285 + Figure 6 """ c1: Final[float] = 0.00283 c2: Final[float] = -0.0711 @@ -170,11 +170,10 @@ def _copper_magneto_resistivity(resistivity: float, field: float) -> float: References ---------- - - J. Simon, E. S. Drexler, and R. P. Reed, *NIST Monograph 177*, "Properties of - Copper and Copper Alloys at Cryogenic Temperatures", - U.S. Government Printing Office, February 1992. - https://nvlpubs.nist.gov/nistpubs/Legacy/MONO/nistmonograph177.pdf - Equation 8-7 + [1] J. Simon, E. S. Drexler, and R. P. Reed, *NIST Monograph 177*, "Properties of + Copper and Copper Alloys at Cryogenic Temperatures", U.S. Government Printing Office, + February 1992. https://nvlpubs.nist.gov/nistpubs/Legacy/MONO/nistmonograph177.pdf + Equation 8-7 """ p9: Final[float] = 1.553e-8 @@ -194,8 +193,7 @@ def _copper_electrical_resistivity( temperature: float, field: float, rrr: float, fluence: float ) -> float: """Calculates the electrical resistivity of cryogenic copper with temperature, RRR, - magnetic field, and fluence dependence [Ω·m]. - + magnetic field, and fluence dependence [Ω·m]. Parameters ---------- @@ -219,22 +217,21 @@ def _copper_electrical_resistivity( References ---------- - - J. Simon, E. S. Drexler, and R. P. Reed, *NIST Monograph 177*, "Properties of - Copper and Copper Alloys at Cryogenic Temperatures", - U.S. Government Printing Office, February 1992. - https://nvlpubs.nist.gov/nistpubs/Legacy/MONO/nistmonograph177.pdf - Equation 8-1 - - - J. G. Hust, A. B. Lankford, NBSIR 84-3007 "THERMAL CONDUCTIVITY OF ALUMINUM, - COPPER, IRON, AND TUNGSTEN FOR TEMPERATURES FROM 1 K TO THE MELTING POINT", 1984 - https://nvlpubs.nist.gov/nistpubs/Legacy/IR/nbsir84-3007.pdf - - - M. Kovari, 09/11/2012, internal notes (Excel / Mathcad), Technology Program, - WP12, PEX, Super-X Divertor for DEMO. - - - M. Nakagawa et al., "High-dose neutron-irradiation effects in fcc metals at - 4.6 K", *Phys. Rev. B*, 16, 5285 (1977). https://doi.org/10.1103/PhysRevB.16.5285 - Figure 6 + [1] J. Simon, E. S. Drexler, and R. P. Reed, *NIST Monograph 177*, "Properties of + Copper and Copper Alloys at Cryogenic Temperatures", U.S. Government Printing Office, + February 1992. https://nvlpubs.nist.gov/nistpubs/Legacy/MONO/nistmonograph177.pdf + Equation 8-1 + + [2] J. G. Hust, A. B. Lankford, NBSIR 84-3007 "THERMAL CONDUCTIVITY OF ALUMINUM, + COPPER, IRON, AND TUNGSTEN FOR TEMPERATURES FROM 1 K TO THE MELTING POINT", 1984 + https://nvlpubs.nist.gov/nistpubs/Legacy/IR/nbsir84-3007.pdf + + [3] M. Kovari, 09/11/2012, internal notes (Excel / Mathcad), Technology Program, + WP12, PEX, Super-X Divertor for DEMO. + + [4] M. Nakagawa et al., "High-dose neutron-irradiation effects in fcc metals at + 4.6 K", *Phys. Rev. B*, 16, 5285 (1977). https://doi.org/10.1103/PhysRevB.16.5285 + Figure 6 """ rho_rrr = _copper_rrr_resistivity(temperature, rrr) rho_irr = _copper_irradiation_resistivity(fluence) @@ -264,16 +261,18 @@ def _nb3sn_specific_heat_capacity(temperature: float) -> float: References ---------- - - EFDA Material Data Compilation for Superconductor Simulation, P. Bauer, - H. Rajainmaki, E. Salpietro, EFDA CSU, Garching, 04/18/07. - - - ITER DRG1 Annex, Superconducting Material Database, Article 5, - N 11 FDR 42 01-07-05 R 0.1. - - V.D. Arp, Stability and Thermal Quenches in Force-Cooled Superconducting Cables, - Superconducting MHD Magnet Design Conf., MIT, pp 142-157, 1980. - - G.S. Knapp, S.D. Bader, Z. Fisk, Phonon properties of A-15 superconductors - obtained from heat capacity measurements, Phys. Rev. B, 13(9), pp 3783-3789, 1976. - https://doi.org/10.1103/PhysRevB.13.3783 + [1] EFDA Material Data Compilation for Superconductor Simulation, P. Bauer, + H. Rajainmaki, E. Salpietro, EFDA CSU, Garching, 04/18/07. + + [2] ITER DRG1 Annex, Superconducting Material Database, Article 5, + N 11 FDR 42 01-07-05 R 0.1. + + [3] V.D. Arp, Stability and Thermal Quenches in Force-Cooled Superconducting Cables, + Superconducting MHD Magnet Design Conf., MIT, pp 142-157, 1980. + + [4] G.S. Knapp, S.D. Bader, Z. Fisk, Phonon properties of A-15 superconductors + obtained from heat capacity measurements, Phys. Rev. B, 13(9), pp 3783-3789, 1976. + https://doi.org/10.1103/PhysRevB.13.3783 """ gamma: Final[float] = 0.1 # [J/K²/kg] (Grueneisen) beta: Final[float] = 0.001 # [J/K⁴/kg] (Debye) @@ -294,7 +293,7 @@ def _nb3sn_specific_heat_capacity(temperature: float) -> float: def _quench_integrals( t_he_peak: float, - t_max: float, + temp_quench_max: float, field: float, rrr: float, fluence: float, @@ -309,8 +308,8 @@ def _quench_integrals( ---------- t_he_peak: t_he_peak: Lower temperature bound of integration [K]. - t_max: - t_max: Upper temperature bound of integration [K]. + temp_quench_max: + temp_quench_max: Upper temperature bound of integration [K]. field: field: Magnetic field [T]. rrr: @@ -327,7 +326,7 @@ def _quench_integrals( Notes ----- - Integrals assume temperature-dependent material models are defined for the entire - range [t_he_peak, t_max]. + range [t_he_peak, temp_quench_max]. - Helium is assumed to be at constant pressure throughout the quench (i.e. some PRV in the cooling system) """ @@ -340,8 +339,8 @@ def _quench_integrals( isc = 0.0 for xi, wi in zip(GAUSS_LEG_NODES, GAUSS_LEG_WEIGHTS, strict=False): - ti = 0.5 * (xi + 1.0) * (t_max - t_he_peak) + t_he_peak - dti = 0.5 * wi * (t_max - t_he_peak) + ti = 0.5 * (xi + 1.0) * (temp_quench_max - t_he_peak) + t_he_peak + dti = 0.5 * wi * (temp_quench_max - t_he_peak) nu_cu = _copper_electrical_resistivity(ti, field, rrr, fluence) factor = dti / nu_cu @@ -357,13 +356,13 @@ def _quench_integrals( def calculate_quench_protection_current_density( tau_discharge: float, - peak_field: float, + b_peak: float, f_cu: float, f_he: float, t_he_peak: float, - t_max: float, + temp_quench_max: float, cu_rrr: float, - detection_time: float, + t_quench_detection: float, fluence: float, ) -> float: """Calculates the current density limited by the protection limit. @@ -380,20 +379,20 @@ def calculate_quench_protection_current_density( ---------- tau_discharge: tau_discharge: Quench discharge time constant [s]. - peak_field: - peak_field: Magnetic field at the peak point [T]. + b_peak: + b_peak: Magnetic field at the peak point [T]. f_cu: f_cu: Fraction of conductor cross-section that is copper. f_he: f_he: Fraction of cable occupied by helium. t_he_peak: t_he_peak: Peak helium temperature at quench initiation [K]. - t_max: - t_max: Maximum allowed conductor temperature during quench [K]. + temp_quench_max: + temp_quench_max: Maximum allowed conductor temperature during quench [K]. cu_rrr: cu_rrr: Residual resistivity ratio of copper. - detection_time: - detection_time: Detection time delay [s]. + t_quench_detection: + t_quench_detection: Detection time delay [s]. fluence: fluence: Neutron fluence [n/m²]. @@ -418,13 +417,297 @@ def calculate_quench_protection_current_density( warn("Fluence values out of range [0.0, 1.5e23]; kludging.", stacklevel=2) fluence = np.clip(fluence, 0.0, 1.5e23) - i_he, i_cu, i_sc = _quench_integrals(t_he_peak, t_max, peak_field, cu_rrr, fluence) + i_he, i_cu, i_sc = _quench_integrals( + t_he_peak, temp_quench_max, b_peak, cu_rrr, fluence + ) f_cond = 1.0 - f_he # Fraction of the cable XS area that is not helium f_cu_cable = f_cond * f_cu # Fraction of the cable XS that is copper f_sc_cable = f_cond * (1.0 - f_cu) # Fraction of the cable XS that is superconductor - factor = 1.0 / (0.5 * tau_discharge + detection_time) + factor = 1.0 / (0.5 * tau_discharge + t_quench_detection) total_integral = f_he * i_he + f_cu_cable * i_cu + f_sc_cable * i_sc return np.sqrt(factor * f_cu_cable * total_integral) + + +def plot_quench_time_evolution( + tau_discharge: float, + b_peak: float, + f_cu: float, + f_he: float, + t_he_peak: float, + temp_quench_max: float, + cu_rrr: float, + t_quench_detection: float, + fluence: float, + j_operating: float, + n_points: int = 500, + axes_1: plt.Axes = None, + axes_2: plt.Axes = None, + show: bool = False, +) -> None: + """Plots the time evolution of the quench model hotspot temperature and current. + + Visualises the adiabatic hotspot temperature rise and exponentially decaying + current during a quench, highlighting the quench detection time. + + Parameters + ---------- + tau_discharge: + Quench discharge time constant [s]. + b_peak: + Magnetic field at the peak point [T]. + f_cu: + Fraction of conductor cross-section that is copper. + f_he: + Fraction of cable occupied by helium. + t_he_peak: + Peak helium temperature at quench initiation [K]. + temp_quench_max: + Maximum allowed conductor temperature during quench [K]. + cu_rrr: + Residual resistivity ratio of copper. + t_quench_detection: + Detection time delay [s]. + fluence: + Neutron fluence [n/m²]. + j_operating: + Operating current density [A/m²] to compare against the quench protection limit. + n_points: + Number of time points for the plot. + axes_1: + Optional axis for the current density panel. + axes_2: + Optional axis for the hotspot temperature panel. + show: + Whether to display the plot with Matplotlib. Defaults to False to avoid + GUI backend warnings in non-interactive environments. + """ + figure = None + if axes_1 is None and axes_2 is None: + figure, (axes_1, axes_2) = plt.subplots(2, 1, sharex=True) + elif axes_1 is None or axes_2 is None: + msg = "Both axes_1 and axes_2 must be provided together, or neither." + raise ValueError(msg) + + fluence = np.clip(fluence, 0.0, 1.5e23) + + j_max = calculate_quench_protection_current_density( + tau_discharge, + b_peak, + f_cu, + f_he, + t_he_peak, + temp_quench_max, + cu_rrr, + t_quench_detection, + fluence, + ) + + fluence_1e23 = 1e23 + j_max_1e23 = calculate_quench_protection_current_density( + tau_discharge, + b_peak, + f_cu, + f_he, + t_he_peak, + temp_quench_max, + cu_rrr, + t_quench_detection, + fluence_1e23, + ) + + # Time axis: from 0 to ~4 time constants + t_end = 4.0 * tau_discharge + times = np.linspace(0.0, t_end, n_points) + + # Current density decays exponentially after detection + j_profile_required = np.where( + times < t_quench_detection, + j_max, + j_max * np.exp(-(times - t_quench_detection) / tau_discharge), + ) + + j_profile_required_1e23 = np.where( + times < t_quench_detection, + j_max_1e23, + j_max_1e23 * np.exp(-(times - t_quench_detection) / tau_discharge), + ) + + j_profile_real = np.where( + times < t_quench_detection, + j_operating, + j_operating * np.exp(-(times - t_quench_detection) / tau_discharge), + ) + + # Adiabatic hotspot temperature: integrate heat balance over time + # T(t) is found by inverting: integral_{T0}^{T(t)} [sum(rho*cp)] / rho_cu dT = integral_0^t J^2 dt + # We accumulate the "MIIT" (integral J^2 dt) and map it to temperature via the precomputed integral. + pressure = 6e5 + f_cond = 1.0 - f_he + f_cu_cable = f_cond * f_cu + f_sc_cable = f_cond * (1.0 - f_cu) + + def _build_cum_integral(fluence_val): + n_temp = 300 + temp_array = np.linspace(t_he_peak, temp_quench_max, n_temp) + cum_integral = np.zeros(n_temp) + for k in range(1, n_temp): + t_lo = temp_array[k - 1] + t_hi = temp_array[k] + val = 0.0 + for xi, wi in zip(GAUSS_LEG_NODES, GAUSS_LEG_WEIGHTS, strict=False): + ti = 0.5 * (xi + 1.0) * (t_hi - t_lo) + t_lo + dti = 0.5 * wi * (t_hi - t_lo) + nu_cu = _copper_electrical_resistivity(ti, b_peak, cu_rrr, fluence_val) + he_props = FluidProperties.of("He", temperature=ti, pressure=pressure) + integrand = ( + f_he * he_props.specific_heat_const_p * he_props.density + + f_cu_cable * _copper_specific_heat_capacity(ti) * COPPER_DENSITY + + f_sc_cable * _nb3sn_specific_heat_capacity(ti) * NB3SN_DENSITY + ) / nu_cu + val += dti * integrand + cum_integral[k] = cum_integral[k - 1] + val + return temp_array, cum_integral + + # Build a temperature lookup: cumulative integral from t_he_peak to T + temp_array, cum_integral = _build_cum_integral(fluence) + temp_array_1e23, cum_integral_1e23 = _build_cum_integral(fluence_1e23) + + # Numerically integrate J² dt over time to get MIIT at each time step + dt = times[1] - times[0] + miit_required = np.cumsum(j_profile_required**2) * dt + miit_required_1e23 = np.cumsum(j_profile_required_1e23**2) * dt + miit_real = np.cumsum(j_profile_real**2) * dt + + # Map MIIT -> temperature via inverse interpolation of cum_integral * f_cu_cable + scaled_integral = f_cu_cable * cum_integral + scaled_integral_1e23 = f_cu_cable * cum_integral_1e23 + hotspot_temp_required = np.interp(miit_required, scaled_integral, temp_array) + hotspot_temp_required_1e23 = np.interp( + miit_required_1e23, scaled_integral_1e23, temp_array_1e23 + ) + hotspot_temp_real = np.interp(miit_real, scaled_integral, temp_array) + + # --- Current density panel --- + axes_1.plot( + times, + j_profile_required, + color="darkorange", + linewidth=2, + label=f"Max allowed current density for protection (fluence = {fluence:.2e} n/m²)", + ) + axes_1.plot( + times, + j_profile_required_1e23, + color="darkorange", + linewidth=2, + linestyle="--", + label="Max allowed current density for protection (fluence = 1e23 n/m²)", + ) + axes_1.plot( + times, + j_profile_real, + color="blue", + linewidth=2, + label="Operating current density", + ) + axes_1.axvline( + t_quench_detection, + color="crimson", + linestyle="--", + linewidth=1.5, + label=f"Detection time ({t_quench_detection:.1f} s)", + ) + axes_1.axvspan( + 0, t_quench_detection, alpha=0.08, color="crimson", label="Pre-detection phase" + ) + axes_1.set_ylabel("Current density [A/m²]") + axes_1.legend(fontsize=9) + axes_1.grid(True, alpha=0.3) + axes_1.set_title( + "TF Coil Quench Protection: Current Density and Hotspot Temperature Evolution" + ) + + # --- Temperature panel --- + axes_2.plot( + times, + hotspot_temp_required, + color="darkorange", + linewidth=2, + label=f"Hotspot temperature at protection limit (fluence = {fluence:.2e} n/m²)", + ) + axes_2.plot( + times, + hotspot_temp_required_1e23, + color="darkorange", + linewidth=2, + linestyle="--", + label="Hotspot temperature at protection limit (fluence = 1e23 n/m²)", + ) + axes_2.plot( + times, + hotspot_temp_real, + color="blue", + linewidth=2, + label="Operating hotspot temperature", + ) + + axes_2.axvline( + t_quench_detection, + color="crimson", + linestyle="--", + linewidth=1.5, + label=f"$t_{{\\text{{detect}}}}$ ({t_quench_detection:.2f} s)", + ) + axes_2.axvspan(0, t_quench_detection, alpha=0.08, color="crimson") + axes_2.axhline( + temp_quench_max, + color="grey", + linestyle=":", + linewidth=1.5, + label=f"$T_{{\\text{{max}}}}$ = {temp_quench_max} K", + ) + axes_2.set_xlabel("Time [s]") + axes_2.set_ylabel("Temperature [K]") + axes_2.legend(fontsize=9) + axes_2.grid(True, alpha=0.3) + + # Mark tau_discharge after detection time with vertical and horizontal lines + tau_time = t_quench_detection + tau_discharge + tau_j = j_max * np.exp( + -1 + ) # current density at t = t_quench_detection + tau_discharge + tau_temp = float(np.interp(tau_time, times, hotspot_temp_required)) + + for ax, val, label in [ + (axes_1, tau_j, f"$J$ at $\\tau_{{\\text{{discharge}}}}$ ({tau_j:.2e} A/m²)"), + (axes_2, tau_temp, f"$T$ at $\\tau_{{\\text{{discharge}}}}$ ({tau_temp:.1f} K)"), + ]: + ax.axvline( + tau_time, + color="forestgreen", + linestyle="--", + linewidth=1.5, + label=f"$t_{{\\text{{detect}}}} + \\tau_{{\\text{{discharge}}}}$ ({tau_time:.2f} s)", + ) + ax.axhline( + val, + color="forestgreen", + linestyle=":", + linewidth=1.5, + label=label, + ) + axes_1.legend(fontsize=9) + axes_1.minorticks_on() + axes_2.legend(fontsize=9) + axes_2.minorticks_on() + + if figure is not None: + figure.tight_layout() + else: + plt.tight_layout() + + if show: + plt.show() diff --git a/process/models/tfcoil/superconducting.py b/process/models/tfcoil/superconducting.py index 1a9d7283d..98f51e42e 100644 --- a/process/models/tfcoil/superconducting.py +++ b/process/models/tfcoil/superconducting.py @@ -594,6 +594,20 @@ def output_tf_superconductor_info(self): tfcoil_variables.j_tf_wp_quench_heat_max, "OP ", ) + po.ovarre( + self.outfile, + "Max allowed fast neutron fluence on TF coil (n/m²)", + "(nflutfmax)", + self.data.constraints.nflutfmax, + "OP ", + ) + po.ovarre( + self.outfile, + "Residual resistivity ratio of TF coil copper", + "(rrr_tf_cu)", + tfcoil_variables.rrr_tf_cu, + "OP ", + ) @staticmethod def calculate_superconductor_temperature_margin(