VV stress under quench model review

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}
 $$
 
 

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)",

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

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,18 +7389,22 @@ 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
     ----------
     1. ITOH, Yasuyuki & Utoh, Hiroyasu & SAKAMOTO, Yoshiteru & Hiwatari, Ryoji. (2020).
     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 = (

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

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'
 

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

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')

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