From aec213e73f5af6bb480f4fa24794213e4a0cda76 Mon Sep 17 00:00:00 2001 From: apearce Date: Wed, 30 Oct 2024 11:00:24 +0000 Subject: [PATCH 1/8] =?UTF-8?q?=F0=9F=90=9B=20update=20VV=20thickness=20to?= =?UTF-8?q?=20be=20double=20shell=20thickness=20in=20quench=20model?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/sctfcoil.py | 4 ++-- source/fortran/build_variables.f90 | 4 ++++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/process/sctfcoil.py b/process/sctfcoil.py index 0b70aa3b..aced6c75 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -2337,7 +2337,7 @@ def vv_stress_on_quench(self): # TODO: is this the correct current? I_op=sctfcoil_module.tfc_current / tfcoil_variables.n_tf_turn, # VV properties - d_vv=build_variables.d_vv_in, + d_vv=build_variables.d_vv_shell_thickness, ) def tf_field_and_force(self): @@ -7248,7 +7248,7 @@ def vv_stress_on_quench( :param taud: the discharge time of the TF coil when quench occurs :param I_op: the 'normal' operating current of the TF coil - :param d_vv: the thickness of the vacuum vessel + :param d_vv: the thickness of the vacuum vessel shell :returns: the maximum stress experienced by the vacuum vessel diff --git a/source/fortran/build_variables.f90 b/source/fortran/build_variables.f90 index 843ad5f6..bf6957cc 100644 --- a/source/fortran/build_variables.f90 +++ b/source/fortran/build_variables.f90 @@ -78,6 +78,9 @@ module build_variables real(dp) :: d_vv_bot !! vacuum vessel underside thickness (TF coil / shield) (m) + real(dp) :: d_vv_shell_thickness + !! vacuum vessel double walled shell thickness (m) + real(dp) :: f_avspace !! F-value for stellarator radial space check (`constraint equation 83`) @@ -321,6 +324,7 @@ subroutine init_build_variables d_vv_out = 0.07D0 d_vv_top = 0.07D0 d_vv_bot = 0.07D0 + d_vv_shell_thickness = 0.12D0 f_avspace = 1.0D0 fcspc = 0.6D0 fseppc = 3.5D8 From 3d3b25e0856fe6ca65b8bcf344a8d4c27c00418b Mon Sep 17 00:00:00 2001 From: apearce Date: Wed, 30 Oct 2024 11:13:03 +0000 Subject: [PATCH 2/8] =?UTF-8?q?=E2=9C=8F=EF=B8=8F=20update=20doc=20string?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/sctfcoil.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/process/sctfcoil.py b/process/sctfcoil.py index aced6c75..ad663761 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -7257,8 +7257,8 @@ def vv_stress_on_quench( The theta1 quantity for the TF coil and VV is not very meaningful. The impact of it of the inductance is rather small. Generally, the paper seems to suggest the TF coil is between 40 and 60, as this is the range they calculate - the surrogates over. No range is provided for the VV but the example using - JA DEMO is 1 degree suggesting the quantity will be very small. + the surrogates over. The thickness of the VV considers an ITER like design and + only the outer and inner shells that which act of conductuve structural material. References ---------- From f38afc8ac59c849d34a64d0cc08e915bdddfe69d Mon Sep 17 00:00:00 2001 From: apearce Date: Wed, 30 Oct 2024 11:39:29 +0000 Subject: [PATCH 3/8] =?UTF-8?q?=F0=9F=90=9B=20include=20size=20tf=20cases?= =?UTF-8?q?=20in=20coil=20case=20cross-sections?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/sctfcoil.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/process/sctfcoil.py b/process/sctfcoil.py index ad663761..4a2d6d6f 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -2331,8 +2331,9 @@ def vv_stress_on_quench(self): # Area of the radial plate taken to be the area of steel in the WP # TODO: value clipped due to #1883 S_rp=numpy.clip(sctfcoil_module.a_tf_steel, 0, None), - # TODO: Does this calculation of Scc exclude the area of the case down the side? - S_cc=sctfcoil_module.a_case_front + sctfcoil_module.a_case_nose, + S_cc=sctfcoil_module.a_case_front + + sctfcoil_module.a_case_nose + + 2.0 * sctfcoil_module.t_lat_case_av, taud=tfcoil_variables.tdmptf, # TODO: is this the correct current? I_op=sctfcoil_module.tfc_current / tfcoil_variables.n_tf_turn, From c4ad18c7c7af4b166f97a4435c94a3abaf87d87f Mon Sep 17 00:00:00 2001 From: apearce Date: Mon, 4 Nov 2024 16:29:14 +0000 Subject: [PATCH 4/8] =?UTF-8?q?=E2=9C=A8=20add=20integral=20for=20theta=20?= =?UTF-8?q?factor=20into=20VV=20quench=20calculaiton?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/sctfcoil.py | 67 ++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 66 insertions(+), 1 deletion(-) diff --git a/process/sctfcoil.py b/process/sctfcoil.py index 4a2d6d6f..24f35bbd 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -3,6 +3,7 @@ import logging import copy import numba +import numpy as np from process.fortran import rebco_variables from process.fortran import global_variables @@ -7186,6 +7187,69 @@ def _inductance_factor(H, Ri, Ro, Rm, theta1): ) +def lambda_term(tau, omega): + """ + words + """ + p = 1.0 - omega**2.0 + + if p < 0: + integral = (1.0 / np.sqrt(np.abs(p))) * np.arcsin( + (1.0 + omega * tau) / (tau + omega) + ) + else: + integral = (1.0 / np.sqrt(np.abs(p))) * np.log( + (2.0 * (1.0 + tau * omega - np.sqrt(p * (1 - tau**2.0)))) / (tau + omega) + ) + + return integral + + +def _theta_factor_integral(Ro_vv, Ri_vv, Rm_vv, H_vv, theta1_vv): + """ + words + """ + theta2 = np.pi / 2.0 + theta1_vv + a = (Ro_vv - Ri_vv) / 2.0 + Rbar = (Ro_vv + Ri_vv) / 2.0 + A = Rbar / a + delta = (Rbar - Rm_vv) / a + kappa = H_vv / a + iota = (1.0 + delta) / kappa + + denom = np.cos(theta1_vv) + np.sin(theta1_vv) - 1.0 + + R1 = H_vv * ((np.cos(theta1_vv) + iota * (np.sin(theta1_vv) - 1.0)) / denom) + R2 = H_vv * ((np.cos(theta1_vv) - 1.0 + iota * np.sin(theta1_vv)) / denom) + R3 = H_vv * (1 - delta) / kappa + + Rc1 = (H_vv / kappa) * (A + 1.0) - R1 + Rc2 = Rc1 + (R1 - R2) * np.cos(theta1_vv) + Rc3 = Rc2 + Zc2 = (R1 - R2) * np.sin(theta1_vv) + Zc3 = Zc2 + R2 - R3 + + tau = np.array( + [ + [np.cos(theta1_vv), np.cos(theta1_vv + theta2), -1.0], + [1.0, np.cos(theta1_vv), np.cos(theta1_vv + theta2)], + ] + ) + + omega = np.array([Rc1 / R1, Rc2 / R2, Rc3 / R3]) + + # Assume up down symmetry and let Zc6 = - Zc3 + chi1 = (Zc3 + np.abs(-Zc3)) / Ri_vv + chi2 = 0.0 + + for k in range(len(omega)): + chi2 = chi2 + np.abs( + lambda_term(tau[1, k], omega[k]) - lambda_term(tau[0, k], omega[k]) + ) + + return (chi1 + chi2) / np.pi + + def vv_stress_on_quench( # TF shape H_coil, @@ -7268,8 +7332,9 @@ def vv_stress_on_quench( Plasma and Fusion Research. 15. 1405078-1405078. 10.1585/pfr.15.1405078. """ # Poloidal loop resistance (PLR) in ohms + theta_vv = _theta_factor_integral(Ro_vv, Ri_vv, Rm_vv, H_vv, theta1_vv) plr_coil = ((0.5 * ccl_length_coil) / (n_tf * (S_cc + S_rp))) * 1e-6 - plr_vv = ((0.84 / d_vv) * 0.94) * 1e-6 + plr_vv = ((0.84 / d_vv) * theta_vv) * 1e-6 # relevant self-inductances in henry (H) coil_structure_self_inductance = ( From a87fa0f67c11e408d134cd164ed5f39cd6175e86 Mon Sep 17 00:00:00 2001 From: apearce Date: Tue, 5 Nov 2024 14:38:16 +0000 Subject: [PATCH 5/8] =?UTF-8?q?=F0=9F=8E=A8=20add=20doc=20strings?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/sctfcoil.py | 27 +++++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/process/sctfcoil.py b/process/sctfcoil.py index 24f35bbd..b908368d 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -7187,9 +7187,17 @@ def _inductance_factor(H, Ri, Ro, Rm, theta1): ) +@staticmethod def lambda_term(tau, omega): - """ - words + """Lambda Term + + Author: Alexander Pearce, UKAEA + + :param tau: tau_{s,k} = (R_{s,k} - R_{c,k}) / R_k + :param omega: omega_k = R_{c,k}/R_k + + The lammbda fucntion used inegral in inductance calcuation found + in Y. Itoh et al. The full form of the functions are given in appendix A. """ p = 1.0 - omega**2.0 @@ -7205,9 +7213,20 @@ def lambda_term(tau, omega): return integral +@staticmethod def _theta_factor_integral(Ro_vv, Ri_vv, Rm_vv, H_vv, theta1_vv): - """ - words + """Theta Factor Integral + Author: Alexander Pearce, UKAEA + + :param Ro_vv: the radius of the outboard edge of the VV CCL + :param Ri_vv: the radius of the inboard edge of the VV CCL + :param Rm_vv: the radius where the maximum height of the VV CCL occurs + :param H_vv: the maximum height of the VV CCL + :param theta1_vv: the polar angle of the point at which one circular arc is + joined to another circular arc in the approximation to the VV CCL + + The calcuation of the theta factor found in Eq 4 of Y. Itoh et al. The + full form of the integral is given in appendix A. """ theta2 = np.pi / 2.0 + theta1_vv a = (Ro_vv - Ri_vv) / 2.0 From 04e936ac27c75fa14354bddcc6b8bb79a520d824 Mon Sep 17 00:00:00 2001 From: apearce Date: Fri, 22 Nov 2024 17:36:11 +0000 Subject: [PATCH 6/8] =?UTF-8?q?=E2=9C=85=20fix=20unit=20tests=20and=20fix?= =?UTF-8?q?=20bug=20in=20integral?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/sctfcoil.py | 9 ++++++--- source/fortran/tfcoil_variables.f90 | 2 +- tests/unit/test_sctfcoil.py | 4 ++-- 3 files changed, 9 insertions(+), 6 deletions(-) diff --git a/process/sctfcoil.py b/process/sctfcoil.py index b908368d..2d364af5 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -7266,7 +7266,7 @@ def _theta_factor_integral(Ro_vv, Ri_vv, Rm_vv, H_vv, theta1_vv): lambda_term(tau[1, k], omega[k]) - lambda_term(tau[0, k], omega[k]) ) - return (chi1 + chi2) / np.pi + return (chi1 + 2.0 * chi2) / (2.0 * np.pi) def vv_stress_on_quench( @@ -7350,8 +7350,11 @@ def vv_stress_on_quench( Empirical Formulas for Estimating Self and Mutual Inductances of Toroidal Field Coils and Structures. Plasma and Fusion Research. 15. 1405078-1405078. 10.1585/pfr.15.1405078. """ + # Convert angles into radians + theta1_vv_rad = np.pi * (theta1_vv / 180.0) + # Poloidal loop resistance (PLR) in ohms - theta_vv = _theta_factor_integral(Ro_vv, Ri_vv, Rm_vv, H_vv, theta1_vv) + theta_vv = _theta_factor_integral(Ro_vv, Ri_vv, Rm_vv, H_vv, theta1_vv_rad) plr_coil = ((0.5 * ccl_length_coil) / (n_tf * (S_cc + S_rp))) * 1e-6 plr_vv = ((0.84 / d_vv) * theta_vv) * 1e-6 @@ -7397,5 +7400,5 @@ def vv_stress_on_quench( zeta = 1 + ((A_vv - 1) * numpy.log((A_vv + 1) / (A_vv - 1)) / (2 * A_vv)) tresca_stress_vv = zeta * B_vvi * J_vvi * Ri_vv - + print(tresca_stress_vv) return tresca_stress_vv diff --git a/source/fortran/tfcoil_variables.f90 b/source/fortran/tfcoil_variables.f90 index f8b989ed..64c0f950 100644 --- a/source/fortran/tfcoil_variables.f90 +++ b/source/fortran/tfcoil_variables.f90 @@ -1054,7 +1054,7 @@ subroutine init_tfcoil_variables i_cp_joints = -1 cryo_cool_req = 0.0D0 theta1_coil = 45.0D0 - theta1_vv = 1.0D0 + theta1_vv = 1.0D0 ! 1 Deg max_vv_stress = 143.0D6 end subroutine init_tfcoil_variables end module tfcoil_variables diff --git a/tests/unit/test_sctfcoil.py b/tests/unit/test_sctfcoil.py index 0a6e8006..878f3d37 100644 --- a/tests/unit/test_sctfcoil.py +++ b/tests/unit/test_sctfcoil.py @@ -14031,7 +14031,7 @@ def test_vv_stress_on_quench(): d_vv=0.12, # for 6.6 restistance -> lambda2 = 2.1 ) ) - == 56835032.21809308 + == 57045874.69917925 ) @@ -14094,4 +14094,4 @@ def test_vv_stress_on_quench_integration(sctfcoil, monkeypatch): sctfcoil.vv_stress_on_quench() - assert pytest.approx(sctfcoil_module.vv_stress_quench) == 56834395.24352395 + assert pytest.approx(sctfcoil_module.vv_stress_quench) == 56893800.120420754 From 31f72d74d7eb084444deeaf4115422dea0117e24 Mon Sep 17 00:00:00 2001 From: apearce Date: Fri, 22 Nov 2024 17:37:21 +0000 Subject: [PATCH 7/8] =?UTF-8?q?=F0=9F=94=A5=20remove=20left=20over=20print?= =?UTF-8?q?=20statement?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- process/sctfcoil.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/process/sctfcoil.py b/process/sctfcoil.py index 2d364af5..9d5d3002 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -7400,5 +7400,5 @@ def vv_stress_on_quench( zeta = 1 + ((A_vv - 1) * numpy.log((A_vv + 1) / (A_vv - 1)) / (2 * A_vv)) tresca_stress_vv = zeta * B_vvi * J_vvi * Ri_vv - print(tresca_stress_vv) + return tresca_stress_vv From 4e15e5c7c586affa6e74d95715ad694dddf24cc9 Mon Sep 17 00:00:00 2001 From: apearce Date: Wed, 27 Nov 2024 15:22:26 +0000 Subject: [PATCH 8/8] =?UTF-8?q?=E2=9C=85=20fix=20unit=20test=20missing=20v?= =?UTF-8?q?ariable?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- tests/unit/test_sctfcoil.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/tests/unit/test_sctfcoil.py b/tests/unit/test_sctfcoil.py index 878f3d37..fa8d8e54 100644 --- a/tests/unit/test_sctfcoil.py +++ b/tests/unit/test_sctfcoil.py @@ -14059,8 +14059,9 @@ def test_vv_stress_on_quench_integration(sctfcoil, monkeypatch): monkeypatch.setattr(sctfcoil_module, "a_tf_steel", 0.55) # Section 3 # Sum from Section 3 - monkeypatch.setattr(sctfcoil_module, "a_case_front", 0.47) - monkeypatch.setattr(sctfcoil_module, "a_case_nose", 0.47) + monkeypatch.setattr(sctfcoil_module, "a_case_front", 0.42) + monkeypatch.setattr(sctfcoil_module, "a_case_nose", 0.42) + monkeypatch.setattr(sctfcoil_module, "t_lat_case_av", 0.05) monkeypatch.setattr(build_variables, "vgap_xpoint_divertor", 0.05) # Baseline 2018 monkeypatch.setattr(build_variables, "shldtth", 0.3) # Baseline 2018