diff --git a/CMakeLists.txt b/CMakeLists.txt index ec1238130c..316f46cde4 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -81,7 +81,6 @@ LIST(APPEND PROCESS_SRCS reinke_module.f90 sctfcoil.f90 plasma_profiles.f90 - superconductors.f90 final_module.f90 cost_variables.f90 divertor_variables.f90 diff --git a/process/pfcoil.py b/process/pfcoil.py index b0495cabef..85e26539c4 100644 --- a/process/pfcoil.py +++ b/process/pfcoil.py @@ -11,12 +11,12 @@ from process.fortran import maths_library as ml from process.fortran import process_output as op from process.fortran import numerics -from process.fortran import superconductors as sc from process.fortran import rebco_variables as rcv from process.fortran import constraint_variables as ctv from process import maths_library as pml from process.utilities.f2py_string_patch import f2py_compatible_to_string from process import fortran as ft +import process.superconductors as superconductors import math import numpy as np import numba @@ -175,7 +175,6 @@ def pfcoil(self): # N.B. Problems here if k=ncls(group) is greater than 2. for j in range(pfv.ngrp): - if pfv.ipfloc[j] == 1: # PF coil is stacked on top of the Central Solenoid for k in pfv.ncls[j]: @@ -333,7 +332,6 @@ def pfcoil(self): ) else: - # Conventional aspect ratio scaling nfxf0 = 0 ngrp0 = 0 @@ -554,7 +552,6 @@ def pfcoil(self): for ii in range(pfv.ngrp): for ij in range(pfv.ncls[ii]): - if pfv.ipfloc[ii] == 1: # PF coil is stacked on top of the Central Solenoid dx = 0.5e0 * bv.ohcth @@ -588,7 +585,6 @@ def pfcoil(self): pfv.zh[i] = pfv.zpf[i] - dz else: - # Other coils. N.B. Current density RJCONPF[i] is defined in # routine INITIAL for these coils. area = abs(pfv.ric[i] * 1.0e6 / pfv.rjconpf[i]) @@ -624,7 +620,6 @@ def pfcoil(self): for ii in range(pfv.ngrp): iii = ii for ij in range(pfv.ncls[ii]): - # Peak field if ij == 0: @@ -1023,7 +1018,6 @@ def ohcalc(self): sgn = 1.0e0 pfv.ric[pfv.nohc - 1] = sgn * 1.0e-6 * pfv.cohbop * pfv.areaoh else: - sgn = -1.0e0 pfv.ric[pfv.nohc - 1] = sgn * 1.0e-6 * pfv.coheof * pfv.areaoh @@ -1104,7 +1098,6 @@ def ohcalc(self): # Stress ==> cross-sectional area of supporting steel to use if pfv.ipfres == 0: - # Superconducting coil # New calculation from M. N. Wilson for hoop stress @@ -1192,7 +1185,6 @@ def ohcalc(self): ) if pfv.ipfres == 0: - # Allowable coil overall current density at EOF # (superconducting coils only) @@ -1306,7 +1298,6 @@ def peakb(self, i, ii, it): jj = 0 for iii in range(pfv.ngrp): for jjj in range(pfv.ncls[iii]): - jj = jj + 1 # Radius, z-coordinate and current for each coil if iii == ii - 1: @@ -1419,32 +1410,27 @@ def bfmax(self, rj, a, b, h): ) if beta > 3.0: - b1 = constants.rmu0 * rj * (b - a) f = (3.0 / beta) ** 2 bfmax = f * b0 * (1.007 + (alpha - 1.0) * 0.0055) + (1.0 - f) * b1 elif beta > 2.0: - rat = (1.025 - (beta - 2.0) * 0.018) + (alpha - 1.0) * ( 0.01 - (beta - 2.0) * 0.0045 ) bfmax = rat * b0 elif beta > 1.0: - rat = (1.117 - (beta - 1.0) * 0.092) + (alpha - 1.0) * (beta - 1.0) * 0.01 bfmax = rat * b0 elif beta > 0.75: - rat = (1.30 - 0.732 * (beta - 0.75)) + (alpha - 1.0) * ( 0.2 * (beta - 0.75) - 0.05 ) bfmax = rat * b0 else: - rat = (1.65 - 1.4 * (beta - 0.5)) + (alpha - 1.0) * ( 0.6 * (beta - 0.5) - 0.20 ) @@ -1769,7 +1755,6 @@ def induct(self, output): ] if bv.iohcl != 0: - # Central Solenoid self inductance a = pfv.rohc # mean radius of coil b = 2.0e0 * pfv.zh[pfv.nohc - 1] # length of coil @@ -2704,7 +2689,6 @@ def waveform(self): pfv.waves[nplas - 1, it] = 1.0e0 for ic in range(pfv.nohc): - # Find where the peak current occurs # Beginning of pulse, t = tramp if (abs(pfv.curpfs[ic]) >= abs(pfv.curpfb[ic])) and ( @@ -2794,7 +2778,7 @@ def deltaj_nbti(temperature): :return: difference in current density :rtype: float """ - jcrit0, t = sc.jcrit_nbti(temperature, bmax, c0, bc20m, tc0m) + jcrit0, _ = superconductors.jcrit_nbti(temperature, bmax, c0, bc20m, tc0m) if ml.variable_error(jcrit0): # superconductors.jcrit_nbti has failed. print(f"superconductors.jcrit_nbti: {bmax=} {temperature=} {jcrit0=}") @@ -2809,7 +2793,7 @@ def deltaj_wst(temperature): :return: difference in current density :rtype: float """ - jcrit0, b, t = sc.wstsc(temperature, bmax, strain, bc20m, tc0m) + jcrit0, _, _ = superconductors.wstsc(temperature, bmax, strain, bc20m, tc0m) if ml.variable_error(jcrit0): # superconductors.wstsc has failed. print(f"deltaj_wst: {bmax=} {temperature=} {jcrit0=}") @@ -2824,7 +2808,9 @@ def deltaj_gl_nbti(temperature): :return: difference in current density :rtype: float """ - jcrit0, b, t = sc.gl_nbti(temperature, bmax, strain, bc20m, tc0m) + jcrit0, _, _ = superconductors.gl_nbti( + temperature, bmax, strain, bc20m, tc0m + ) if ml.variable_error(jcrit0): # GL_Nbti has failed. print(f"deltaj_GL_nbti: {bmax=} {temperature=} {jcrit0=}") @@ -2839,7 +2825,9 @@ def deltaj_gl_rebco(temperature): :return: difference in current density :rtype: float """ - jcrit0, b, t = sc.gl_rebco(temperature, bmax, strain, bc20m, tc0m) + jcrit0, _, _ = superconductors.gl_rebco( + temperature, bmax, strain, bc20m, tc0m + ) if ml.variable_error(jcrit0): # superconductors.GL_REBCO has failed. print(f"deltaj_gl_REBCO: {bmax=} {temperature=} {jcrit0=}") @@ -2854,7 +2842,9 @@ def deltaj_hijc_rebco(temperature): :return: difference in current density :rtype: float """ - jcrit0, b, t = sc.hijc_rebco(temperature, bmax, strain, bc20m, tc0m) + jcrit0, _, _ = superconductors.hijc_rebco( + temperature, bmax, strain, bc20m, tc0m + ) if ml.variable_error(jcrit0): # superconductors.GL_REBCO has failed. print(f"deltaj_hijc_REBCO: {bmax=} {temperature=} {jcrit0=}") @@ -2870,7 +2860,7 @@ def deltaj_hijc_rebco(temperature): # jcritsc returned by superconductors.itersc is the critical current density in the # superconductor - not the whole strand, which contains copper - jcritsc, bcrit, tcrit = sc.itersc(thelium, bmax, strain, bc20m, tc0m) + jcritsc, _, _ = superconductors.itersc(thelium, bmax, strain, bc20m, tc0m) jcritstr = jcritsc * (1.0e0 - fcu) elif isumat == 2: @@ -2885,23 +2875,22 @@ def deltaj_hijc_rebco(temperature): jstrand = jwp / (1.0e0 - fhe) - jcritstr, tmarg = sc.bi2212(bmax, jstrand, thelium, fhts) + jcritstr, tmarg = superconductors.bi2212(bmax, jstrand, thelium, fhts) jcritsc = jcritstr / (1.0e0 - fcu) - tcrit = thelium + tmarg elif isumat == 3: # NbTi data bc20m = 15.0e0 tc0m = 9.3e0 c0 = 1.0e10 - jcritsc, tcrit = sc.jcrit_nbti(thelium, bmax, c0, bc20m, tc0m) + jcritsc, _ = superconductors.jcrit_nbti(thelium, bmax, c0, bc20m, tc0m) jcritstr = jcritsc * (1.0e0 - fcu) elif isumat == 4: # As (1), but user-defined parameters bc20m = bcritsc tc0m = tcritsc - jcritsc, bcrit, tcrit = sc.itersc(thelium, bmax, strain, bc20m, tc0m) + jcritsc, _, _ = superconductors.itersc(thelium, bmax, strain, bc20m, tc0m) jcritstr = jcritsc * (1.0e0 - fcu) elif isumat == 5: @@ -2912,12 +2901,12 @@ def deltaj_hijc_rebco(temperature): # jcritsc returned by superconductors.itersc is the critical current density in the # superconductor - not the whole strand, which contains copper - jcritsc, bcrit, tcrit = sc.wstsc(thelium, bmax, strain, bc20m, tc0m) + jcritsc, _, _ = superconductors.wstsc(thelium, bmax, strain, bc20m, tc0m) jcritstr = jcritsc * (1.0e0 - fcu) elif isumat == 6: # "REBCO" 2nd generation HTS superconductor in CrCo strand - jcritsc, validity = sc.jcrit_rebco(thelium, bmax, 0) + jcritsc, _ = superconductors.jcrit_rebco(thelium, bmax) jcritstr = jcritsc * (1.0e0 - fcu) # The CS coil current at EOF @@ -2931,7 +2920,7 @@ def deltaj_hijc_rebco(temperature): # Durham Ginzburg-Landau critical surface model for Nb-Ti bc20m = tfv.b_crit_upper_nbti tc0m = tfv.t_crit_nbti - jcritsc, bcrit, tcrit = sc.gl_nbti(thelium, bmax, strain, bc20m, tc0m) + jcritsc, _, _ = superconductors.gl_nbti(thelium, bmax, strain, bc20m, tc0m) jcritstr = jcritsc * (1.0e0 - fcu) # The CS coil current at EOF @@ -2941,7 +2930,7 @@ def deltaj_hijc_rebco(temperature): # Durham Ginzburg-Landau critical surface model for REBCO bc20m = 429e0 tc0m = 185e0 - jcritsc, bcrit, tcrit = sc.gl_rebco(thelium, bmax, strain, bc20m, tc0m) + jcritsc, _, _ = superconductors.gl_rebco(thelium, bmax, strain, bc20m, tc0m) # A0 calculated for tape cross section already jcritstr = jcritsc * (1.0e0 - fcu) @@ -2954,7 +2943,9 @@ def deltaj_hijc_rebco(temperature): # Hazelton experimental data + Zhai conceptual model for REBCO bc20m = 138 tc0m = 92 - jcritsc, bcrit, tcrit = sc.hijc_rebco(thelium, bmax, strain, bc20m, tc0m) + jcritsc, _, _ = superconductors.hijc_rebco( + thelium, bmax, strain, bc20m, tc0m + ) # A0 calculated for tape cross section already jcritstr = jcritsc * (1.0e0 - fcu) @@ -2995,14 +2986,20 @@ def deltaj_hijc_rebco(temperature): ttestp = ttest + delt if isumat in [1, 4]: - jcrit0, b, t = sc.itersc(ttest, bmax, strain, bc20m, tc0m) + jcrit0, _, _ = superconductors.itersc( + ttest, bmax, strain, bc20m, tc0m + ) if (abs(jsc - jcrit0) <= jtol) and ( abs((jsc - jcrit0) / jsc) <= 0.01 ): break - jcritm, b, t = sc.itersc(ttestm, bmax, strain, bc20m, tc0m) - jcritp, b, t = sc.itersc(ttestp, bmax, strain, bc20m, tc0m) + jcritm, _, _ = superconductors.itersc( + ttestm, bmax, strain, bc20m, tc0m + ) + jcritp, _, _ = superconductors.itersc( + ttestp, bmax, strain, bc20m, tc0m + ) # Kludge to avoid divide by 0 if jcritm == jcritp: @@ -3021,7 +3018,9 @@ def deltaj_hijc_rebco(temperature): deltaj_nbti, x1, x2, 100e0 ) tmarg = current_sharing_t - thelium - jcrit0, t = sc.jcrit_nbti(current_sharing_t, bmax, c0, bc20m, tc0m) + jcrit0, _ = superconductors.jcrit_nbti( + current_sharing_t, bmax, c0, bc20m, tc0m + ) if ml.variable_error( current_sharing_t ): # current sharing secant solver has failed. @@ -3039,7 +3038,9 @@ def deltaj_hijc_rebco(temperature): deltaj_wst, x1, x2, 100e0 ) tmarg = current_sharing_t - thelium - jcrit0, b, t = sc.wstsc(current_sharing_t, bmax, strain, bc20m, tc0m) + jcrit0, _, _ = superconductors.wstsc( + current_sharing_t, bmax, strain, bc20m, tc0m + ) if ml.variable_error( current_sharing_t ): # current sharing secant solver has failed. @@ -3049,7 +3050,7 @@ def deltaj_hijc_rebco(temperature): # Temperature margin: An alternative method using secant solver elif isumat == 6: - current_sharing_t = sc.current_sharing_rebco(bmax, jsc) + current_sharing_t = superconductors.current_sharing_rebco(bmax, jsc) tmarg = current_sharing_t - thelium tfv.temp_margin = tmarg @@ -3063,7 +3064,9 @@ def deltaj_hijc_rebco(temperature): deltaj_gl_nbti, x1, x2, 100e0 ) tmarg = current_sharing_t - thelium - jcrit0, b, t = sc.gl_nbti(current_sharing_t, bmax, strain, bc20m, tc0m) + jcrit0, _, _ = superconductors.gl_nbti( + current_sharing_t, bmax, strain, bc20m, tc0m + ) if ml.variable_error( current_sharing_t ): # current sharing secant solver has failed. @@ -3081,7 +3084,9 @@ def deltaj_hijc_rebco(temperature): deltaj_gl_rebco, x1, x2, 100e0 ) tmarg = current_sharing_t - thelium - jcrit0, b, t = sc.gl_rebco(current_sharing_t, bmax, strain, bc20m, tc0m) + jcrit0, _, _ = superconductors.gl_rebco( + current_sharing_t, bmax, strain, bc20m, tc0m + ) if ml.variable_error( current_sharing_t ): # current sharing secant solver has failed. @@ -3098,7 +3103,9 @@ def deltaj_hijc_rebco(temperature): deltaj_hijc_rebco, x1, x2, 100e0 ) tmarg = current_sharing_t - thelium - jcrit0, b, t = sc.hijc_rebco(current_sharing_t, bmax, strain, bc20m, tc0m) + jcrit0, _, _ = superconductors.hijc_rebco( + current_sharing_t, bmax, strain, bc20m, tc0m + ) if ml.variable_error( current_sharing_t ): # current sharing secant solver has failed. diff --git a/process/sctfcoil.py b/process/sctfcoil.py index eb95fe41ed..f4a7854068 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -6,7 +6,6 @@ from process.fortran import rebco_variables from process.fortran import global_variables -from process.fortran import superconductors from process.fortran import tfcoil_variables from process.fortran import physics_variables from process.fortran import build_variables @@ -19,6 +18,8 @@ from process.fortran import numerics from process.fortran import divertor_variables +import process.superconductors as superconductors + from process.utilities.f2py_string_patch import f2py_compatible_to_string @@ -118,7 +119,7 @@ def supercon_croco(self, aturn, bmax, iop, thelium, output: bool): jcritsc: float = 0.0 # Find critical current density in superconducting strand, jcritstr - jcritsc, _ = superconductors.jcrit_rebco(thelium, bmax, int(output)) + jcritsc, _ = superconductors.jcrit_rebco(thelium, bmax) # tfcoil_variables.acstf : Cable space - inside area (m2) # Set new rebco_variables.croco_od # allowing for scaling of rebco_variables.croco_od @@ -142,8 +143,7 @@ def supercon_croco(self, aturn, bmax, iop, thelium, output: bool): sctfcoil_module.conductor_jacket_fraction = ( sctfcoil_module.conductor_jacket_area / sctfcoil_module.conductor_area ) - superconductors.croco( - jcritsc, + ( sctfcoil_module.croco_strand_area, sctfcoil_module.croco_strand_critical_current, sctfcoil_module.conductor_copper_area, @@ -158,10 +158,13 @@ def supercon_croco(self, aturn, bmax, iop, thelium, output: bool): sctfcoil_module.conductor_rebco_area, sctfcoil_module.conductor_rebco_fraction, sctfcoil_module.conductor_critical_current, + ) = superconductors.croco( + jcritsc, sctfcoil_module.conductor_area, rebco_variables.croco_od, rebco_variables.croco_thick, ) + rebco_variables.coppera_m2 = iop / sctfcoil_module.conductor_copper_area icrit = sctfcoil_module.conductor_critical_current @@ -645,9 +648,7 @@ def supercon( # jcritsc returned by superconductors.itersc is the critical current density in the # superconductor - not the whole strand, which contains copper - jcritsc, bcrit, tcrit = superconductors.itersc( - thelium, bmax, strain, bc20m, tc0m - ) + jcritsc, _, _ = superconductors.itersc(thelium, bmax, strain, bc20m, tc0m) jcritstr = jcritsc * (1.0e0 - fcu) # Critical current in cable icrit = jcritstr * acs * fcond @@ -663,7 +664,6 @@ def supercon( jcritstr, tmarg = superconductors.bi2212(bmax, jstrand, thelium, fhts) jcritsc = jcritstr / (1.0e0 - fcu) - tcrit = thelium + tmarg # Critical current in cable icrit = jcritstr * acs * fcond @@ -671,7 +671,7 @@ def supercon( bc20m = 15.0e0 tc0m = 9.3e0 c0 = 1.0e10 - jcritsc, tcrit = superconductors.jcrit_nbti(thelium, bmax, c0, bc20m, tc0m) + jcritsc, _ = superconductors.jcrit_nbti(thelium, bmax, c0, bc20m, tc0m) jcritstr = jcritsc * (1.0e0 - fcu) # Critical current in cable icrit = jcritstr * acs * fcond @@ -685,9 +685,7 @@ def supercon( error_handling.report_error(261) strain = numpy.sign(strain) * 0.5e-2 - jcritsc, bcrit, tcrit = superconductors.itersc( - thelium, bmax, strain, bc20m, tc0m - ) + jcritsc, _, _ = superconductors.itersc(thelium, bmax, strain, bc20m, tc0m) jcritstr = jcritsc * (1.0e0 - fcu) # Critical current in cable icrit = jcritstr * acs * fcond @@ -703,9 +701,7 @@ def supercon( # jcritsc returned by superconductors.itersc is the critical current density in the # superconductor - not the whole strand, which contains copper - jcritsc, bcrit, tcrit = superconductors.wstsc( - thelium, bmax, strain, bc20m, tc0m - ) + jcritsc, _, _ = superconductors.wstsc(thelium, bmax, strain, bc20m, tc0m) jcritstr = jcritsc * (1.0e0 - fcu) # Critical current in cable icrit = jcritstr * acs * fcond @@ -718,9 +714,7 @@ def supercon( elif isumat == 7: # Durham Ginzburg-Landau Nb-Ti parameterisation bc20m = tfcoil_variables.b_crit_upper_nbti tc0m = tfcoil_variables.t_crit_nbti - jcritsc, bcrit, tcrit = superconductors.gl_nbti( - thelium, bmax, strain, bc20m, tc0m - ) + jcritsc, _, _ = superconductors.gl_nbti(thelium, bmax, strain, bc20m, tc0m) jcritstr = jcritsc * (1.0e0 - fcu) # Critical current in cable icrit = jcritstr * acs * fcond @@ -734,9 +728,7 @@ def supercon( error_handling.report_error(261) strain = numpy.sign(strain) * 0.7e-2 - jcritsc, bcrit, tcrit = superconductors.gl_rebco( - thelium, bmax, strain, bc20m, tc0m - ) + jcritsc, _, _ = superconductors.gl_rebco(thelium, bmax, strain, bc20m, tc0m) # A0 calculated for tape cross section already jcritstr = jcritsc * (1.0e0 - fcu) # Critical current in cable (copper added at this stage in HTS cables) @@ -756,7 +748,7 @@ def supercon( # 'high current density' as per parameterisation described in Wolf, # and based on Hazelton experimental data and Zhai conceptual model; # see subroutine for full references - jcritsc, bcrit, tcrit = superconductors.hijc_rebco( + jcritsc, _, _ = superconductors.hijc_rebco( thelium, bmax, strain, bc20m, tc0m ) @@ -821,85 +813,83 @@ def supercon( # Issue #483 to be on the safe side, check the fractional as well as the absolute error if isumat in (1, 4): - jcrit0, b, t = superconductors.itersc( + jcrit0, _, _ = superconductors.itersc( ttest, bmax, strain, bc20m, tc0m ) if (abs(jsc - jcrit0) <= jtol) and ( abs((jsc - jcrit0) / jsc) <= 0.01 ): break - jcritm, b, t = superconductors.itersc( + jcritm, _, _ = superconductors.itersc( ttestm, bmax, strain, bc20m, tc0m ) - jcritp, b, t = superconductors.itersc( + jcritp, _, _ = superconductors.itersc( ttestp, bmax, strain, bc20m, tc0m ) elif isumat == 3: - jcrit0, t = superconductors.jcrit_nbti(ttest, bmax, c0, bc20m, tc0m) + jcrit0, _ = superconductors.jcrit_nbti(ttest, bmax, c0, bc20m, tc0m) if (abs(jsc - jcrit0) <= jtol) and ( abs((jsc - jcrit0) / jsc) <= 0.01 ): break - jcritm, t = superconductors.jcrit_nbti( - ttestm, bmax, c0, bc20m, tc0m - ) - jcritp, t = superconductors.jcrit_nbti( + jcritm, _ = super.jcrit_nbti(ttestm, bmax, c0, bc20m, tc0m) + jcritp, _ = superconductors.jcrit_nbti( ttestp, bmax, c0, bc20m, tc0m ) elif isumat == 5: - jcrit0, b, t = superconductors.wstsc( + jcrit0, _, _ = superconductors.wstsc( ttest, bmax, strain, bc20m, tc0m ) if (abs(jsc - jcrit0) <= jtol) and ( abs((jsc - jcrit0) / jsc) <= 0.01 ): break - jcritm, b, t = superconductors.wstsc( + jcritm, _, _ = superconductors.wstsc( ttestm, bmax, strain, bc20m, tc0m ) - jcritp, b, t = superconductors.wstsc( + jcritp, _, _ = superconductors.wstsc( ttestp, bmax, strain, bc20m, tc0m ) elif isumat == 7: - jcrit0, b, t = superconductors.gl_nbti( + jcrit0, _, _ = superconductors.gl_nbti( ttest, bmax, strain, bc20m, tc0m ) if (abs(jsc - jcrit0) <= jtol) and ( abs((jsc - jcrit0) / jsc) <= 0.01 ): break - jcritm, b, t = superconductors.gl_nbti( + jcritm, _, _ = superconductors.gl_nbti( ttestm, bmax, strain, bc20m, tc0m ) - jcritp, b, t = superconductors.gl_nbti( + jcritp, _, _ = superconductors.gl_nbti( ttestp, bmax, strain, bc20m, tc0m ) elif isumat == 8: - jcrit0, b, t = superconductors.gl_rebco( + jcrit0, _, _ = superconductors.gl_rebco( ttest, bmax, strain, bc20m, tc0m ) if (abs(jsc - jcrit0) <= jtol) and ( abs((jsc - jcrit0) / jsc) <= 0.01 ): break - jcritm, b, t = superconductors.gl_rebco( + jcritm, _, _ = superconductors.gl_rebco( ttestm, bmax, strain, bc20m, tc0m ) - jcritp, b, t = superconductors.gl_rebco( + jcritp, _, _ = superconductors.gl_rebco( ttestp, bmax, strain, bc20m, tc0m ) elif isumat == 9: - jcrit0, b, t = superconductors.hijc_rebco( + jcrit0, _, _ = superconductors.hijc_rebco( ttest, bmax, strain, bc20m, tc0m ) if (abs(jsc - jcrit0) <= jtol) and ( abs((jsc - jcrit0) / jsc) <= 0.01 ): break - jcritm, b, t = superconductors.hijc_rebco( + jcritm, _, _ = superconductors.hijc_rebco( ttestm, bmax, strain, bc20m, tc0m ) - jcritp, b, t = superconductors.hijc_rebco( + jcritp, _, _ = superconductors.hijc_rebco( ttestp, bmax, strain, bc20m, tc0m ) diff --git a/process/stellarator.py b/process/stellarator.py index e359cffc3b..6edefb4faf 100644 --- a/process/stellarator.py +++ b/process/stellarator.py @@ -25,13 +25,13 @@ constraint_variables, rebco_variables, maths_library, - superconductors, profiles_module, physics_functions_module, neoclassics_module, impurity_radiation_module, current_drive_module, ) +import process.superconductors as superconductors import process.physics_functions as physics_funcs from process.coolprop_interface import FluidProperties diff --git a/process/superconductors.py b/process/superconductors.py new file mode 100644 index 0000000000..964914befb --- /dev/null +++ b/process/superconductors.py @@ -0,0 +1,707 @@ +import logging +import numpy as np + +from process.fortran import error_handling as eh, rebco_variables +from process.maths_library import secant_solve + +logger = logging.getLogger(__name__) + + +def jcrit_rebco(temperature, b): + """Critical current density for "REBCO" 2nd generation HTS superconductor + temperature : input real : superconductor temperature (K) + b : input real : Magnetic field at superconductor (T) + jcrit : output real : Critical current density in superconductor (A/m2) + + Will return a negative number if the temperature is greater than Tc0, the + zero-field critical temperature. + """ + tc0 = 90.0 # (K) + birr0 = 132.5 # (T) + a = 1.82962e8 # scaling constant + # exponents + p = 0.5875 + q = 1.7 + alpha = 1.54121 + beta = 1.96679 + oneoveralpha = 1 / alpha + + validity = True + + if (temperature < 4.2) or (temperature > 72.0): + validity = False + if temperature < 65: + if (b < 0.0) or (b > 15.0): + validity = False + else: + if (b < 0.0) or (b > 11.5): + validity = False + + if not validity: + logger.warning( + f"""jcrit_rebco: input out of range + temperature: {temperature} + Field: {b} + """ + ) + + if temperature < tc0: + # Normal case + birr = birr0 * (1 - temperature / tc0) ** alpha + else: + # If temp is greater than critical temp, ensure result is real but negative. + birr = birr0 * (1 - temperature / tc0) + + if b < birr: + # Normal case + factor = (b / birr) ** p * (1 - b / birr) ** q + jcrit = (a / b) * (birr**beta) * factor + else: + # Field is too high + # Ensure result is real but negative, and varies with temperature. + # tcb = critical temperature at field b + tcb = tc0 * (1 - (b / birr0) ** oneoveralpha) + jcrit = -(temperature - tcb) + + return jcrit, validity + + +def current_sharing_rebco(bfield, j): + """Current sharing temperature for "REBCO" 2nd generation HTS superconductor + b : input real : Magnetic field at superconductor (T) + j : input real : Current density in superconductor (A/m2) + current_sharing_t : output real : Current sharing temperature (K) + """ + + def deltaj_rebco(temperature): + jcritical, _ = jcrit_rebco(temperature, bfield) + return jcritical - j + + current_sharing_t, _, _ = secant_solve(deltaj_rebco, 4, 20, 1e7) + return current_sharing_t + + +def itersc(thelium, bmax, strain, bc20max, tc0max): + """Implementation of ITER Nb3Sn critical surface implementation + author: R Kemp, CCFE, Culham Science Centre + author: P J Knight, CCFE, Culham Science Centre + thelium : input real : Coolant/SC temperature (K) + bmax : input real : Magnetic field at conductor (T) + strain : input real : Strain in superconductor + bc20max : input real : Upper critical field (T) for superconductor + at zero temperature and strain + tc0max : input real : Critical temperature (K) at zero field and strain + jcrit : output real : Critical current density in superconductor (A/m2) + bcrit : output real : Critical field (T) + tcrit : output real : Critical temperature (K) + This routine calculates the critical current density and + temperature in the superconducting TF coils using the + ITER Nb3Sn critical surface model. + $J_C(B,T,\\epsilon)$ Parameterization for ITER Nb3Sn production, + L. Bottura, CERN-ITER Collaboration Report, Version 2, April 2nd 2008 + (distributed by Arnaud Devred, ITER, 10th April 2008) + ITER Nb3Sn critical surface parameterization (2MMF7J) (2008), + https://user.iter.org/?uid=2MMF7J&action=get_document + ITER DDD 11-7: Magnets - conductors (2NBKXY) (2009), + https://user.iter.org/?uid=2NBKXY&action=get_document + """ + + csc = 19922.0e6 # scaling constant C + p = 0.63 # low field exponent p + q = 2.1 # high field exponent q + ca1 = 44.48 # strain fitting constant C_{a1} + ca2 = 0.0 # strain fitting constant C_{a2} + eps0a = 0.00256 # epsilon_{0,a} + diter = 0.82 # ITER strand diameter (mm) + cuiter = 0.5 # ITER strand copper fraction + + # $\epsilon_{sh}$ + epssh = (ca2 * eps0a) / (np.sqrt(ca1**2 - ca2**2)) + + # Strain function $s(\epsilon)$ + # 0.83 < s < 1.0, for -0.005 < strain < 0.005 + strfun = np.sqrt(epssh**2 + eps0a**2) - np.sqrt( + (strain - epssh) ** 2 + eps0a**2 + ) + strfun = strfun * ca1 - ca2 * strain + strfun = 1.0 + (1 / (1.0 - ca1 * eps0a)) * strfun + + # $B^*_{C2} (0,\epsilon)$ + bc20eps = bc20max * strfun + + # $T^*_C (0,\epsilon)$ + tc0eps = tc0max * strfun ** (1 / 3) + + # Reduced temperature, restricted to be < 1 + # Should remain < 1 for thelium < 0.94*tc0max (i.e. 15 kelvin for i_tf_sc_mat=1) + + if thelium / tc0eps >= 1.0: + eh.fdiags[0] = thelium + eh.fdiags[1] = tc0eps + eh.report_error(159) + + t = min(thelium / tc0eps, 0.9999) + + # Reduced magnetic field at zero temperature + # Should remain < 1 for bmax < 0.83*bc20max (i.e. 27 tesla for i_tf_sc_mat=1) + + if bmax / bc20eps >= 1.0: + eh.fdiags[0] = bmax + eh.fdiags[1] = bc20eps + eh.report_error(160) + + bzero = min(bmax / bc20eps, 0.9999) + + # Critical temperature (K) + tcrit = tc0eps * (1.0 - bzero) ** (1 / 1.52) # bzero must be < 1 to avoid NaNs + + # Critical field (T) + bcrit = bc20eps * (1.0 - t**1.52) + + # Reduced magnetic field, restricted to be < 1 + if bmax / bcrit >= 1.0: + eh.fdiags[0] = bmax + eh.fdiags[1] = bcrit + eh.report_error(161) + + bred = min(bmax / bcrit, 0.9999) + + # Critical current density in superconductor (A/m2) + # ITER parameterization is for the current in a single strand, + # not per unit area, so scalefac converts to current density + + scalefac = np.pi * (0.5 * diter) ** 2 * (1.0 - cuiter) + + jc1 = (csc / bmax) * strfun + jc2 = (1.0 - t**1.52) * (1.0 - t**2) # t must be < 1 to avoid NaNs + jc3 = bred**p * (1.0 - bred) ** q # bred must be < 1 to avoid NaNs + + jcrit = jc1 * jc2 * jc3 / scalefac + + return jcrit, bcrit, tcrit + + +def jcrit_nbti(temperature, bmax, c0, bc20max, tc0max): + """Critical current density in a NbTi superconductor strand + author: P J Knight, CCFE, Culham Science Centre + temperature : input real : SC temperature (K) + bmax : input real : Magnetic field at conductor (T) + c0 : input real : Scaling constant (A/m2) + bc20max : input real : Upper critical field (T) for superconductor + at zero temperature and strain + tc0max : input real : Critical temperature (K) at zero field and strain + jcrit : output real : Critical current density in superconductor (A/m2) + tcrit : output real : Critical temperature (K) + This routine calculates the critical current density and + temperature in superconducting TF coils using NbTi + as the superconductor. + """ + + bratio = bmax / bc20max + + if bmax < bc20max: + # Critical temperature (K) + tcrit = tc0max * (1.0 - bratio) ** 0.59 + else: + # Allow bmax > bc20max but set error flag + # Fudge to give real (negative) value if bratio < 1 + tcrit = tc0max * (1.0 - bratio) + + # Allow tbar to be negative but set error flag + tbar = 1.0 - temperature / tcrit + + # Critical current density (A/m2) + jcrit = c0 * (1.0 - bratio) * tbar + + # if ((temperature > tcrit).or.(bmax > bc20max))then + # write(*,*)'jcrit_nbti: out of range: ', 'bmax =', bmax, ' bc20max =', bc20max, & + # ' temperature =',temperature, ' tcrit =',tcrit + # end if + + return jcrit, tcrit + + +def bi2212(bmax, jstrand, tsc, fhts): + """Fitted parameterization to Bi-2212 superconductor properties + author: P J Knight, CCFE, Culham Science Centre + author: M Kovari, CCFE, Culham Science Centre + bmax : input real : Magnetic field at conductor (T) + jstrand : input real : Current density in strand (A/m2) + tsc : input real : Superconductor temperature (K) + fhts : input real : Adjustment factor (<= 1) to account for strain, + radiation damage, fatigue or AC losses + jcrit : output real : Critical current density in strand (A/m2) + tmarg : output real : Temperature margin (K) + This routine calculates the critical current density and + the temperature margin for Bi-2212 superconductor in the TF coils + using a fit by M. Kovari to measurements described in the reference, + specifically from the points shown in Figure 6. +

Bi-2212 (Bi2Sr2CaCu2O8-x) + is a first-generation high temperature superconductor; it still needs + to be operated below about 10K, but remains superconducting at much + higher fields at that temperature than Nb3Sn etc. + The model's range of validity is T < 20K, adjusted field + b < 104 T, B > 6 T. + A transformative superconducting magnet technology for fields well + above 30 T using isotropic round wire multifilament + Bi2Sr2CaCu2O8-x conductor, D. C. Larbalestier et al., preprint, + 9th April 2013 + """ + b = bmax / np.exp(-0.168 * (tsc - 4.2)) + + # Engineering (i.e. strand) critical current density (A/m2) + + jcrit = fhts * (1.175e9 * np.exp(-0.02115 * b) - 1.288e8) + + # Temperature margin (K) + # Simple inversion of above calculation, using actual current density + # in strand instead of jcrit + + tmarg = ( + 1.0 + / 0.168 + * np.log(np.log(1.175e9 / (jstrand / fhts + 1.288e8)) / (0.02115 * bmax)) + + 4.2 + - tsc + ) + + # Check if ranges of validity have been violated + + if (tsc > 20.0) or (bmax < 6.0) or (b > 104.0): + eh.fdiags[0] = tsc + eh.fdiags[1] = bmax + eh.fdiags[2] = b + eh.report_error(106) + + return jcrit, tmarg + + +def gl_nbti(thelium, bmax, strain, bc20max, t_c0): + """Author: S B L Chislett-McDonald Durham University + Category: subroutine + + Critical current density of the superconductor in an ITER + Nb-Ti strand based on the Ginzburg-Landau theory of superconductivity + + \\begin{equation} + J_{c,TS}(B,T,\\epsilon_{I}) = A(\\epsilon_{I}) \\left[T_{c}(\\epsilon_{I})*(1-t^2)\\right]^2\\left + [B_{c2}(\\epsilon_I)*(1-t^\\nu)\\right]^{n-3}b^{p-1}(1-b)^q~. + \\end{equation} + + - \\( \\thelium \\) -- Coolant/SC temperature [K] + - \\( \\bmax \\) -- Magnetic field at conductor [T] + - \\( \\epsilon_{I} \\) -- Intrinsic strain in superconductor [\\%] + - \\( \\B_{c2}(\\epsilon_I) \\) -- Strain dependent upper critical field [T] + - \\( \\b \\) -- Reduced field = bmax / \\B_{c2}(\\epsilon_I)*(1-t^\\nu) [unitless] + - \\( \\T_{c}(\\epsilon_{I}) \\) -- Strain dependent critical temperature (K) + - \\( \\t \\) -- Reduced temperature = thelium / \\T_{c}(\\epsilon_{I}) [unitless] + - \\( \\A(\\epsilon_{I}) \\) -- Strain dependent Prefactor [A / ( m\\(^2\\) K\\(^-2) T\\(^n-3))] + - \\( \\J_{c,TS} \\) -- Critical current density in superconductor [A / m\\(^-2\\)] + """ + + A_0 = 1102e6 + p = 0.49 + q = 0.56 + n = 1.83 + v = 1.42 + c2 = -0.0025 + c3 = -0.0003 + c4 = -0.0001 + em = -0.002e-2 + u = 0.0 + w = 2.2 + + epsilon_I = strain - em + + strain_func = ( + 1 + c2 * (epsilon_I) ** 2 + c3 * (epsilon_I) ** 3 + c4 * (epsilon_I) ** 4 + ) + + T_e = t_c0 * strain_func ** (1 / w) + + t_reduced = thelium / T_e + + A_e = A_0 * strain_func ** (u / w) + + # Critical Field + bcrit = bc20max * (1 - t_reduced**v) * strain_func + + b_reduced = bmax / bcrit + + # Critical temperature (K) + tcrit = T_e + + # Critical current density (A/m2) + if b_reduced <= 1.0: + jcrit = ( + A_e + * (T_e * (1 - t_reduced**2)) ** 2 + * bcrit ** (n - 3) + * b_reduced ** (p - 1) + * (1 - b_reduced) ** q + ) + else: # Fudge to yield negative single valued function of Jc for fields above Bc2 + jcrit = ( + A_e + * (T_e * (1 - t_reduced**2)) ** 2 + * bcrit ** (n - 3) + * b_reduced ** (p - 1) + * (1 - b_reduced) ** 1.0 + ) + + return jcrit, bcrit, tcrit + + +def gl_rebco(thelium, bmax, strain, bc20max, t_c0): + """Author: S B L Chislett-McDonald Durham University + Category: subroutine + + Critical current density of a SuperPower REBCO tape based on measurements by P. Branch + at Durham University + https://git.ccfe.ac.uk/process/process/uploads/e98c6ea13da782cdc6fe16daea92078a/20200707_Branch-Osamura-Hampshire_-_accepted_SuST.pdf + and fit to state-of-the-art measurements at 4.2 K published in SuST + http://dx.doi.org/10.1088/0953-2048/24/3/035001 + + \\begin{equation} + J_{c,TS}(B,T,\\epsilon_{I}) = A(\\epsilon_{I}) \\left[T_{c}(\\epsilon_{I})*(1-t^2)\\right]^2\\left + [B_{c2}(\\epsilon_I)*(1-t)^s\\right]^{n-3}b^{p-1}(1-b)^q~. + \\end{equation} + + - \\( \\thelium \\) -- Coolant/SC temperature [K] + - \\( \\bmax \\) -- Magnetic field at conductor [T] + - \\( \\epsilon_{I} \\) -- Intrinsic strain in superconductor [\\%] + - \\( \\B_{c2}(\\epsilon_I) \\) -- Strain dependent upper critical field [T] + - \\( \\b \\) -- Reduced field = bmax / \\B_{c2}(\\epsilon_I)*(1-t^\\nu) [unitless] + - \\( \\T_{c}(\\epsilon_{I}) \\) -- Strain dependent critical temperature (K) + - \\( \\t \\) -- Reduced temperature = thelium / \\T_{c}(\\epsilon_{I}) [unitless] + - \\( \\A(epsilon_{I}) \\) -- Strain dependent Prefactor [A / ( m\\(^2\\) K\\(^-2) T\\(^n-3))] + - \\( \\J_{c,TS} \\) -- Critical current density in superconductor [A / m\\(^-2\\)] + - \\( \\epsilon_{m} \\) -- Strain at which peak in J_c occurs [\\%] + """ + # critical current density prefactor + A_0 = 2.95e2 + # flux pinning field scaling parameters + p = 0.32 + q = 2.50 + n = 3.33 + # temperatute scaling parameter + s = 5.27 + # strain scaling parameters + c2 = -0.0191 + c3 = 0.0039 + c4 = 0.00103 + em = 0.058 + # strain conversion parameters + u = 0.0 + w = 2.2 + + epsilon_I = strain - em + + strain_func = ( + 1 + c2 * (epsilon_I) ** 2 + c3 * (epsilon_I) ** 3 + c4 * (epsilon_I) ** 4 + ) + + T_e = t_c0 * strain_func ** (1 / w) + + t_reduced = thelium / T_e + + A_e = A_0 * strain_func ** (u / w) + + # Critical Field + bcrit = bc20max * (1 - t_reduced) ** s * strain_func + + b_reduced = bmax / bcrit + + # Critical temperature (K) + tcrit = T_e + + # Critical current density (A/m2) + jcrit = ( + A_e + * (T_e * (1 - t_reduced**2)) ** 2 + * bcrit ** (n - 3) + * b_reduced ** (p - 1) + * (1 - b_reduced) ** q + ) + + return jcrit, bcrit, tcrit + + +def hijc_rebco(thelium, bmax, strain, bc20max, t_c0): + """Implementation of High Current Density REBCO tape + author: R Chapman, UKAEA + thelium : input real : SC temperature (K) + bmax : input real : Magnetic field at conductor (T) + strain : input real : Strain in superconductor + bc20max : input real : Upper critical field (T) for superconductor + at zero temperature and strain + t_c0 : input real : Critical temperature (K) at zero field and strain + jcrit : output real : Critical current density in superconductor (A/m2) + bcrit : output real : Critical field (T) + tcrit : output real : Critical temperature (K) + + Returns the critical current of a REBCO tape based on a critical surface + (field, temperature) parameterization. Based in part on the parameterization + described in: M. J. Wolf, N. Bagrets, W. H. Fietz, C. Lange and K. Weiss, + "Critical Current Densities of 482 A/mm2 in HTS CrossConductors at 4.2 K and 12 T," + in IEEE Transactions on Applied Superconductivity, vol. 28, no. 4, pp. 1-4, + June 2018, Art no. 4802404, doi: 10.1109/TASC.2018.2815767. And on the experimental + data presented here: "2G HTS Wire Development at SuperPower", Drew W. Hazelton, + February 16, 2017 https://indico.cern.ch/event/588810/contributions/2473740/ + The high Ic parameterization is a result of modifications based on Ic values + observed in: "Conceptual design of HTS magnets for fusion nuclear science facility", + Yuhu Zhai, Danko van der Laan, Patrick Connolly, Charles Kessel, 2021, + https://doi.org/10.1016/j.fusengdes.2021.112611 + The parameter A is transformed into a function A(T) based on a Newton polynomial fit + considering A(4.2 K) = 2.2e8, A(20 K) = 2.3e8 and A(65 K) = 3.5e8. These values were + selected manually. A good fit to the pubished data can be seen in the 4-10 T range + but the fit deviates at very low or very high field. + """ + + a = 1.4 + b = 2.005 + # critical current density prefactor + A_0 = 2.2e8 + # flux pinning field scaling parameters + p = 0.39 + q = 0.9 + # strain conversion parameters + u = 33450.0 + v = -176577.0 + + # Critical Field (T) + # B_crit(T) calculated using temperature and critical temperature + bcrit = bc20max * (1.0 - thelium / t_c0) ** a + + # Critical temperature (K) + # scaled to match behaviour in GL_REBCO routine, + # ONLY TO BE USED until a better suggestion is received + tcrit = 0.999965 * t_c0 + + # finding A(T); constants based on a Newton polynomial fit to pubished data + A_t = A_0 + (u * thelium**2) + (v * thelium) + + # Critical current density (A/m2) + + jcrit = (A_t / bmax) * bcrit**b * (bmax / bcrit) ** p * (1 - bmax / bcrit) ** q + + # Jc times HTS area: default area is width 4mm times HTS layer thickness 1 um, + # divided by the tape area to provide engineering Jc per tape, then multiplied by fraction 0.4 + # to reach the level of current density expected in the space where the tapes are wound in A/m^2! + jcrit = ( + jcrit + * (rebco_variables.tape_width * rebco_variables.rebco_thickness) + / (rebco_variables.tape_width * rebco_variables.tape_thickness) + * 0.4 + ) + + return jcrit, bcrit, tcrit + + +def wstsc(temperature, bmax, strain, bc20max, tc0max): + """Implementation of WST Nb3Sn critical surface implementation + author: J Morris, CCFE, Culham Science Centre + temperature : input real : SC temperature (K) + bmax : input real : Magnetic field at conductor (T) + strain : input real : Strain in superconductor + bc20max : input real : Upper critical field (T) for superconductor + at zero temperature and strain + tc0max : input real : Critical temperature (K) at zero field and strain + jcrit : output real : Critical current density in superconductor (A/m2) + bcrit : output real : Critical field (T) + tcrit : output real : Critical temperature (K) + This routine calculates the critical current density and + temperature in the superconducting TF coils using the + WST Nb3Sn critical surface model. + V. Corato et al, "Common operating values for DEMO magnets design for 2016", + https://scipub.euro-fusion.org/wp-content/uploads/eurofusion/WPMAGREP16_16565_submitted.pdf + """ + + # Scaling constant C [AT/mm2] + csc = 83075.0 + # Low field exponent p + p = 0.593 + # High field exponent q + q = 2.156 + # Strain fitting constant C_{a1} + ca1 = 50.06 + # Strain fitting constant C_{a2} + ca2 = 0.0 + # epsilon_{0,a} + eps0a = 0.00312 + + # $\epsilon_{sh}$ + epssh = (ca2 * eps0a) / (np.sqrt(ca1**2 - ca2**2)) + + # Strain function $s(\epsilon)$ + # 0.83 < s < 1.0, for -0.005 < strain < 0.005 + strfun = np.sqrt(epssh**2 + eps0a**2) - np.sqrt( + (strain - epssh) ** 2 + eps0a**2 + ) + strfun = strfun * ca1 - ca2 * strain + strfun = 1.0 + (1.0 / (1.0 - ca1 * eps0a)) * strfun + if strfun < 0: + logger.warning(f"Strain function < 0: {strfun} for {strain = }") + + # $B^*_{C2} (0,\epsilon)$ + bc20eps = bc20max * strfun + + # $T^*_C (0,\epsilon)$ + tc0eps = tc0max * strfun ** (1.0 / 3.0) + + # Reduced temperature + # Should remain < 1 for temperature < 0.94*tc0max (i.e. 15 kelvin for i_tf_sc_mat=1) + + if temperature / tc0eps >= 1.0: + eh.fdiags[0] = temperature + eh.fdiags[1] = tc0eps + eh.report_error(159) + + # t = min(thelium/tc0eps, 0.9999D0) + t = temperature / tc0eps + + # Reduced magnetic field at zero temperature + # Should remain < 1 for bmax < 0.83*bc20max (i.e. 27 tesla for i_tf_sc_mat=1) + + if bmax / bc20eps >= 1.0: + eh.fdiags[0] = bmax + eh.fdiags[1] = bc20eps + eh.report_error(160) + + # bzero = min(bmax/bc20eps, 0.9999D0) + bzero = bmax / bc20eps + + if bzero < 1.0: + # Critical temperature (K) + tcrit = tc0eps * (1.0 - bzero) ** (1.0 / 1.52) + else: + # Allow bzero > 1, fudge to give real (negative) value of tcrit + # This generates a real (negative) and continuous (but not differentiable) + # function of bzero. + tcrit = tc0eps + + # Critical field (T). Negative if normalised temperature t>1 + if t > 0.0: + bcrit = bc20eps * (1.0 - t**1.52) + else: + # Allow t<0, fudge to give real value of bcrit + bcrit = bc20eps * (1.0 - t) + + # Reduced magnetic field, restricted to be < 1 + if bmax / bcrit >= 1.0: + eh.fdiags[0] = bmax + eh.fdiags[1] = bcrit + eh.report_error(161) + + # bred = min(bmax/bcrit, 0.9999D0) + bred = bmax / bcrit + + if (bred > 0.0) and (bred < 1.0): + jc3 = bred**p * (1.0 - bred) ** q # bred must be < 1 to avoid NaNs + else: + # Allow bred > 1 or <0, fudge to give real (negative) value of jc3 + # This generates a real (negative) and continuous (but not differentiable) + # function of bred. + jc3 = bred * (1.0 - bred) + if not np.isfinite(jc3): + raise RuntimeError("jc3 jcrit is NaN.") + + # Critical current density in superconductor (A/m2) + jc1 = (csc / bmax) * strfun + + if t > 0.0: + jc2 = (1.0 - t**1.52) * (1.0 - t**2) + else: + # Allow t<0, fudge to give real value of jc2 + # This generates a real and continuous (but not differentiable) function of t. + jc2 = (1.0 - t) * (1.0 - t**2) + + # jc3 = bred**p * (1.0D0-bred)**q ! bred must be < 1 to avoid NaNs + + # scale from mm2 to m2 + scalefac = 1.0e6 + + jcrit = jc1 * jc2 * jc3 * scalefac + if not np.isfinite(jcrit): + raise RuntimeError("jcrit is NaN.") + + return jcrit, bcrit, tcrit + + +def croco(jcritsc, conductor_area, croco_od, croco_thick): + """'CroCo' (cross-conductor) strand and cable design for + 'REBCO' 2nd generation HTS superconductor + Updated 13/11/18 using data from Lewandowska et al 2018. + """ + d = croco_od + # d = conductor_width / 3.0d0 - thwcndut * ( 2.0d0 / 3.0d0 ) + + croco_id = d - 2.0 * croco_thick # scaling * 5.4d-3 + if croco_id <= 0.0: + logger.warning("Negitive inner croco diameter") + + # Define the scaling factor for the input REBCO variable + # Ratio of new croco inner diameter and fixed base line value + scaling = croco_id / 5.4e-3 + tape_width = scaling * 3.75e-3 + # Properties of a single strand + tape_thickness = ( + rebco_variables.rebco_thickness + + rebco_variables.copper_thick + + rebco_variables.hastelloy_thickness + ) + stack_thickness = np.sqrt(croco_id**2 - tape_width**2) + tapes = stack_thickness / tape_thickness + + copper_area = ( + np.pi * croco_thick * d + - np.pi * croco_thick**2 + + rebco_variables.copper_thick + * tape_width + * tapes # copper tube # copper in tape + ) + hastelloy_area = rebco_variables.hastelloy_thickness * tape_width * tapes + solder_area = np.pi / 4.0 * croco_id**2 - stack_thickness * tape_width + + rebco_area = rebco_variables.rebco_thickness * tape_width * tapes + croco_strand_area = np.pi / 4.0 * d**2 + croco_strand_critical_current = jcritsc * rebco_area + + # Conductor properties + # conductor%number_croco = conductor%acs*(1.0-cable_helium_fraction-copper_bar)/croco_strand_area + conductor_critical_current = croco_strand_critical_current * 6.0 + # Area of core = area of strand + conductor_copper_bar_area = croco_strand_area + conductor_copper_area = copper_area * 6.0 + conductor_copper_bar_area + conductor_copper_fraction = conductor_copper_area / conductor_area + + # Helium area is set by the user. + # conductor_helium_area = cable_helium_fraction * conductor_acs + conductor_helium_area = np.pi / 2.0 * d**2 + conductor_helium_fraction = conductor_helium_area / conductor_area + + conductor_hastelloy_area = hastelloy_area * 6.0 + conductor_hastelloy_fraction = conductor_hastelloy_area / conductor_area + + conductor_solder_area = solder_area * 6.0 + conductor_solder_fraction = conductor_solder_area / conductor_area + + conductor_rebco_area = rebco_area * 6.0 + conductor_rebco_fraction = conductor_rebco_area / conductor_area + + return ( + croco_strand_area, + croco_strand_critical_current, + conductor_copper_area, + conductor_copper_fraction, + conductor_copper_bar_area, + conductor_hastelloy_area, + conductor_hastelloy_fraction, + conductor_helium_area, + conductor_helium_fraction, + conductor_solder_area, + conductor_solder_fraction, + conductor_rebco_area, + conductor_rebco_fraction, + conductor_critical_current, + ) diff --git a/source/fortran/main_module.f90 b/source/fortran/main_module.f90 index 6f40f37611..23abe16672 100644 --- a/source/fortran/main_module.f90 +++ b/source/fortran/main_module.f90 @@ -684,7 +684,6 @@ subroutine runtests use maths_library, only: nearly_equal, binomial, test_secant_solve use process_output, only: ocmmnt, ovarre ! use pfcoil_module, only: brookscoil - use superconductors, only: test_quench ! use reinke_module, only: test_reinke implicit none real(dp) :: fshift, xf, enpa,ftherm,fpp,cdeff, ampperwatt @@ -696,7 +695,6 @@ subroutine runtests call ovarre(nout,'Binomial coefficients C(5,4): 5', '(binomial(5,4))', binomial(5,4)) call ovarre(nout,'Binomial coefficients C(5,5): 1', '(binomial(5,5))', binomial(5,5)) - call test_quench() ! call brookscoil(nout) Moved to pytest call test_secant_solve() ! Disabled for ease of #1542 - Tim diff --git a/source/fortran/superconductors.f90 b/source/fortran/superconductors.f90 deleted file mode 100755 index 73c76d10df..0000000000 --- a/source/fortran/superconductors.f90 +++ /dev/null @@ -1,1905 +0,0 @@ -module superconductors - !! Module containing superconducter critical surfaces and conductor data -#ifndef dp - use, intrinsic :: iso_fortran_env, only: dp=>real64 -#endif - implicit none -contains - -! ------------------------------------------------------------------- - -subroutine jcrit_rebco(temperature, b, jcrit, validity, iprint) - - !! Critical current density for "REBCO" 2nd generation HTS superconductor - !! temperature : input real : superconductor temperature (K) - !! b : input real : Magnetic field at superconductor (T) - !! jcrit : output real : Critical current density in superconductor (A/m2) - ! Will return a negative number if the temperature is greater than Tc0, the - ! zero-field critical temperature. - use constants, only: pi - implicit none - - ! Arguments - real(dp), intent(in) :: temperature, b - real(dp), intent(out) :: jcrit - logical, intent(out) :: validity - integer, intent(in) :: iprint - - ! Local variables - real(dp) :: birr, factor, tcb - - ! Parameters - real(dp), parameter :: tc0 = 90.0d0 ! (K) - real(dp), parameter :: birr0 = 132.5d0 ! (T) - real(dp), parameter :: a = 1.82962d8 ! scaling constant - real(dp), parameter :: p = 0.5875d0 ! exponent - real(dp), parameter :: q = 1.7d0 ! exponent - real(dp), parameter :: alpha =1.54121d0 ! exponent - real(dp), parameter :: beta = 1.96679d0 ! exponent - real(dp), parameter :: oneoveralpha = 1d0 / alpha ! inverse - - validity = .true. - if((temperature<4.2d0).or.(temperature>72.0d0) )validity = .false. - if(temperature<65)then - if((b<0.0d0).or.(b>15.0d0))validity = .false. - else - if((b<0.0d0).or.(b>11.5d0))validity = .false. - endif - - if((iprint==1).and.(.not.validity))then - write(*,*)'ERROR in subroutine jcrit_rebco: input out of range' - write(*,*)'temperature = ', temperature, ' Field = ', b - endif - - if(temperature= 1.0D0) then - fdiags(1) = thelium ; fdiags(2) = tc0eps - call report_error(159) - end if - t = min(thelium/tc0eps, 0.9999D0) - - ! Reduced magnetic field at zero temperature - ! Should remain < 1 for bmax < 0.83*bc20max (i.e. 27 tesla for i_tf_sc_mat=1) - - if (bmax/bc20eps >= 1.0D0) then - fdiags(1) = bmax ; fdiags(2) = bc20eps - call report_error(160) - end if - bzero = min(bmax/bc20eps, 0.9999D0) - - ! Critical temperature (K) - tcrit = tc0eps * (1.0D0 - bzero)**(1.0D0/1.52D0) ! bzero must be < 1 to avoid NaNs - - ! Critical field (T) - bcrit = bc20eps * (1.0D0 - t**1.52D0) - - ! Reduced magnetic field, restricted to be < 1 - if (bmax/bcrit >= 1.0D0) then - fdiags(1) = bmax ; fdiags(2) = bcrit - call report_error(161) - end if - bred = min(bmax/bcrit, 0.9999D0) - - ! Critical current density in superconductor (A/m2) - ! ITER parameterization is for the current in a single strand, - ! not per unit area, so scalefac converts to current density - - scalefac = pi * (0.5D0*diter)**2 * (1.0D0-cuiter) - - jc1 = (csc/bmax)*strfun - jc2 = (1.0D0-t**1.52D0) * (1.0D0-t**2) ! t must be < 1 to avoid NaNs - jc3 = bred**p * (1.0D0-bred)**q ! bred must be < 1 to avoid NaNs - - jcrit = jc1 * jc2 * jc3 / scalefac - -end subroutine itersc - -! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -subroutine bi2212(bmax,jstrand,tsc,fhts,jcrit,tmarg) - - !! Fitted parameterization to Bi-2212 superconductor properties - !! author: P J Knight, CCFE, Culham Science Centre - !! author: M Kovari, CCFE, Culham Science Centre - !! bmax : input real : Magnetic field at conductor (T) - !! jstrand : input real : Current density in strand (A/m2) - !! tsc : input real : Superconductor temperature (K) - !! fhts : input real : Adjustment factor (<= 1) to account for strain, - !! radiation damage, fatigue or AC losses - !! jcrit : output real : Critical current density in strand (A/m2) - !! tmarg : output real : Temperature margin (K) - !! This routine calculates the critical current density and - !! the temperature margin for Bi-2212 superconductor in the TF coils - !! using a fit by M. Kovari to measurements described in the reference, - !! specifically from the points shown in Figure 6. - !!

Bi-2212 (Bi2Sr2CaCu2O8-x) - !! is a first-generation high temperature superconductor; it still needs - !! to be operated below about 10K, but remains superconducting at much - !! higher fields at that temperature than Nb3Sn etc. - !! The model's range of validity is T < 20K, adjusted field - !! b < 104 T, B > 6 T. - !! A transformative superconducting magnet technology for fields well - !! above 30 T using isotropic round wire multifilament - !! Bi2Sr2CaCu2O8-x conductor, D. C. Larbalestier et al., preprint, - !! 9th April 2013 - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use error_handling, only: fdiags, report_error - implicit none - - ! Arguments - - real(dp), intent(in) :: bmax, jstrand, tsc, fhts - real(dp), intent(out) :: jcrit, tmarg - - ! Local variables - - real(dp) :: b - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Adjusted field (T) - - b = bmax / exp(-0.168D0*(tsc-4.2D0)) - - ! Engineering (i.e. strand) critical current density (A/m2) - - jcrit = fhts * (1.175D9*exp(-0.02115D0*b) - 1.288D8) - - ! Temperature margin (K) - ! Simple inversion of above calculation, using actual current density - ! in strand instead of jcrit - - tmarg = 1.0D0/0.168D0 * & - log( log(1.175D9/(jstrand/fhts + 1.288D8)) / (0.02115*bmax) ) & - + 4.2D0 - tsc - - ! Check if ranges of validity have been violated - - if ((tsc > 20.0D0).or.(bmax < 6.0D0).or.(b > 104.0D0)) then - fdiags(1) = tsc ; fdiags(2) = bmax ; fdiags(3) = b - call report_error(106) - write(*,*) 'Warning in routine BI2212:' - write(*,*) 'Range of validity of the HTS Bi-2212 model has been violated:' - write(*,*) ' S/C temperature (K) = ',tsc, ' (should be < 20 K)' - write(*,*) 'Field at conductor (T) = ',bmax, ' (should be > 6 T)' - write(*,*) ' Adjusted field (T) = ',b, ' (should be < 104 T)' - write(*,*) ' ' - end if - -end subroutine bi2212 -!------------------------------------------------------------------ -subroutine jcrit_nbti(temperature,bmax,c0,bc20max,tc0max,jcrit,tcrit) - - !! Critical current density in a NbTi superconductor strand - !! author: P J Knight, CCFE, Culham Science Centre - !! temperature : input real : SC temperature (K) - !! bmax : input real : Magnetic field at conductor (T) - !! c0 : input real : Scaling constant (A/m2) - !! bc20max : input real : Upper critical field (T) for superconductor - !! at zero temperature and strain - !! tc0max : input real : Critical temperature (K) at zero field and strain - !! jcrit : output real : Critical current density in superconductor (A/m2) - !! tcrit : output real : Critical temperature (K) - !! This routine calculates the critical current density and - !! temperature in superconducting TF coils using NbTi - !! as the superconductor. - !! None - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use constants, only: pi - implicit none - - ! Arguments - real(dp), intent(in) :: temperature, bmax, c0, bc20max, tc0max - real(dp), intent(out) :: jcrit, tcrit - - ! Local variables - real(dp) :: bratio, tbar - !----------------------------------- - bratio = bmax/bc20max - - if (bmax < bc20max) then - ! Critical temperature (K) - tcrit = tc0max * (1.0D0 - bratio)**0.59D0 - else - ! Allow bmax > bc20max but set error flag - ! Fudge to give real (negative) value if bratio < 1 - tcrit = tc0max * (1.0D0 - bratio) - end if - - ! Allow tbar to be negative but set error flag - tbar = 1.0D0 - temperature/tcrit - - ! Critical current density (A/m2) - jcrit = c0 * (1.0D0 - bratio) * tbar - - ! if ((temperature > tcrit).or.(bmax > bc20max))then - ! write(*,*)'jcrit_nbti: out of range: ', 'bmax =', bmax, ' bc20max =', bc20max, & - ! ' temperature =',temperature, ' tcrit =',tcrit - ! end if - -end subroutine jcrit_nbti -!-------------------------------------------------------------------- -subroutine GL_nbti(thelium,bmax,strain,bc20max,t_c0,jcrit,bcrit,tcrit) - - !! Author: S B L Chislett-McDonald Durham University - !! Category: subroutine - !! - !! Critical current density of the superconductor in an ITER - !! Nb-Ti strand based on the Ginzburg-Landau theory of superconductivity - !! - !! \begin{equation} - !! J_{c,TS}(B,T,\epsilon_{I}) = A(\epsilon_{I}) \left[T_{c}(\epsilon_{I})*(1-t^2)\right]^2\left - !! [B_{c2}(\epsilon_I)*(1-t^\nu)\right]^{n-3}b^{p-1}(1-b)^q~. - !! \end{equation} - !! - !! - \( \thelium \) -- Coolant/SC temperature [K] - !! - \( \bmax \) -- Magnetic field at conductor [T] - !! - \( \\epsilon_{I} \) -- Intrinsic strain in superconductor [\%] - !! - \( \B_{c2}(\epsilon_I) \) -- Strain dependent upper critical field [T] - !! - \( \b \) -- Reduced field = bmax / \B_{c2}(\epsilon_I)*(1-t^\nu) [unitless] - !! - \( \T_{c}(\epsilon_{I}) \) -- Strain dependent critical temperature (K) - !! - \( \t \) -- Reduced temperature = thelium / \T_{c}(\epsilon_{I}) [unitless] - !! - \( \A(epsilon_{I}) \) -- Strain dependent Prefactor [A / ( m\(^2\) K\(^-2) T\(^n-3))] - !! - \( \J_{c,TS} \) -- Critical current density in superconductor [A / m\(^-2\)] - - implicit none - - ! Arguments - real(dp), intent(in) :: thelium, bmax, strain, bc20max, t_c0 - real(dp), intent(out) :: jcrit, tcrit, bcrit - - ! Local variables - real(dp) :: strain_func, T_e, A_e, b_reduced, t_reduced, epsilon_I - !critical current density prefactor (strand non-copper J_c) - real(dp), parameter :: A_0 = 1102D6 - !flux pinning field scaling parameters - real(dp), parameter :: p = 0.49D0 - real(dp), parameter :: q = 0.56D0 - real(dp), parameter :: n = 1.83D0 - !temperatute scaling parameter - real(dp), parameter :: v = 1.42D0 - !strain scaling parameters - real(dp), parameter :: c2 = -0.0025D0 - real(dp), parameter :: c3 = -0.0003D0 - real(dp), parameter :: c4 = -0.0001D0 - real(dp), parameter :: em = -0.002D-2 - !strain conversion parameters - real(dp), parameter :: u = 0.0D0 - real(dp), parameter :: w = 2.2D0 - - epsilon_I = strain - em - - strain_func = 1 + c2*(epsilon_I)**2 + c3*(epsilon_I)**3 + c4*(epsilon_I)**4 - - T_e = t_c0 * strain_func**(1 / w) - - t_reduced = thelium/T_e - - A_e = A_0 * strain_func**(u / w) - - ! Critical Field - bcrit = bc20max * (1 - t_reduced**v) * strain_func - - b_reduced = bmax/bcrit - - ! Critical temperature (K) - tcrit = T_e - - ! Critical current density (A/m2) - if (b_reduced <= 1.0D0) then - jcrit = A_e * (T_e*(1-t_reduced**2))**2 * bcrit**(n-3) * b_reduced**(p-1) * (1 - b_reduced)**q - else !Fudge to yield negative single valued function of Jc for fields above Bc2 - jcrit = A_e * (T_e*(1-t_reduced**2))**2 * bcrit**(n-3) * b_reduced**(p-1) * (1 - b_reduced)**1.0 - end if - -end subroutine GL_nbti - -subroutine wstsc(temperature,bmax,strain,bc20max,tc0max,jcrit,bcrit,tcrit) - - !! Implementation of WST Nb3Sn critical surface implementation - !! author: J Morris, CCFE, Culham Science Centre - !! temperature : input real : SC temperature (K) - !! bmax : input real : Magnetic field at conductor (T) - !! strain : input real : Strain in superconductor - !! bc20max : input real : Upper critical field (T) for superconductor - !! at zero temperature and strain - !! tc0max : input real : Critical temperature (K) at zero field and strain - !! jcrit : output real : Critical current density in superconductor (A/m2) - !! bcrit : output real : Critical field (T) - !! tcrit : output real : Critical temperature (K) - !! This routine calculates the critical current density and - !! temperature in the superconducting TF coils using the - !! WST Nb3Sn critical surface model. - !! V. Corato et al, "Common operating values for DEMO magnets design for 2016", - !! https://scipub.euro-fusion.org/wp-content/uploads/eurofusion/WPMAGREP16_16565_submitted.pdf - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use error_handling, only: fdiags, report_error - use maths_library, only: variable_error - implicit none - - ! Arguments - real(dp), intent(in) :: temperature, bmax, strain, bc20max, tc0max - real(dp), intent(out) :: jcrit, bcrit, tcrit - - ! Local variables - - ! Scaling constant C [AT/mm2] - real(dp), parameter :: csc = 83075.0D0 - ! Low field exponent p - real(dp), parameter :: p = 0.593D0 - ! High field exponent q - real(dp), parameter :: q = 2.156D0 - ! Strain fitting constant C_{a1} - real(dp), parameter :: ca1 = 50.06D0 - ! Strain fitting constant C_{a2} - real(dp), parameter :: ca2 = 0.0D0 - ! epsilon_{0,a} - real(dp), parameter :: eps0a = 0.00312D0 - - !real(dp), parameter :: epsmax = -1.1D-3 ! epsilon_{max} (not used) - - real(dp) :: bred, epssh, t, bc20eps, & - tc0eps, bzero, strfun, jc1, jc2, jc3, scalefac - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! $\epsilon_{sh}$ - epssh = (ca2*eps0a)/(sqrt(ca1**2 - ca2**2)) - - ! Strain function $s(\epsilon)$ - ! 0.83 < s < 1.0, for -0.005 < strain < 0.005 - strfun = sqrt(epssh**2 + eps0a**2) - sqrt((strain-epssh)**2 + eps0a**2) - strfun = strfun*ca1 - ca2*strain - strfun = 1.0D0 + (1.0D0/(1.0D0 - ca1*eps0a))*strfun - if(strfun<0.d0)write(*,*)'subroutine wstsc: strfun<0.d0. strfun =', strfun, ', strain = ',strain - - ! $B^*_{C2} (0,\epsilon)$ - bc20eps = bc20max*strfun - - ! $T^*_C (0,\epsilon)$ - tc0eps = tc0max * strfun**(1.0D0/3.0D0) - - ! Reduced temperature - ! Should remain < 1 for temperature < 0.94*tc0max (i.e. 15 kelvin for i_tf_sc_mat=1) - - if (temperature/tc0eps >= 1.0D0) then - fdiags(1) = temperature ; fdiags(2) = tc0eps - call report_error(159) - end if - ! t = min(thelium/tc0eps, 0.9999D0) - t = temperature/tc0eps - - ! Reduced magnetic field at zero temperature - ! Should remain < 1 for bmax < 0.83*bc20max (i.e. 27 tesla for i_tf_sc_mat=1) - - if (bmax/bc20eps >= 1.0D0) then - fdiags(1) = bmax ; fdiags(2) = bc20eps - call report_error(160) - end if - - ! bzero = min(bmax/bc20eps, 0.9999D0) - bzero = bmax/bc20eps - - if (bzero < 1.0d0) then - ! Critical temperature (K) - tcrit = tc0eps * (1.0D0 - bzero)**(1.0D0/1.52D0) - else - ! Allow bzero > 1, fudge to give real (negative) value of tcrit - ! This generates a real (negative) and continuous (but not differentiable) - ! function of bzero. - tcrit = tc0eps - end if - - - - ! Critical field (T). Negative if normalised temperature t>1 - if(t>0.0d0)then - bcrit = bc20eps * (1.0D0 - t**1.52D0) - else - ! Allow t<0, fudge to give real value of bcrit - bcrit = bc20eps * (1.0D0 - t) - end if - - ! Reduced magnetic field, restricted to be < 1 - if (bmax/bcrit >= 1.0D0) then - fdiags(1) = bmax ; fdiags(2) = bcrit - call report_error(161) - end if - ! bred = min(bmax/bcrit, 0.9999D0) - bred = bmax/bcrit - - if ((bred>0.0d0).and.(bred < 1.0d0)) then - jc3 = bred**p * (1.0D0-bred)**q ! bred must be < 1 to avoid NaNs - else - ! Allow bred > 1 or <0, fudge to give real (negative) value of jc3 - ! This generates a real (negative) and continuous (but not differentiable) - ! function of bred. - jc3 = bred * (1.0D0-bred) - if(variable_error(jc3))then - write(*,'(a24, 8(a12,es12.3))')'jc3 jcrit is NaN.',' bred=',bred, ' bmax=',bmax, ' bcrit=',bcrit, ' t=',t - stop 1 - end if - end if - - ! Critical current density in superconductor (A/m2) - jc1 = (csc/bmax)*strfun - - if(t>0.0d0)then - jc2 = (1.0D0-t**1.52D0) * (1.0D0-t**2) - else - ! Allow t<0, fudge to give real value of jc2 - ! This generates a real and continuous (but not differentiable) function of t. - jc2 = (1.0D0-t) * (1.0D0-t**2) - end if - - ! jc3 = bred**p * (1.0D0-bred)**q ! bred must be < 1 to avoid NaNs - - ! scale from mm2 to m2 - scalefac = 1.0D6 - - jcrit = jc1 * jc2 * jc3*scalefac - if(variable_error(jcrit))then - write(*,'(a24, 8(a12,es12.3))')'WST jcrit is NaN.',' jc1=',jc1, ' jc2=',jc2, ' jc3=',jc3, ' t=',t - write(*,'(a24, 8(a12,es12.3))')'T=',T,' bmax=',bmax,' strain=',strain,' bc20max=',bc20max, & - ' tc0max=',tc0max,'jcrit=',jcrit,' bcrit=',bcrit,' tcrit=', tcrit - stop 1 - end if - -end subroutine wstsc -!-------------------------------------------------------------------------- - -subroutine GL_REBCO(thelium,bmax,strain,bc20max,t_c0,jcrit,bcrit,tcrit) !SCM added 13/06/18 - - !! Author: S B L Chislett-McDonald Durham University - !! Category: subroutine - !! - !! Critical current density of a SuperPower REBCO tape based on measurements by P. Branch - !! at Durham University - !! https://git.ccfe.ac.uk/process/process/uploads/e98c6ea13da782cdc6fe16daea92078a/20200707_Branch-Osamura-Hampshire_-_accepted_SuST.pdf - !! and fit to state-of-the-art measurements at 4.2 K published in SuST - !! http://dx.doi.org/10.1088/0953-2048/24/3/035001 - !! - !! \begin{equation} - !! J_{c,TS}(B,T,\epsilon_{I}) = A(\epsilon_{I}) \left[T_{c}(\epsilon_{I})*(1-t^2)\right]^2\left - !! [B_{c2}(\epsilon_I)*(1-t)^s\right]^{n-3}b^{p-1}(1-b)^q~. - !! \end{equation} - !! - !! - \( \thelium \) -- Coolant/SC temperature [K] - !! - \( \bmax \) -- Magnetic field at conductor [T] - !! - \( \\epsilon_{I} \) -- Intrinsic strain in superconductor [\%] - !! - \( \B_{c2}(\epsilon_I) \) -- Strain dependent upper critical field [T] - !! - \( \b \) -- Reduced field = bmax / \B_{c2}(\epsilon_I)*(1-t^\nu) [unitless] - !! - \( \T_{c}(\epsilon_{I}) \) -- Strain dependent critical temperature (K) - !! - \( \t \) -- Reduced temperature = thelium / \T_{c}(\epsilon_{I}) [unitless] - !! - \( \A(epsilon_{I}) \) -- Strain dependent Prefactor [A / ( m\(^2\) K\(^-2) T\(^n-3))] - !! - \( \J_{c,TS} \) -- Critical current density in superconductor [A / m\(^-2\)] - !! - \( \\epsilon_{m} \) -- Strain at which peak in J_c occurs [\%] - - - implicit none - - ! Arguments - real(dp), intent(in) :: thelium, bmax, strain, bc20max, t_c0 - real(dp), intent(out) :: jcrit, tcrit, bcrit - - ! Local variables - real(dp) :: strain_func, T_e, A_e, b_reduced, t_reduced, epsilon_I - !critical current density prefactor - real(dp), parameter :: A_0 = 2.95D2 - !flux pinning field scaling parameters - real(dp), parameter :: p = 0.32D0 - real(dp), parameter :: q = 2.50D0 - real(dp), parameter :: n = 3.33D0 - !temperatute scaling parameter - real(dp), parameter :: s = 5.27D0 - !strain scaling parameters - real(dp), parameter :: c2 = -0.0191D0 - real(dp), parameter :: c3 = 0.0039D0 - real(dp), parameter :: c4 = 0.00103D0 - real(dp), parameter :: em = 0.058D0 - !strain conversion parameters - real(dp), parameter :: u = 0.0D0 - real(dp), parameter :: w = 2.2D0 - - epsilon_I = strain - em - - strain_func = 1 + c2*(epsilon_I)**2 + c3*(epsilon_I)**3 + c4*(epsilon_I)**4 - - T_e = t_c0 * strain_func**(1 / w) - - t_reduced = thelium/T_e - - A_e = A_0 * strain_func**(u / w) - - ! Critical Field - bcrit = bc20max * (1 - t_reduced)**s * strain_func - - b_reduced = bmax/bcrit - - ! Critical temperature (K) - tcrit = T_e - - ! Critical current density (A/m2) - jcrit = A_e * (T_e*(1-t_reduced**2))**2 * bcrit**(n-3) * b_reduced**(p-1) * (1 - b_reduced)**q - -end subroutine GL_REBCO - -!-------------------------------------------------------------------------- - -subroutine HIJC_REBCO(thelium,bmax,strain,bc20max,t_c0,jcrit,bcrit,tcrit) - - !! Implementation of High Current Density REBCO tape - !! author: R Chapman, UKAEA - !! thelium : input real : SC temperature (K) - !! bmax : input real : Magnetic field at conductor (T) - !! strain : input real : Strain in superconductor - !! bc20max : input real : Upper critical field (T) for superconductor - !! at zero temperature and strain - !! t_c0 : input real : Critical temperature (K) at zero field and strain - !! jcrit : output real : Critical current density in superconductor (A/m2) - !! bcrit : output real : Critical field (T) - !! tcrit : output real : Critical temperature (K) - !! - !! Returns the critical current of a REBCO tape based on a critical surface - !! (field, temperature) parameterization. Based in part on the parameterization - !! described in: M. J. Wolf, N. Bagrets, W. H. Fietz, C. Lange and K. Weiss, - !! "Critical Current Densities of 482 A/mm2 in HTS CrossConductors at 4.2 K and 12 T," - !! in IEEE Transactions on Applied Superconductivity, vol. 28, no. 4, pp. 1-4, - !! June 2018, Art no. 4802404, doi: 10.1109/TASC.2018.2815767. And on the experimental - !! data presented here: "2G HTS Wire Development at SuperPower", Drew W. Hazelton, - !! February 16, 2017 https://indico.cern.ch/event/588810/contributions/2473740/ - !! The high Ic parameterization is a result of modifications based on Ic values - !! observed in: "Conceptual design of HTS magnets for fusion nuclear science facility", - !! Yuhu Zhai, Danko van der Laan, Patrick Connolly, Charles Kessel, 2021, - !! https://doi.org/10.1016/j.fusengdes.2021.112611 - !! The parameter A is transformed into a function A(T) based on a Newton polynomial fit - !! considering A(4.2 K) = 2.2e8, A(20 K) = 2.3e8 and A(65 K) = 3.5e8. These values were - !! selected manually. A good fit to the pubished data can be seen in the 4-10 T range - !! but the fit deviates at very low or very high field. - ! - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - use rebco_variables, only: tape_thickness, tape_width, rebco_thickness - - implicit none - - ! Arguments - real(dp), intent(in) :: thelium, bmax, strain, bc20max, t_c0 - real(dp), intent(out) :: jcrit, tcrit, bcrit - - ! Local variables - real(dp) :: A_t - real(dp), parameter :: a = 1.4D0 - real(dp), parameter :: b = 2.005D0 - !critical current density prefactor - real(dp), parameter :: A_0 = 2.2D8 - !flux pinning field scaling parameters - real(dp), parameter :: p = 0.39D0 - real(dp), parameter :: q = 0.9D0 - !strain conversion parameters - real(dp), parameter :: u = 33450.0D0 - real(dp), parameter :: v = -176577.0D0 - - ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - - ! Critical Field (T) - ! B_crit(T) calculated using temperature and critical temperature - bcrit = bc20max * (1.0D0 - thelium/t_c0)**a - - ! Critical temperature (K) - ! scaled to match behaviour in GL_REBCO routine, - ! ONLY TO BE USED until a better suggestion is received - tcrit = 0.999965D0 * t_c0 - - ! finding A(T); constants based on a Newton polynomial fit to pubished data - A_t = A_0 + ( u * thelium**2 ) + ( v * thelium ) - - ! Critical current density (A/m2) - - jcrit = ( A_t / bmax ) * bcrit**b * ( bmax / bcrit )**p * ( 1 - bmax/bcrit )**q - - ! Jc times HTS area: default area is width 4mm times HTS layer thickness 1 um, - ! divided by the tape area to provide engineering Jc per tape, then multiplied by fraction 0.4 - ! to reach the level of current density expected in the space where the tapes are wound in A/m^2! - jcrit = jcrit * (tape_width * rebco_thickness) / (tape_width * tape_thickness) * 0.4D0 - -end subroutine HIJC_REBCO - -!---------------------------------------------------------------- - -subroutine croco(jcritsc, croco_strand_area, croco_strand_critical_current, & - conductor_copper_area, conductor_copper_fraction, conductor_copper_bar_area, & - conductor_hastelloy_area, conductor_hastelloy_fraction, conductor_helium_area, & - conductor_helium_fraction, conductor_solder_area, conductor_solder_fraction, & - conductor_rebco_area, conductor_rebco_fraction, conductor_critical_current, & - conductor_area, croco_od,croco_thick) - - !! "CroCo" (cross-conductor) strand and cable design for - !! "REBCO" 2nd generation HTS superconductor - ! Updated 13/11/18 using data from Lewandowska et al 2018. - - use rebco_variables, only: copper_area, copper_thick, croco_id, & - hastelloy_area, hastelloy_thickness, rebco_area, solder_area, & - stack_thickness, tape_thickness, tape_width, tapes, rebco_thickness - use resistive_materials, only: volume_fractions, supercon_strand - use constants, only: pi - implicit none - real(dp), intent(in) ::jcritsc - real(dp) :: d, scaling, croco_od, croco_thick - - ! conductor - real(dp), intent(inout) :: conductor_copper_area, conductor_copper_fraction - real(dp), intent(inout) :: conductor_copper_bar_area - real(dp), intent(inout) :: conductor_hastelloy_area, conductor_hastelloy_fraction - real(dp), intent(inout) :: conductor_helium_area, conductor_helium_fraction - real(dp), intent(inout) :: conductor_solder_area, conductor_solder_fraction - real(dp), intent(inout) :: conductor_rebco_area, conductor_rebco_fraction - real(dp), intent(inout) :: conductor_critical_current - real(dp), intent(in) :: conductor_area - - ! croco_strand - real(dp), intent(inout) :: croco_strand_area - real(dp), intent(inout) :: croco_strand_critical_current - - - ! Define local alias - d = croco_od - !d = conductor_width / 3.0d0 - thwcndut * ( 2.0d0 / 3.0d0 ) - - croco_id = d - 2.0d0 * croco_thick !scaling * 5.4d-3 - if (croco_id <= 0.0d0) then - write(*,*) 'Warning: negitive inner croco diameter!' - write(*,*)'croco_id =', croco_id, ',croco_thick = ', croco_thick, ', croco_od =', croco_od - end if - ! Define the scaling factor for the input REBCO variable - ! Ratio of new croco inner diameter and fixed base line value - scaling = croco_id / 5.4d-3 - tape_width = scaling * 3.75d-3 - ! Properties of a single strand - tape_thickness = rebco_thickness + copper_thick + hastelloy_thickness - stack_thickness = sqrt(croco_id**2 - tape_width**2) - tapes = stack_thickness / tape_thickness - - copper_area = pi * croco_thick * d - pi * croco_thick**2 & ! copper tube - + copper_thick*tape_width*tapes ! copper in tape - hastelloy_area = hastelloy_thickness * tape_width * tapes - solder_area = pi / 4.0d0 * croco_id**2 - stack_thickness * tape_width - - rebco_area = rebco_thickness * tape_width * tapes - croco_strand_area = pi / 4.0d0 * d**2 - croco_strand_critical_current = jcritsc * rebco_area - - ! Conductor properties - !conductor%number_croco = conductor%acs*(1d0-cable_helium_fraction-copper_bar)/croco_strand_area - conductor_critical_current = croco_strand_critical_current * 6.0d0 - ! Area of core = area of strand - conductor_copper_bar_area = croco_strand_area - conductor_copper_area = copper_area * 6.0d0 + conductor_copper_bar_area - conductor_copper_fraction = conductor_copper_area / conductor_area - - ! Helium area is set by the user. - !conductor_helium_area = cable_helium_fraction * conductor_acs - conductor_helium_area = pi / 2.0d0 * d**2 - conductor_helium_fraction = conductor_helium_area / conductor_area - - conductor_hastelloy_area = hastelloy_area * 6.0d0 - conductor_hastelloy_fraction = conductor_hastelloy_area / conductor_area - - conductor_solder_area = solder_area * 6.0d0 - conductor_solder_fraction = conductor_solder_area / conductor_area - - conductor_rebco_area = rebco_area * 6.0d0 - conductor_rebco_fraction = conductor_rebco_area / conductor_area - -end subroutine croco - -subroutine copper_properties(T,B, copper_resistivity, copper_rrr, copper_cp) - ! Review of ROXIE's Material Properties Database for Quench Simulation, - ! Author: Giulio Manfreda, December 2011 - ! https://espace.cern.ch/roxie/Documentation/Materials.pdf - ! Different models use different definitions for residual resistivity ratio RRR. - ! CUDI: resistivity at 290 K / 4 K. - ! The range of validity of this t is between 4 K and 300 K. - implicit none - - real(dp), intent(in) :: T, B ! temperature, field - real(dp), intent(out) :: copper_resistivity - real(dp), intent(in) :: copper_rrr - real(dp), intent(out) :: copper_cp - real(dp):: bracket, logt, sum - ! Fitting constants: resistivity - real(dp), parameter:: t5 = 2.32547d9 - real(dp), parameter:: t3 = 9.57137d5 - real(dp), parameter:: t1 = 1.62735d2 - real(dp), parameter:: mr = 5.0d-11 ! ohm.m/T - real(dp), parameter:: a = 1.7 ! ohm.m - ! Fitting constants: specific heat (p.13) - ! Checked against - ! http://cryogenics.nist.gov/MPropsMAY/OFHC%20Copper/OFHC_Copper_rev1.htm - real(dp), parameter::a0 = -1.91844d0 - real(dp), parameter::a1 = -0.15973d0 - real(dp), parameter::a2 = 8.61013d0 - real(dp), parameter::a3 = -18.996d0 - real(dp), parameter::a4 = 21.9661d0 - real(dp), parameter::a5 = -12.7328d0 - real(dp), parameter::a6 = 3.54322d0 - real(dp), parameter::a7 = -0.3797d0 - - ! page 5: Copper resistivity is computed in CUDI with the t function similar - ! to the one of McAshan [McA88] - ! Note this formula is much quicker to evaluate than the NIST formula. - - bracket = 1d0 / (t5/T**5 + t3/T**3 + t1/T) - - copper_resistivity = 1.d-8 * (a/copper_rrr + bracket) + mr*B - - ! Specific heat - ! NIST typical polynomial interpolation, equation 4 page 3 - logt = log10(T) - sum = a0 + a1*logt + a2*logt**2 + a3*logt**3 + a4*logt**4 + a5*logt**5 + a6*logt**6 +a7*logt**7 - copper_cp = 10**sum - -end subroutine copper_properties -! ------------------------------------------------------------------------- -subroutine hastelloy_properties(temperature, hastelloy_cp) - implicit none - - real(dp), intent(in) :: temperature ! temperature - real(dp), intent(out) :: hastelloy_cp - real(dp) :: T - - T = temperature - if(temperature>300d0) T=300.0d0 - - ! Reinhard Heller: obtained by fitting data published in : - ! J. Lu, E. S. Choi, and H. D. Zhou, - ! "Physical properties of Hastelloy C-276 at cryogenic temperatures," - ! J. Appl. Phys. 103(6) 2008 064908 - - if(T<42.2571d0)then - hastelloy_cp = 0.60796d0 + 0.15309d0*T - 0.00237d0*T**2 + 6.76732d-4*T**3 - else if(T.ge.42.2571d0)then - hastelloy_cp = -147.06251d0 + 5.43432d0*T - 0.01937d0*T**2 + 2.71669d-5*T**3 & - -8.12438d-9*T**4 - endif - - ! if(T<48.7135d0)then - ! hastelloy%resistivity = 1.23386d-6 - 1.40462d-9*T + 1.09943d-10*T**2 - & - ! 3.81875d-12*T**3 + 6.3866d-14*T**4 - 4.10322d-16*T**5 - ! else if(T.ge.48.7135d0)then - ! hastelloy%resistivity = 1.22641d-6 + 1.19188d-10*T - ! endif -end subroutine hastelloy_properties -! -------------------------------------------------------------------------- -subroutine solder_properties(T, solder_cp) - use constants, only: pi - implicit none - - real(dp), intent(in) :: T ! temperature - real(dp), intent(out) :: solder_cp - - ! Reinhard Heller: obtained by fitting data - ! Material Database from Cryodata Software Package, CRYOCOMP, version 3.0, Florence, SC - - if(T<20.0d0)then - solder_cp = 4.88028d0 - 2.92865d0*T + 0.52736d0*T**2 - 0.01861d0*T**3 + 2.36019d-4*T**4 - else if((T.ge.20.0d0).and.(T<300.0d0))then - solder_cp = -10.15269d0 + 3.70087d0*T - 0.02947d0*T**2 + 1.02933d-4*T**3 - 1.29518d-7*T**4 - else - solder_cp = 181.29d0 - endif -end subroutine solder_properties -! ------------------------------------------------------------------------- -subroutine jacket_properties(T, jacket_cp) - implicit none - - real(dp), intent(in) :: T ! temperature - real(dp), intent(out) :: jacket_cp - real(dp):: logt, sum - ! Fitting constants: specific heat (p.13) - ! http://cryogenics.nist.gov/MPropsMAY/304Stainless/304Stainless_rev.htm - real(dp), parameter::a0 = 22.0061 - real(dp), parameter::a1 = -127.5528 - real(dp), parameter::a2 = 303.647 - real(dp), parameter::a3 = -381.0098 - real(dp), parameter::a4 = 274.0328 - real(dp), parameter::a5 = -112.9212 - real(dp), parameter::a6 = 24.7593 - real(dp), parameter::a7 = -2.239153 - - logt = log10(T) - if(T>300) logt = log10(300d0) - sum = a0 + a1*logt + a2*logt**2 + a3*logt**3 + a4*logt**4 + a5*logt**5 + & - a6*logt**6 + a7*logt**7 - jacket_cp = 10**sum - -end subroutine jacket_properties -! ---------------------------------------------------------------------------- -subroutine helium_properties(T, helium_cp_density) - implicit none - - ! Isobaric Data for P = 0.60000 MPa - ! http://webbook.nist.gov/chemistry/fluid/ - ! See very approximate fits in quench_data.xlsx (Issue #522) - - real(dp), intent(in) :: T ! temperature - real(dp), intent(out) :: helium_cp_density - ! Cp x density J/K/m3 - if((T>=4d0).and.(T<10d0))then - helium_cp_density = 12.285d3*T**3 - 309.92d3*T**2 + 2394.6d3*T - 5044.8d3 - else if(T>=10d0)then - ! This probably works OK for arbitrarily high temperature - helium_cp_density = 1745.1d3*T**(-1.031d0) - else - write(*,*)'temperature is below the range of helium data fit (>4 K): ', T - end if - -end subroutine helium_properties - -! --------------------------------------------------------------------------- -!##################################################################### -! -! Sn-40Pb PROPERTIES PACKAGE -! ------------------------ -! Contains functions for the calculation of the thermo-physical -! properties of Sn-40Pb -! -!##################################################################### - -!##################################################################### -real function dSn40Pb(T) - !##################################################################### - ! - ! Density of Sn-40Pb - ! - ! Range: 0 <= T <= inf K - ! - ! References - ! ---------- - ! Cryocomp v3.0 - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! dSn40Pb x density Kg/m**3 - ! - ! - ! Author : L.Bottura at Cryosoft - ! Version: 1.0 June 2015 - ! - !##################################################################### - implicit none - ! external variables - real T - ! fit variables - - ! local variables - dSn40Pb = 8400.0 - - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function dSn40Pb -!##################################################################### -real function cSn40Pb(T) - !##################################################################### - ! - ! Specific heat of Sn-40Pb - ! - ! Range: 1 <= T <= 300 K - ! - ! References - ! ---------- - ! MISSING - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! cSn40Pb x specific heat J/Kg K - ! - ! - ! Author : L.Bottura at Cryosoft - ! Version: 1.0 June 2015 - ! - !##################################################################### - implicit none - ! external variables - real T - ! fit variables - real Tmin,Tmax,T0 - real AA,BB,CC,DD,a,b,c,d,na,nb,nc,nd - - data Tmin / 1.0 /, Tmax / 300.0/ - data T0 / 3.96249314 / - - save - ! local variables - real TT - - ! - TT=T - TT=min(TT,Tmax) - TT=max(TT,Tmin) - - if(TT.ge.1.and.TT.le.T0) then - - AA = -0.07848998 - BB = 0.2150118 - CC = -0.42863307 - DD = 0.06224209 - a = 1.89412709 - b = 556.143928 - c = 0.20145037 - d = 0.00086726 - na = -0.10237772 - nb = -0.06338306 - nc = 0.65668303 - nd = 0.31793398 - - else - - AA = -0.0002679 - BB = -0.39223444 - CC = 1.11558347 - DD = -1.06506974 - a = 18.7914061 - b = 43.8271896 - c = 6.10665608 - d = 3.14878533 - na = 4.50255658 - nb = 1.22694566 - nc = 2.24090761 - nd = 3.38478332 - - endif - - cSn40Pb = AA*TT**1/(1+TT/a)**na + BB*TT**2/(1+TT/b)**nb + & - CC*TT**3/(1+TT/c)**nc + DD*TT**4/(1+TT/d)**nd - - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function cSn40Pb -!##################################################################### -real function kSn40Pb(T) - !##################################################################### - ! - ! Thermal conductivity of Sn-40Pb - ! - ! Range: 2 <= T <= 300 K - ! - ! References - ! ---------- - ! MISSING - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! kSn40Pb x thermal conductivity W/m K - ! - ! - ! Author : L.Bottura at Cryosoft - ! Version: 1.0 June 2015 - ! - !##################################################################### - implicit none - ! external variables - real T - ! fit variables - real Tmin,Tmax,T0 - real rho0,rho273,RRR - real alpha,alpha2,beta - real L0,K0 - real m,n,arg1,arg2 - real W0,Wi,Wi0 - real l - - data Tmin / 2.0/, Tmax / 300.0/ - data T0 / 22.702953 / - data rho273 / 1.50E-07 / , RRR / 26.4214669 / - data alpha2 / 2.1156885E-08 / - data L0 / 2.443E-08 / , K0 / 5.74E+01 / - data m / 3.092523 /, n / 3.0255307 / - data l / 0.20328615 / - save - ! local variables - real TT - - ! - TT=T - TT=min(TT,Tmax) - TT=max(TT,Tmin) - - rho0 = rho273/(RRR-1) - beta = rho0/L0 - arg1 = (m-n)*(m+l) - arg2 = (TT-T0)/T0 - alpha = alpha2*(beta/n/alpha2)**arg1 - W0 = beta/TT - Wi = alpha*TT**n - Wi0 = Wi+W0 - - if(TT.le.T0) then - - kSn40Pb = 1/Wi0 - - else - kSn40Pb = 1/Wi0 + K0*(1-exp(-arg2)) - endif - - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function kSn40Pb -!##################################################################### -real function rSn40Pb(T) - !##################################################################### - ! - ! Electrical resistivity of Sn-40Pb - ! - ! Range: 0 <= T <= 300 K - ! - ! References - ! ---------- - ! MISSING - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! rSn40Pb x resistivity Ohm m - ! - ! Author : L.Bottura at CryoSoft - ! Version: 1.0 June 2015 - ! - !##################################################################### - implicit none - ! external variables - real T - ! fit variables - real RRR,rho273,rho0,rhoi,rhoi0 - real p1,p2,p3,p4,p5,p6,p7 - real arg - real Tmin,Tmax - - data RRR / 45 /, rho273 / 1.50E-07 / - data p1 / 3.98307546E-14 /, p2 / 3.47905728 /, & - p3 / 3.6576564E+9 /, p4 / -1.10 /, & - p5 / 0.0 /, p6 / 1.0 / , & - p7 / 0.333 / - data Tmin / 0.0/, Tmax / 300.0/ - save - ! local variables - real TT - - ! - TT=T - TT=min(TT,Tmax) - TT=max(TT,Tmin) - - rho0 = rho273/(RRR-1) - arg = (p5/TT)**p6 - rhoi = p1*TT**p2/(1.+p1*p3*TT**(p2+p4)*exp(-arg)) - rhoi0 = (p7*rhoi*rho0)/(rhoi+rho0) - - rSn40Pb = rho0+rhoi+rhoi0 - - ! - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function rSn40Pb - -!##################################################################### -! -! HASTELLOY C276 PROPERTIES PACKAGE -! --------------------------------- -! Contains functions for the calculation of the thermo-physical -! properties of Hastelloy C276 (Haynes superalloy) -! -!##################################################################### - -!##################################################################### -real function dHastelloyC276(T) - !##################################################################### - ! - ! Density of Hastelloy C276 - ! - ! Range: 0 <= T <= inf K - ! - ! References - ! ---------- - ! J.Lu, E.S. Choi, H.D. Zhou, Physical Properties of Hastelloy(R) - ! C-276 at cryogenic temperatures, J. Appl. Phys, 103, 064908, 2008 - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! dHastelloyC276 x density Kg/m**3 - ! - ! - ! Author : L.Bottura at Cryosoft - ! Version: 1.0 December 2012 - ! - !##################################################################### - implicit none - ! external variables - real T - ! fit variables - - ! local variables - - dHastelloyC276 = 8890.0 - - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function dHastelloyC276 -!##################################################################### -real(dp) function cHastelloyC276(T) - !##################################################################### - ! - ! Specific heat (not density) of Hastelloy C276 - ! - ! Range: 1 <= T <= 300 K - ! - ! References - ! ---------- - ! J.Lu, E.S. Choi, H.D. Zhou, Physical Properties of Hastelloy(R) - ! C-276 at cryogenic temperatures, J. Appl. Phys, 103, 064908, 2008 - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! cHastelloyC276 x specific heat J/Kg K - ! - ! - ! Author : L.Bottura at Cryosoft - ! Version: 1.0 December 2012 - ! - !##################################################################### - implicit none - ! external variables - real(dp) T - ! fit variables - real(dp) AA,BB,CC,DD,a,b,c,d,na,nb,nc,nd - real(dp) Tmin,Tmax - data AA / 18.46314493 / , BB / 1298.042986 / ,& - CC /-1105.076534 / , DD / 2.226310361 / - data a / 2.634486768 / , b / 13.0796954 / ,& - c / 23.29856632 / , d / 0.932886135 / - data na / 14.010908 / , nb / 12.68689825 / ,& - nc / 14.64960608 / , nd / 2.626448071 / - data Tmin / 1.0/, Tmax / 300.0/ - save - ! local variables - real(dp) TT - - TT=T - TT=min(TT,Tmax) - TT=max(TT,Tmin) - - cHastelloyC276 = AA*TT**na/(a+TT)**na + BB*TT**nb/(b+TT)**nb + & - CC*TT**nc/(c+TT)**nc + DD*TT**nd/(d+TT)**nd - - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function cHastelloyC276 -!##################################################################### -real function kHastelloyC276(T) - !##################################################################### - ! - ! Thermal conductivity of Hastelloy C276 - ! - ! Range: 2 <= T <= 300 K - ! - ! References - ! ---------- - ! J.Lu, E.S. Choi, H.D. Zhou, Physical Properties of Hastelloy(R) - ! C-276 at cryogenic temperatures, J. Appl. Phys, 103, 064908, 2008 - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! kHastelloyC276 x thermal conductivity W/m K - ! - ! - ! Author : L.Bottura at Cryosoft - ! Version: 1.0 December 2012 - ! - !##################################################################### - implicit none - ! external variables - real T - ! fit variables - real p1,p2,p3,p4 - real Tmin,Tmax - data p1 / 0.2093 / , p2 / 125.7155 / ,& - p3 / 60.0344 / , p4 / 2.1669 / - data Tmin / 2.0/, Tmax / 400.0/ - save - ! local variables - real TT - - TT=T - TT=min(TT,Tmax) - TT=max(TT,Tmin) - kHastelloyC276 = p1*TT*(1+(TT/p2)**p4) / (1+(TT/p3)**p4) - - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function kHastelloyC276 -!##################################################################### -real function rHastelloyC276(T) - !##################################################################### - ! - ! Electrical resistivity of Hastelloy C276 - ! - ! Range: 0 <= T <= 400 K - ! - ! References - ! ---------- - ! J.Lu, E.S. Choi, H.D. Zhou, Physical Properties of Hastelloy(R) - ! C-276 at cryogenic temperatures, J. Appl. Phys, 103, 064908, 2008 - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! rHastelloyC276 x resistivity Ohm m - ! - ! Author : L.Bottura at CryoSoft - ! Version: 1.0 October 2012 - ! - !##################################################################### - implicit none - ! external variables - real T - ! fit variables - real A,B,C - real Tmin,Tmax - data A / 1.22634e-6 / , B / 1.26121e-10 / ,& - C / -2.2417e-14 / - data Tmin / 0.0/, Tmax / 400.0/ - save - ! local variables - real TT - - TT=T - TT=min(TT,Tmax) - TT=max(TT,Tmin) - rHastelloyC276 = A + B*TT + C*TT*TT - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function rHastelloyC276 -!##################################################################### -! -! COPPER PROPERTIES PACKAGE -! ------------------------- -! -! Contains functions for the calculation of the thermo-physical -! properties of pure Copper (Cu) -! -!##################################################################### - -!##################################################################### -real function dCu(T) - !##################################################################### - ! - ! Density of Copper - ! - ! Range: 0 <= T <= inf K - ! - ! References - ! ---------- - ! http://en.wikipedia.org/wiki/Copper - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! dCu x density Kg/m**3 - ! - ! - ! Author : L.Bottura at Cryosoft - ! Version: 1.0 October 2012 - ! - !##################################################################### - implicit none - ! external variables - real T - ! fit variables - - ! local variables - - ! - dCu = 8960.0 - ! - - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function dCu -!##################################################################### -real function cCu(T) - !##################################################################### - ! - ! Specific heat of Copper - ! - ! Range: 1 <= T <= 1000 K - ! - ! References - ! ---------- - ! CryoComp version 3.0 - ! Handbook on Materials for S.C. Machinery, NBS Boulder (yellow - ! book), 1977 - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! cCu x specific heat J/Kg K - ! - ! Author : L.Bottura & C. Rosso at Cryosoft - ! Version: 3.5 October 2012 - ! - !##################################################################### - implicit none - ! external variables - real T - ! fit variables - real T0,T1 - real Tmin,Tmax - real a1,a2,a3 - real b0,b1,b2,b3,b4 - real AA,CC,DD,a,c,d,na,nc,nd - data T0 / 10.4529369 / , T1 / 48.26583891 / - data a1 / 0.01188007 / , a2 / -0.00050323 / , & - & a3 / 0.00075762 / - data b0 / -5.03283229 / , b1 / 1.27426834 / , & - & b2 / -0.11610718 / , b3 / 0.00522402 / ,& - & b4 / -5.2996E-05 / - data AA / -65.07570094 / , & - & CC / 624.7552517 / , DD / 0.529512119 / - data a / 1.833505318 / , & - & c / 16.55124429 / , d / -0.000101401 / - data na / 0.518553624 / ,& - & nc / 2.855560719 / , nd / 2.983928329 / - data Tmin / 1.0/, Tmax / 1000.0/ - save - - ! local variables - real TT - - ! - TT=T - TT=max(TT,Tmin) - TT=min(TT,Tmax) - - if (TT.le.t0) then - cCu = a1*TT + a2*TT**2 + a3*TT**3 - elseif (TT.gt.t0 .and. TT.le.t1) then - cCu = b0 + b1*TT+ b2*TT**2 + b3*TT**3 + b4*TT**4 - elseif(TT.gt.t1) then - cCu = AA*TT /(a+TT)**na +& - & CC*TT**3/(c+TT)**nc + DD*TT**4/(d+TT)**nd - endif - - ! - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function cCu -!##################################################################### -real function kCu(T,B,RRR) - !##################################################################### - ! - ! Thermal conductivity of Copper - ! - ! Range: 0.1 <= T <= 1000 K, 0 <= B <= 30 T, 1.5 <= RRR <= 3000 - ! - ! References - ! ---------- - ! CryoComp version 3.0 - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! B x magnetic field T - ! RRR x residual resistivity ratio - - ! kCu x thermal conductivity W/m K - ! - ! Author : L.Bottura & C. Rosso at CryoSoft - ! Version: 3.5 October 2012 - ! - !##################################################################### - implicit none - ! external variables - real T,B,RRR - ! fit variables - real Tmin,Tmax - parameter(Tmin = 0.1, Tmax = 1000.0) - real Bmin,Bmax - parameter(Bmin = 0.0, Bmax = 30.0) - real RRRmin,RRRmax - parameter(RRRmin = 1.5, RRRmax = 3000.0) - real rho273 - parameter(rho273=1.54e-8) - real L0 - parameter(L0=2.443e-8) - real p1,p2,p3,p4,p5,p6 - parameter(p1=1.754e-8, p2=2.763, p3=1102.,& - & p4=-0.165 , p5=70.0 , p6=1.765) - ! local variables - real TT,BB,R - real rhozero,beta,p7,arga,argb,argc,argd - real wt,wi0,wi,w0,wc,magr - !real magrCu - - ! - TT=T - TT=max(TT,Tmin) - TT=min(TT,Tmax) - - BB=B - BB=max(BB,Bmin) - BB=min(BB,Bmax) - - R =RRR - R =max(R,RRRmin) - R =min(R,RRRmax) - - ! - rhozero = rho273/(R-1.) - - beta = rhozero/L0 - - p7 = 0.838/(beta/0.0003)**0.1661 - arga = min((log(TT/470.0)/0.70)**2,30.0) - argb = min((log(TT/ 87.0)/0.45)**2,30.0) - argc = min((log(TT/ 21.0)/0.50)**2,30.0) - argd = min((p5/TT)**p6,30.0) - - wc = -0.00012*log(TT/420.0)*exp(-arga) -& - & 0.00016*log(TT/ 73.0)*exp(-argb) -& - & 0.00002*log(TT/ 18.0)*exp(-argc) - w0 = beta/TT - wi = p1*TT**p2/(1.+p1*p3*TT**(p2+p4)*exp(-argd))+wc - wi0 = p7*wi*w0/(wi+w0) - wt = w0+wi+wi0 - - magr = magrCu(TT,BB,R) - - kCu = 1.0/(wt*magr) - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function kCu -!##################################################################### -real function rCu(T,B,RRR) - !##################################################################### - ! - ! Electrical resistivity of Copper - ! - ! Range: 0.1 <= T <= 1000 K, 0 <= B <= 30 T, 1.5 <= RRR <= 3000 - ! - ! References - ! ---------- - ! CryoComp version 3.0 - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! B x magnetic field T - ! RRR x residual resistivity ratio - - ! rCu x resistivity Ohm m - ! - ! Author : L.Bottura & C. Rosso at CryoSoft - ! Version: 3.5 October 2012 - ! - !##################################################################### - implicit none - ! external variables - real T,B,RRR - ! fit variables - real Tmin,Tmax - parameter(Tmin = 0.1, Tmax = 1000.0) - real Bmin,Bmax - parameter(Bmin = 0.0, Bmax = 30.0) - real RRRmin,RRRmax - parameter(RRRmin = 1.5, RRRmax = 3000.0) - real rho273 - parameter(rho273=1.54e-8) - real p1,p2,p3,p4,p5,p6,p7 - parameter(p1=0.1171e-16, p2=4.49, p3=3.841e10,& - & p4=-1.14, p5=50.0, p6=6.428, p7=0.4531) - ! local variables - real TT,BB,R,rhozero,arg,rhoi,rhoi0,rho0,magr - !real magrCu - - TT=T - TT=max(TT,Tmin) - TT=min(TT,Tmax) - - BB=B - BB=max(BB,Bmin) - BB=min(BB,Bmax) - - R =RRR - R =max(R,RRRmin) - R =min(R,RRRmax) - - ! resistivity at absolute zero - rhozero = rho273/(R-1.0) - - ! resistivity at the temperature t - arg = min((p5/TT)**p6,30.0) - rhoi = p1*TT**p2/(1.+p1*p3*TT**(p2+p4)*exp(-arg)) - rhoi0 = p7*rhoi*rhozero/(rhoi+rhozero) - rho0 = rhozero+rhoi+rhoi0 - - ! transverse magneto-resistance factor - magr = magrCu(TT,BB,R) - - ! resistivity - rCu = magr * rho0 - - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function rCu -!##################################################################### -! -! Auxiliary functions and calculations -! -!##################################################################### - -!##################################################################### -real function magrCu(T,B,RRR) - !##################################################################### - ! - ! Magneto-resistivity factor of Copper, given as the ratio between - ! resistivity in transverse magnetic field B to resistivity at zero - ! field - ! - ! Range: 0.1 <= T <= 1000 K, 0 <= B <= 30 T, 1.5 <= RRR <= 3000 - ! - ! References - ! ---------- - ! CryoComp version 3.0 - ! - ! variable I/O meaning units - ! -------------------------------------------------------------------- - ! T x absolute temperature K - ! B x magnetic field T - ! RRR x residual resistivity ratio - - ! magrCu x magneto-resistivity factor - - ! - ! Author : L.Bottura and C. Rosso at Cryosoft - ! Version: 1.1 October 2012 - ! - !##################################################################### - implicit none - ! external variables - real T,B,RRR - ! fit variables - real Tmin,Tmax - parameter(Tmin = 0.1, Tmax = 1000.0) - real RRRmin,RRRmax - parameter(RRRmin = 1.5, RRRmax = 3000.0) - real Bmin,Bmax - parameter(Bmin = 0.0, Bmax = 30.0) - real brrmin,brrmax - parameter(brrmin = 0.0, brrmax = 40.0e3) - real rho273,rhorrr - parameter(rho273=1.54e-8,rhorrr=2.37e-8) - real p1,p2,p3,p4,p5,p6,p7 - parameter(p1=0.1171e-16, p2=4.49, p3=3.841e10,& - & p4=-1.14, p5=50.0, p6=6.428, p7=0.4531) - real a1,a2,a3,a4 - parameter(a1=0.382806e-03, a2=0.132407e+01, a3=0.167634e-02,& - & a4=0.789953e+00) - ! local variables - real TT,BB,R,rhozero,rhoice,arg,rhoi,rhoi0,rho0,brr,magr - - TT=T - TT=max(TT,Tmin) - TT=min(TT,Tmax) - - BB=B - BB=max(BB,Bmin) - BB=min(BB,Bmax) - - R =RRR - R =max(R,RRRmin) - R =min(R,RRRmax) - - ! resistivity at absolute zero - rhozero = rho273/(R-1.0) - - ! fit for the resistivity at the ice temperature - rhoice = rho273 + rhorrr/R - - ! resistivity at the temperature t - arg = min((p5/TT)**p6,30.0) - rhoi = p1*TT**p2/(1.+p1*p3*TT**(p2+p4)*exp(-arg)) - rhoi0 = p7*rhoi*rhozero/(rhoi+rhozero) - rho0 = rhozero+rhoi+rhoi0 - - ! product of field and residual resistivity ratio - brr = BB*rhoice/rho0 - brr = max(brr,brrmin) - brr = min(brr,brrmax) - - ! fit for the transverse magneto-resistance increase - if(brr .gt. 1.) then - magr = a1*brr**a2/(1.0+a3*brr**a4) - else - magr = 0.0 - endif - - ! transverse magneto-resistance factor - magrCu = magr+1.0 - - return -! SJP Issue #835 -! For Intel compliance add "end function" -end function magrCu - -!----------------------------------------------------------------- - - -end module superconductors diff --git a/tests/unit/test_superconductors.py b/tests/unit/test_superconductors.py new file mode 100644 index 0000000000..3fa3df587c --- /dev/null +++ b/tests/unit/test_superconductors.py @@ -0,0 +1,192 @@ +import pytest +from typing import NamedTuple, Any + +import process.superconductors as superconductors + + +class IterscParam(NamedTuple): + thelium: Any = None + + bmax: Any = None + + strain: Any = None + + bc20max: Any = None + + tc0max: Any = None + + expected_jcrit: Any = None + + expected_bcrit: Any = None + + expected_tcrit: Any = None + + +@pytest.mark.parametrize( + "iterscparam", + ( + IterscParam( + thelium=4.75, + bmax=13.008974843466492, + strain=0.001601605753441172, + bc20max=32.969999999999999, + tc0max=16.059999999999999, + expected_jcrit=692348194.774593, + expected_bcrit=27.092853296363597, + expected_tcrit=11.338458919718571, + ), + IterscParam( + thelium=6.2510000000000003, + bmax=13.008974843466492, + strain=0.001601605753441172, + bc20max=32.969999999999999, + tc0max=16.059999999999999, + expected_jcrit=495889332.08959526, + expected_bcrit=24.442648486388464, + expected_tcrit=11.338458919718571, + ), + ), +) +def test_itersc(iterscparam): + """ + Automatically generated Regression Unit Test for itersc. + + This test was generated using data from tests/regression/scenarios/large-tokamak/IN.DAT. + + :param iterscparam: the data used to mock and assert in this test. + :type iterscparam: iterscparam + + :param monkeypatch: pytest fixture used to mock module/class variables + :type monkeypatch: _pytest.monkeypatch.monkeypatch + """ + + jcrit, bcrit, tcrit = superconductors.itersc( + thelium=iterscparam.thelium, + bmax=iterscparam.bmax, + strain=iterscparam.strain, + bc20max=iterscparam.bc20max, + tc0max=iterscparam.tc0max, + ) + + assert jcrit == pytest.approx(iterscparam.expected_jcrit) + + assert bcrit == pytest.approx(iterscparam.expected_bcrit) + + assert tcrit == pytest.approx(iterscparam.expected_tcrit) + + +class JcritNbtiParam(NamedTuple): + temperature: Any = None + + bmax: Any = None + + c0: Any = None + + bc20max: Any = None + + tc0max: Any = None + + expected_jcrit: Any = None + + expected_tcrit: Any = None + + +@pytest.mark.parametrize( + "jcritnbtiparam", + ( + JcritNbtiParam( + temperature=4.75, + bmax=8.0517923638507547, + c0=10000000000, + bc20max=15, + tc0max=9.3000000000000007, + expected_jcrit=906668274.04561484, + expected_tcrit=5.9060082696285683, + ), + JcritNbtiParam( + temperature=6, + bmax=8.0517923638507547, + c0=10000000000, + bc20max=15, + tc0max=9.3000000000000007, + expected_jcrit=-73718607.547511846, + expected_tcrit=5.9060082696285683, + ), + ), +) +def test_jcrit_nbti(jcritnbtiparam): + """ + Automatically generated Regression Unit Test for jcrit_nbti. + + This test was generated using data from tests/regression/scenarios/large-tokamak/IN.DAT. + + :param jcritnbtiparam: the data used to mock and assert in this test. + :type jcritnbtiparam: jcritnbtiparam + + :param monkeypatch: pytest fixture used to mock module/class variables + :type monkeypatch: _pytest.monkeypatch.monkeypatch + """ + + jcrit, tcrit = superconductors.jcrit_nbti( + temperature=jcritnbtiparam.temperature, + bmax=jcritnbtiparam.bmax, + c0=jcritnbtiparam.c0, + bc20max=jcritnbtiparam.bc20max, + tc0max=jcritnbtiparam.tc0max, + ) + + assert jcrit == pytest.approx(jcritnbtiparam.expected_jcrit) + + assert tcrit == pytest.approx(jcritnbtiparam.expected_tcrit) + + +def test_jcrit_rebco(): + jcrit_rebco, validity = superconductors.jcrit_rebco(4.75, 7.0) + + assert jcrit_rebco == pytest.approx(55870234414.171684) + assert validity + + +def test_current_sharing_rebco(): + assert superconductors.current_sharing_rebco(7.0, 2e7) == pytest.approx( + 75.76286550648135 + ) + + +def test_bi2212(): + jcrit, tmarg = superconductors.bi2212(7.0, 2e7, 4.75, 0.2) + + assert jcrit == pytest.approx(174017403.16041547) + assert tmarg == pytest.approx(13.750991122745397) + + +def test_gl_nbti(): + jcrit, tcrit, bcrit = superconductors.gl_nbti(4.75, 7.0, 2, 9.5, 13.75) + + assert jcrit == pytest.approx(2551683055.6511745) + assert tcrit == pytest.approx(7.277374792835339) + assert bcrit == pytest.approx(13.662161361675887) + + +def test_wstsc(): + jcrit, bcrit, tcrit = superconductors.wstsc(4.75, 27.0, 0.001, 30.0, 25.0) + + assert jcrit == pytest.approx(195513.0673058944) + assert bcrit == pytest.approx(27.329369840368482) + assert tcrit == pytest.approx(5.170678992915718) + + +def test_gl_rebco(): + jcrit, bcrit, tcrit = superconductors.gl_rebco(4.75, 7.0, 2, 30.0, 25.0) + + assert jcrit == pytest.approx(14527.765708690296) + assert bcrit == pytest.approx(9.439350824747793) + assert tcrit == pytest.approx(24.66989137698065) + + +def test_hijc_rebco(): + jcrit, bcrit, tcrit = superconductors.hijc_rebco(4.75, 7.0, 2, 30.0, 25.0) + + assert jcrit == pytest.approx(44418407.919617616) + assert bcrit == pytest.approx(22.335736687814954) + assert tcrit == pytest.approx(24.999125)