diff --git a/documentation/proc-pages/physics-models/plasma_current/plasma_current.md b/documentation/proc-pages/physics-models/plasma_current/plasma_current.md index feccb438b7..c2711f7d07 100644 --- a/documentation/proc-pages/physics-models/plasma_current/plasma_current.md +++ b/documentation/proc-pages/physics-models/plasma_current/plasma_current.md @@ -511,7 +511,7 @@ It is not possible to derive a general analytic expression for $F$, so it has be The cylindrical equivalent $q$ definition used currently is the one given below from IPDG89[^5] $$ -q_* = \frac{5 \times 10^6a^2B_T}{RI_{\text{p}}}\frac{(1+\kappa^2(1+2\delta^2-1.2\delta^3))}{2} +q_* = \frac{5 \times 10^6a^2B_T}{RI_{\text{p}}}\frac{(1+\kappa_{95}^2(1+2\delta_{95}^2-1.2\delta_{95}^3))}{2} $$ diff --git a/process/hcpb.py b/process/hcpb.py index 300e6cc2ed..b04b81a149 100644 --- a/process/hcpb.py +++ b/process/hcpb.py @@ -765,7 +765,12 @@ def powerflow_calc(self, output: bool): + fwbs_variables.psurffwo + fwbs_variables.pnucblkt ) - primary_pumping_variables.htpmw_fw_blkt = fpump / (1 - fpump) * p_plasma + primary_pumping_variables.htpmw_fw_blkt = ( + primary_pumping_variables.f_p_fw_blkt_pump + * fpump + / (1 - fpump) + * p_plasma + ) # For divertor and shield, mechanical pumping power is a fraction of thermal # power removed by coolant @@ -825,6 +830,13 @@ def powerflow_calc(self, output: bool): primary_pumping_variables.htpmw_fw_blkt, "OP ", ) + po.ovarre( + self.outfile, + "Pumping power for FW and Blanket multiplier factor", + "(f_p_fw_blkt_pump)", + primary_pumping_variables.f_p_fw_blkt_pump, + "IP ", + ) po.ovarre( self.outfile, "Mechanical pumping power for divertor (MW)", diff --git a/process/physics.py b/process/physics.py index 47f59606ea..7bca265167 100644 --- a/process/physics.py +++ b/process/physics.py @@ -3343,7 +3343,7 @@ def calculate_plasma_current( * rminor**2 / (rmajor * plasma_current / bt) * 0.5 - * (1.0 + kappa**2 * (1.0 + 2.0 * triang**2 - 1.2 * triang**3)) + * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) ) # Normalised beta from Troyon beta limit diff --git a/process/sctfcoil.py b/process/sctfcoil.py index 4624d1474e..273fae1944 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -2367,13 +2367,14 @@ 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=np.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, # VV properties - d_vv=build_variables.dr_vv_inboard, + d_vv=build_variables.dr_vv_shells, ) def tf_field_and_force(self): @@ -7199,7 +7200,9 @@ def eyoung_parallel_array(n, eyoung_j_in, a_in, poisson_j_perp_in): return eyoung_j_out, a_out, poisson_j_perp_out -def _inductance_factor(H, ri, ro, rm, theta1): +def _inductance_factor( + H: float, ri: float, ro: float, rm: float, theta1: float +) -> float: """Calculates the inductance factor for a toroidal structure using surrogate 2 #1866. @@ -7232,30 +7235,112 @@ def _inductance_factor(H, ri, ro, rm, theta1): ) +@staticmethod +def lambda_term(tau: float, omega: float) -> float: + """ + The lambda function used inegral in inductance calcuation found + in Y. Itoh et al. The full form of the functions are given in appendix A. + + 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 + :returns: integral + """ + 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 + + +@staticmethod +def _theta_factor_integral( + ro_vv: float, ri_vv: float, rm_vv: float, h_vv: float, theta1_vv: float +) -> float: + """ + 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. + + 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 + :returns: theta factor. + """ + theta2 = np.pi / 2.0 + theta1_vv + a = (ro_vv - ri_vv) / 2.0 + rbar = (ro_vv + ri_vv) / 2.0 + 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) * ((rbar / 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 + 2.0 * chi2) / (2.0 * np.pi) + + def vv_stress_on_quench( # TF shape - H_coil, - ri_coil, - ro_coil, - rm_coil, - ccl_length_coil, - theta1_coil, + H_coil: float, + ri_coil: float, + ro_coil: float, + rm_coil: float, + ccl_length_coil: float, + theta1_coil: float, # VV shape - H_vv, - ri_vv, - ro_vv, - rm_vv, - theta1_vv, + H_vv: float, + ri_vv: float, + ro_vv: float, + rm_vv: float, + theta1_vv: float, # TF properties - n_tf_coils, - n_tf_turn, - s_rp, - s_cc, - taud, - i_op, + n_tf_coils: float, + n_tf_turn: float, + s_rp: float, + s_cc: float, + taud: float, + i_op: float, # VV properties - d_vv, -): + d_vv: float, +) -> float: """Generic model to calculate the Tresca stress of the Vacuum Vessel (VV), experienced when the TF coil quenches. @@ -7295,7 +7380,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 @@ -7304,8 +7389,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 ---------- @@ -7313,9 +7398,13 @@ 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_rad) plr_coil = ((0.5 * ccl_length_coil) / (n_tf_coils * (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 = ( diff --git a/source/fortran/build_variables.f90 b/source/fortran/build_variables.f90 index a71f1f2048..3c8bc900ba 100644 --- a/source/fortran/build_variables.f90 +++ b/source/fortran/build_variables.f90 @@ -78,6 +78,9 @@ module build_variables real(dp) :: dz_vv_lower !! vacuum vessel underside thickness (TF coil / shield) (m) + real(dp) :: dr_vv_shells + !! vacuum vessel double walled shell thicknesses (m) + real(dp) :: f_avspace !! F-value for stellarator radial space check (`constraint equation 83`) @@ -327,6 +330,7 @@ subroutine init_build_variables dr_vv_outboard = 0.07D0 dz_vv_upper = 0.07D0 dz_vv_lower = 0.07D0 + dr_vv_shells = 0.12D0 f_avspace = 1.0D0 fcspc = 0.6D0 fseppc = 3.5D8 diff --git a/source/fortran/constraint_equations.f90 b/source/fortran/constraint_equations.f90 index 9e61ee3dce..9ddb721eb5 100755 --- a/source/fortran/constraint_equations.f90 +++ b/source/fortran/constraint_equations.f90 @@ -2640,7 +2640,9 @@ subroutine constraint_eqn_068(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) !! q95 : input real : safety factor q at 95% flux surface !! aspect : input real : aspect ratio (iteration variable 1) !! rmajor : input real : plasma major radius (m) (iteration variable 3) - use constraint_variables, only: fpsepbqar, psepbqarmax + !! i_q95_fixed : input int : Switch that allows for fixing q95 only in this constraint. + !! q95_fixed : input real : fixed safety factor q at 95% flux surface + use constraint_variables, only: fpsepbqar, psepbqarmax, i_q95_fixed, q95_fixed use physics_variables, only: pdivt, bt, q95, aspect, rmajor implicit none real(dp), intent(out) :: tmp_cc @@ -2649,9 +2651,14 @@ subroutine constraint_eqn_068(tmp_cc, tmp_con, tmp_err, tmp_symbol, tmp_units) character(len=1), intent(out) :: tmp_symbol character(len=10), intent(out) :: tmp_units - tmp_cc = 1.0d0 - fpsepbqar * psepbqarmax / ((pdivt*bt)/(q95*aspect*rmajor)) + if (i_q95_fixed == 1) then + tmp_cc = 1.0d0 - fpsepbqar * psepbqarmax / ((pdivt*bt)/(q95_fixed*aspect*rmajor)) + tmp_err = (pdivt*bt)/(q95_fixed*aspect*rmajor) - psepbqarmax + else + tmp_cc = 1.0d0 - fpsepbqar * psepbqarmax / ((pdivt*bt)/(q95*aspect*rmajor)) + tmp_err = (pdivt*bt)/(q95*aspect*rmajor) - psepbqarmax + end if tmp_con = psepbqarmax - tmp_err = (pdivt*bt)/(q95*aspect*rmajor) - psepbqarmax tmp_symbol = '<' tmp_units = 'MWT/m' diff --git a/source/fortran/constraint_variables.f90 b/source/fortran/constraint_variables.f90 index 753415e4c4..f96f696dea 100644 --- a/source/fortran/constraint_variables.f90 +++ b/source/fortran/constraint_variables.f90 @@ -86,6 +86,10 @@ module constraint_variables !! constraint equation icc = 46 !! iteration variable ixc = 72 + real(dp) :: q95_fixed + !! fixed safety factor q at 95% flux surface + !! (`constraint equation 68`) + real(dp) :: fjohc !! f-value for central solenoid current at end-of-flattop !! (`constraint equation 26`, `iteration variable 38`) @@ -221,6 +225,10 @@ module constraint_variables real(dp) :: gammax !! maximum current drive gamma (`constraint equation 37`) + integer :: i_q95_fixed + !! Switch that allows for fixing q95 only in this constraint equation 68. + !! (`constraint equation 68`) + real(dp) :: pflux_fw_rad_max !! Maximum permitted radiation wall load (MW/m^2) (`constraint equation 67`) @@ -332,6 +340,7 @@ subroutine init_constraint_variables fhldiv = 1.0D0 fiooic = 0.5D0 fipir = 1.0D0 + q95_fixed = 3.0D0 fjohc = 1.0D0 fjohc0 = 1.0D0 fjprot = 1.0D0 @@ -372,6 +381,7 @@ subroutine init_constraint_variables fwalld = 1.0D0 fzeffmax = 1.0D0 gammax = 2.0D0 + i_q95_fixed = 0 pflux_fw_rad_max = 1.0D0 mvalim = 40.0D0 nbshinefmax = 1.0D-3 diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index e13b9df909..4bc55f5a33 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -177,7 +177,8 @@ subroutine parse_input_file(in_file,out_file,show_changes) fcqt, fzeffmax, fstrcase, fhldiv, foh_stress, fwalld, gammax, fjprot, & ft_current_ramp_up, tcycmn, auxmin, zeffmax, f_fw_rad_max, fdtmp, fpoloidalpower, & fnbshinef, freinke, fvvhe, fqval, fq, fmaxvvstress, fbeta_poloidal, fbeta_poloidal_eps, fjohc, & - fflutf, bmxlim, t_burn_min, fbeta_min, fecrh_ignition, fstr_wp, fncycle + fflutf, bmxlim, t_burn_min, fbeta_min, fecrh_ignition, fstr_wp, fncycle, i_q95_fixed, & + q95_fixed use cost_variables, only: ucich, uctfsw, dintrt, ucblbe, uubop, dtlife, & cost_factor_vv, cfind, uccry, fcap0cp, uccase, uuves, cconshtf, conf_mag, & ucbllipb, ucfuel, uumag, ucpfbs, ireactor, uucd, div_umain_time, div_nu, & @@ -257,7 +258,7 @@ subroutine parse_input_file(in_file,out_file,show_changes) use pulse_variables, only: i_pulsed_plant, dtstor, itcycl, istore, bctmp use primary_pumping_variables, only: t_in_bb, t_out_bb, dp_he, p_he, gamma_he, & - dp_fw_blkt, dp_fw, dp_blkt, dp_liq + dp_fw_blkt, dp_fw, dp_blkt, dp_liq, f_p_fw_blkt_pump use scan_module, only: isweep_2, nsweep, isweep, scan_dim, nsweep_2, & sweep_2, sweep, ipnscns, ipnscnv @@ -770,6 +771,9 @@ subroutine parse_input_file(in_file,out_file,show_changes) case ('i_hldiv') call parse_int_variable('i_hldiv', i_hldiv, 0, 2, & 'Switch for user input hldiv') + case ('i_q95_fixed') + call parse_int_variable('i_q95_fixed', i_q95_fixed, 0, 1, & + 'Switch that allows fixed q95 within PsepB/qAR (68 constraint)') case ('fflutf') call parse_real_variable('fflutf', fflutf, 0.001D0, 10.0D0, & 'F-value for neutron fluence on TF coil') @@ -779,6 +783,9 @@ subroutine parse_input_file(in_file,out_file,show_changes) case ('fiooic') call parse_real_variable('fiooic', fiooic, 0.001D0, 10.0D0, & 'F-value for SCTF iop/icrit') + case ('q95_fixed') + call parse_real_variable('q95_fixed', q95_fixed, 0.001D0, 10.0D0, & + 'fixed safety factor q at 95% flux surface') case ('fjprot') call parse_real_variable('fjprot', fjprot, 0.001D0, 10.0D0, & 'F-value for SCTF winding pack J') @@ -2274,7 +2281,9 @@ subroutine parse_input_file(in_file,out_file,show_changes) case ('p_he') call parse_real_variable('p_he', p_he, 0.0D0, 100.0D6, & 'Pressure in FW and blanket coolant at pump exit') - + case ('f_p_fw_blkt_pump') + call parse_real_variable('f_p_fw_blkt_pump', f_p_fw_blkt_pump, 0.0D0, 10.0D0, & + 'Pumping power for FW and Blanket multiplier factor') case ('gamma_he') call parse_real_variable('gamma_he', gamma_he, 1.0D0, 2.0D0, & 'Ratio of specific heats for helium or for another Gas') diff --git a/source/fortran/primary_pumping_variables.f90 b/source/fortran/primary_pumping_variables.f90 index 9ca8823d0c..c06a72f8d2 100644 --- a/source/fortran/primary_pumping_variables.f90 +++ b/source/fortran/primary_pumping_variables.f90 @@ -50,6 +50,9 @@ module primary_pumping_variables !! mechanical pumping power for FW and blanket including heat exchanger and !! pipes (`primary_pumping=3`) [MW] + real(dp) :: f_p_fw_blkt_pump + !! Pumping power for FW and Blanket multiplier factor + contains subroutine init_primary_pumping_variables @@ -68,6 +71,7 @@ subroutine init_primary_pumping_variables dp_blkt = 3.5D3 ! Pa dp_liq = 1.0D7 ! Pa htpmw_fw_blkt = 0.0D0 + f_p_fw_blkt_pump = 1.0D0 end subroutine init_primary_pumping_variables diff --git a/source/fortran/tfcoil_variables.f90 b/source/fortran/tfcoil_variables.f90 index 26bb0e5969..c91793faa2 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_physics.py b/tests/unit/test_physics.py index f44d667a39..63c3a46403 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -947,7 +947,7 @@ class PlasmaCurrentParam(NamedTuple): expected_alphaj=1.9008029008029004, expected_ind_plasma_internal_norm=1.2064840230894305, expected_bp=0.96008591022564971, - expected_qstar=3.869423496255382, + expected_qstar=2.900802902105021, expected_plasma_current=18398455.678867526, ), PlasmaCurrentParam( @@ -974,7 +974,7 @@ class PlasmaCurrentParam(NamedTuple): expected_alphaj=1.9008029008029004, expected_ind_plasma_internal_norm=1.2064840230894305, expected_bp=0.96008591022564971, - expected_qstar=3.869423496255382, + expected_qstar=2.900802902105021, expected_plasma_current=18398455.678867526, ), ), diff --git a/tests/unit/test_sctfcoil.py b/tests/unit/test_sctfcoil.py index a92d0476da..92c1ebd80b 100644 --- a/tests/unit/test_sctfcoil.py +++ b/tests/unit/test_sctfcoil.py @@ -14076,7 +14076,7 @@ def test_vv_stress_on_quench(): d_vv=0.12, # for 6.6 restistance -> lambda2 = 2.1 ) ) - == 56835032.21809308 + == 57045874.69917925 ) @@ -14105,8 +14105,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, "dz_xpoint_divertor", 0.05) # Baseline 2018 monkeypatch.setattr(build_variables, "dz_shld_upper", 0.3) # Baseline 2018 @@ -14142,4 +14143,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