diff --git a/process/physics.py b/process/physics.py index 1de78bf5f3..0357a2c5fc 100644 --- a/process/physics.py +++ b/process/physics.py @@ -2015,7 +2015,7 @@ def culcur( p0, pperim, q0, - qpsi, + q95, rli, rmajor, rminor, @@ -2037,6 +2037,8 @@ def culcur( 5 = Todd empirical scaling I 6 = Todd empirical scaling II 7 = Connor-Hastie model + 8 = Sauter scaling (allowing negative triangularity) Issue #392 + 'Geometric formulas for system codes including the effect of negative triangularity' iprofile : input integer : switch for current profile consistency 0 = use input alphaj, rli 1 = make these consistent with q, q0 @@ -2045,7 +2047,7 @@ def culcur( p0 : input real : central plasma pressure (Pa) pperim : input real : plasma perimeter length (m) q0 : input real : plasma safety factor on axis - qpsi : input real : plasma edge safety factor (= q-bar for icurr=2) + q95 : input real : plasma safety factor at 95% flux (= q-bar for icurr=2) rli : input/output real : plasma normalised internal inductance rmajor : input real : major radius (m) rminor : input real : minor radius (m) @@ -2055,11 +2057,9 @@ def culcur( bp : output real : poloidal field (T) qstar : output real : equivalent cylindrical safety factor (shaped) plascur : output real : plasma current (A) - This routine performs the calculation of the - plasma current, with a choice of formula for the edge - safety factor. It will also make the current profile parameters + This routine calculates the plasma current based on the edge + safety factor q95. It will also make the current profile parameters consistent with the q-profile if required. - AEA FUS 251: A User's Guide to the PROCESS Systems Code J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, unpublished internal Oak Ridge document Y.-K. M. Peng, J. Galambos and P.C. Shipe, 1992, @@ -2082,7 +2082,7 @@ def culcur( # Calculate the function Fq that scales the edge q from the # circular cross-section cylindrical case - # First check for negative triangularity using unsuitable current scaling + # Only the Sauter scaling (icurr=8) is suitable for negative triangularity: if icurr != 8 and triang < 0.0: raise ValueError( @@ -2092,7 +2092,7 @@ def culcur( if icurr == 1: # Peng analytical fit fq = (1.22 - 0.68 * eps) / ((1.0 - eps * eps) ** 2) * sf**2 elif icurr == 2: # Peng scaling for double null divertor; TARTs [STAR Code] - curhat = 1.0e6 * plasc(qpsi, asp, eps, rminor, bt, kappa, triang) / bt + plascur = 1.0e6 * plasc(q95, asp, eps, rminor, bt, kappa, triang) elif icurr == 3: # Simple ITER scaling (simply the cylindrical case) fq = 1.0 elif icurr == 4: # ITER formula (IPDG89) @@ -2130,7 +2130,8 @@ def culcur( w07 = 1.0 # zero squareness - can be modified later if required fq = ( - (1.0 + 1.2 * (kappa - 1.0) + 0.56 * (kappa - 1.0) ** 2) + (4.1e6 / 5.0e6) + * (1.0 + 1.2 * (kappa - 1.0) + 0.56 * (kappa - 1.0) ** 2) * (1.0 + 0.09 * triang + 0.16 * triang**2) * (1.0 + 0.45 * triang * eps) / (1.0 - 0.74 * eps) @@ -2141,26 +2142,23 @@ def culcur( else: raise ValueError(f"Invalid value {icurr=}") - if icurr == 8: - curhat = 4.1e6 * rminor**2 / (rmajor * qpsi) * fq - elif icurr != 2: - curhat = 5.0e6 * rminor**2 / (rmajor * qpsi) * fq + if icurr != 2: + plascur = 5.0e6 * rminor**2 / (rmajor * q95) * fq * bt # == 2 case covered above qstar = ( 5.0e6 * rminor**2 - / (rmajor * curhat) + / (rmajor * plascur / bt) * 0.5 * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) ) - plascur = curhat * bt physics_variables.normalised_total_beta = ( 1.0e8 * physics_variables.beta * rminor * bt / plascur ) bp = bpol( - icurr, plascur, qpsi, asp, eps, bt, kappa, triang, pperim, constants.rmu0 + icurr, plascur, q95, asp, eps, bt, kappa, triang, pperim, constants.rmu0 ) # Ensure current profile consistency, if required diff --git a/process/plasma_geometry.py b/process/plasma_geometry.py index 31e3ea9042..0950fdd129 100644 --- a/process/plasma_geometry.py +++ b/process/plasma_geometry.py @@ -222,24 +222,8 @@ def geomty(self): physics_variables.kappa, physics_variables.triang, ) - - # Poloidal perimeter - physics_variables.pperim = 2.0e0 * (xo * thetao + xi * thetai) - physics_variables.sf = physics_variables.pperim / ( - 2.0e0 * numpy.pi * physics_variables.rminor - ) - - # Volume - physics_variables.vol = physics_variables.cvol * self.xvol( - physics_variables.rmajor, - physics_variables.rminor, - xi, - thetai, - xo, - thetao, - ) - - # Surface area + # Surface area - inboard and outboard. These are not given by Sauter but + # the outboard area is required by DCLL and divertor xsi, xso = self.xsurf( physics_variables.rmajor, physics_variables.rminor, @@ -249,10 +233,45 @@ def geomty(self): thetao, ) physics_variables.sareao = xso - physics_variables.sarea = xsi + xso - # Cross-sectional area - physics_variables.xarea = self.xsecta(xi, thetai, xo, thetao) + # icurr = 8 specifies use of the Sauter geometry as well as plasma current. + if physics_variables.icurr == 8: + ( + physics_variables.pperim, + physics_variables.sf, + physics_variables.sarea, + physics_variables.xarea, + physics_variables.vol, + ) = self.Sauter_geometry( + physics_variables.rminor, + physics_variables.rmajor, + physics_variables.kappa, + physics_variables.triang, + ) + + else: + + # Poloidal perimeter + physics_variables.pperim = 2.0e0 * (xo * thetao + xi * thetai) + physics_variables.sf = physics_variables.pperim / ( + 2.0e0 * numpy.pi * physics_variables.rminor + ) + + # Volume + physics_variables.vol = physics_variables.cvol * self.xvol( + physics_variables.rmajor, + physics_variables.rminor, + xi, + thetai, + xo, + thetao, + ) + + # Cross-sectional area + physics_variables.xarea = self.xsecta(xi, thetai, xo, thetao) + + # Surface area - sum of inboard and outboard. + physics_variables.sarea = xsi + xso def xparam(self, a, kap, tri): """ @@ -534,3 +553,40 @@ def xsect0(self, a, kap, tri): xsect0 = xlo**2 * (thetao - cto * sto) + xli**2 * (thetai - cti * sti) return xsect0 + + def sauter_geometry(self, a, r0, kap, tri): + """ + Plasma geometry based on equations (36) in O. Sauter, Fusion Engineering and Design 112 (2016) 633–645 + 'Geometric formulas for system codes including the effect of negative triangularity' + Author: Michael Kovari, issue #392 + a : input real : plasma minor radius (m) + r0 : input real : plasma major radius (m) + kap : input real : plasma separatrix elongation + tri : input real : plasma separatrix triangularity + """ + w07 = 1 + eps = a / r0 + + # Poloidal perimeter (named Lp in Sauter) + pperim = ( + 2.0e0 + * numpy.pi + * a + * (1 + 0.55 * (kap - 1)) + * (1 + 0.08 * tri**2) + * (1 + 0.2 * (w07 - 1)) + ) + + # A geometric factor + sf = pperim / (2.0e0 * numpy.pi * a) + + # Surface area (named Ap in Sauter) + sarea = 2.0e0 * numpy.pi * r0 * (1 - 0.32 * tri * eps) * pperim + + # Cross-section area (named S_phi in Sauter) + xarea = numpy.pi * a**2 * kap * (1 + 0.52 * (w07 - 1)) + + # Volume + vol = 2.0e0 * numpy.pi * r0 * (1 - 0.25 * tri * eps) * xarea + + return pperim, sf, sarea, xarea, vol diff --git a/process/sctfcoil.py b/process/sctfcoil.py index eb57d9f116..2a6ef7d435 100644 --- a/process/sctfcoil.py +++ b/process/sctfcoil.py @@ -804,7 +804,7 @@ def supercon( jtol = 1.0e4 for lap in range(100): - if ttest <= 0.0e0: + if ttest <= delt: error_handling.idiags[0] = lap error_handling.fdiags[0] = ttest error_handling.report_error(157) diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index d569b594c3..bf27e99fec 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -498,7 +498,7 @@ class CulcurParam(NamedTuple): q0: Any = None - qpsi: Any = None + q95: Any = None rmajor: Any = None @@ -541,7 +541,7 @@ class CulcurParam(NamedTuple): p0=0, pperim=24.081367139525412, q0=1, - qpsi=3.5, + q95=3.5, rmajor=8, rminor=2.6666666666666665, sf=1.4372507312498271, @@ -569,7 +569,7 @@ class CulcurParam(NamedTuple): p0=626431.90482713911, pperim=24.081367139525412, q0=1, - qpsi=3.5, + q95=3.5, rmajor=8, rminor=2.6666666666666665, sf=1.4372507312498271, @@ -616,7 +616,7 @@ def test_culcur(culcurparam, monkeypatch, physics): p0=culcurparam.p0, pperim=culcurparam.pperim, q0=culcurparam.q0, - qpsi=culcurparam.qpsi, + q95=culcurparam.q95, rmajor=culcurparam.rmajor, rminor=culcurparam.rminor, sf=culcurparam.sf,