Current ramp-up #489

process/physics.py

@@ -204,9 +204,8 @@ def rether(alphan, alphat, dene, dlamie, te, ti, zeffai):
     conie = 2.42165e-41 * dlamie * dene**2 * zeffai * profie
 
     return conie * (ti - te) / (te**1.5)
+
 
-
-@nb.jit(nopython=True, cache=True)
 def vscalc(
     csawth,
     eps,
@@ -248,10 +247,23 @@ def vscalc(
     rlpint = rmu0 * rmajor * rli / 2.0
     phiint = rlpint * plascur
 
-    # Start-up resistive component
-    # Uses ITER formula without the 10 V-s add-on
+    # Issue #489, Resistive flux consumption during the current ramp
+    # rho_cr = plasma loop resistance at the (assumed) conditions during the current ramp
+    # Typical temperature and Zeff during the current ramp:
+
+    t_current_ramp = physics_variables.t_cr_ratio * physics_variables.te
+
+    zeff_current_ramp = (
+        physics_variables.z_eff_cr_ratio * (physics_variables.zeff - 1) + 1
+    )
+
+    # Plasma loop resistance at the (assumed) conditions during the current ramp
+    _, rho_cr = plasma_resistance(
+        kappa, rmajor, physics_variables.rminor, t_current_ramp, zeff_current_ramp
+    )
 
-    vsres = gamma * rmu0 * plascur * rmajor
+    # Resistive flux consumption during the current ramp
+    vsres = 0.5 * kappa * gamma * plascur * times_variables.tohs * rho_cr
 
     # Hirshman, Neilson: Physics of Fluids, 29 (1986) p790
     # fit for external inductance
@@ -284,7 +296,7 @@ def vscalc(
     vsbrn = vburn * (t_fusion_ramp + tburn)
     vsstt = vsstt + vsbrn
 
-    return phiint, rlp, vsbrn, vsind, vsres, vsstt
+    return phiint, rlp, vsbrn, vsind, vsres, vsstt, rho_cr
 
 
 @nb.jit(nopython=True, cache=True)
@@ -768,6 +780,36 @@ def tpf(j, triang, sqeps, fit=1):
     raise RuntimeError(f"fit={fit} is not valid. Must be 1, 2, or 3")
 
 
+def plasma_resistance(kappa95, rmajor, rminor, ten, zeff):
+    # Density weighted electron temperature in 10 keV units
+
+    t10 = ten / 10.0
+    rplas = (
+        physics_variables.plasma_res_factor
+        * 2.15e-9
+        * zeff
+        * rmajor
+        / (kappa95 * rminor**2 * t10**1.5)
+    )
+
+    # Neo-classical resistivity enhancement factor
+    # Taken from  N. A. Uckan et al, Fusion Technology 13 (1988) p.411.
+    # The expression is valid for aspect ratios in the range 2.5--4.
+
+    rpfac = 4.3 - 0.6 * rmajor / rminor
+    rplas = rplas * rpfac
+
+    # Check to see if plasma resistance is negative
+    # (possible if aspect ratio is too high)
+
+    if rplas <= 0.0:
+        error_handling.fdiags[0] = rplas
+        error_handling.fdiags[1] = physics_variables.aspect
+        error_handling.report_error(83)
+
+    return rpfac, rplas
+
+
 class Physics:
     def __init__(self, plasma_profile, current_drive):
         self.outfile = constants.nout
@@ -1335,6 +1377,7 @@ def physics(self):
             physics_variables.vsind,
             physics_variables.vsres,
             physics_variables.vsstt,
+            physics_variables.rho_cr,
         ) = vscalc(
             physics_variables.csawth,
             physics_variables.eps,
@@ -1960,34 +2003,9 @@ def phyaux(
 
     @staticmethod
     def pohm(facoh, kappa95, plascur, rmajor, rminor, ten, vol, zeff):
-        # Density weighted electron temperature in 10 keV units
-
-        t10 = ten / 10.0
 
         # Plasma resistance, from loop voltage calculation in IPDG89
-
-        rplas = (
-            physics_variables.plasma_res_factor
-            * 2.15e-9
-            * zeff
-            * rmajor
-            / (kappa95 * rminor**2 * t10**1.5)
-        )
-
-        # Neo-classical resistivity enhancement factor
-        # Taken from  N. A. Uckan et al, Fusion Technology 13 (1988) p.411.
-        # The expression is valid for aspect ratios in the range 2.5--4.
-
-        rpfac = 4.3 - 0.6 * rmajor / rminor
-        rplas = rplas * rpfac
-
-        # Check to see if plasma resistance is negative
-        # (possible if aspect ratio is too high)
-
-        if rplas <= 0.0:
-            error_handling.fdiags[0] = rplas
-            error_handling.fdiags[1] = physics_variables.aspect
-            error_handling.report_error(83)
+        rpfac, rplas = plasma_resistance(kappa95, rmajor, rminor, ten, zeff)
 
         # Ohmic heating power per unit volume
         # Corrected from: pohmpv = (facoh*plascur)**2 * ...
@@ -4203,6 +4221,30 @@ def outplas(self):
                 physics_variables.csawth,
             )
 
+            po.ovarrf(
+                self.outfile,
+                "Temperature during the current ramp as a fraction of the flat-top value",
+                "(t_cr_ratio)",
+                physics_variables.t_cr_ratio,
+                "IP ",
+            )
+
+            po.ovarrf(
+                self.outfile,
+                "(Zeff-1) during the current ramp as a fraction of the flat-top value",
+                "(z_eff_cr_ratio)",
+                physics_variables.z_eff_cr_ratio,
+                "IP ",
+            )
+
+            po.ovarre(
+                self.outfile,
+                "Plasma loop resistance during the current ramp (ohm)",
+                "(rho_cr)",
+                physics_variables.rho_cr,
+                "OP ",
+            )
+
         po.osubhd(self.outfile, "Fuelling :")
         po.ovarre(
             self.outfile,

process/pulse.py

@@ -181,13 +181,13 @@ def burn(self, output: bool):
             )
             po.ovarre(
                 self.outfile,
-                "Required volt-seconds during start-up (Wb)",
+                "Required volt-seconds during current ramp-up (Wb)",
                 "(vssoft)",
                 vssoft,
             )
             po.ovarre(
                 self.outfile,
-                "Available volt-seconds during burn (Wb)",
+                "Available volt-seconds during fusion power ramp and flat-top (Wb)",
                 "(vsmax)",
                 vsmax,
             )

source/fortran/physics_variables.f90

@@ -230,6 +230,15 @@ module physics_variables
   real(dp) :: gamma
   !! Ejima coefficient for resistive startup V-s formula
 
+  real(dp) :: t_cr_ratio 
+  !! Typical temperature during the current ramp as a fraction of the flat-top value
+
+  real(dp) :: z_eff_cr_ratio 
+  !! Typical (Zeff-1) during the current ramp as a fraction of the flat-top value
+
+  real(dp) :: rho_cr 
+  !! Plasma loop resistance during the current ramp (ohm)
+
   real(dp) :: gammaft
   !! ratio of (fast alpha + neutral beam beta) to thermal beta
 
@@ -944,6 +953,8 @@ subroutine init_physics_variables
     fusionrate = 0.0D0
     fvsbrnni = 1.0D0
     gamma = 0.4D0
+    t_cr_ratio = 0.25D0
+    z_eff_cr_ratio = 0.25D0
     gammaft = 0.0D0
     hfac = 0.0D0
     hfact = 1.0D0

source/fortran/times_variables.f90

@@ -78,7 +78,9 @@ subroutine init_times_variables
     tcycle = 0.0D0
     tdown = 0.0D0
     tdwell = 1800.0D0
-    t_fusion_ramp = 10.0D0
+    ! Issue #489 Current ramp-up.  700 s taken from 
+    ! Impact of the plasma operation on the technical requirements in EU-DEMO, M. Siccinio
+    t_fusion_ramp = 700.0D0  
     tim = 0.0D0
     timelabel = (/ 'Start', &
       'BOP  ', &

tests/unit/test_physics.py

@@ -1344,8 +1344,8 @@ class VscalcParam(NamedTuple):
             expected_rlp=1.4075705307248088e-05,
             expected_vsbrn=42.109179697761263,
             expected_vsind=258.97124024420435,
-            expected_vsres=55.488435095110333,
-            expected_vsstt=356.56885503707593,
+            expected_vsres=3.8778020459507383,
+            expected_vsstt=304.9582219879163,
         ),
         VscalcParam(
             csawth=1,
@@ -1363,8 +1363,8 @@ class VscalcParam(NamedTuple):
             expected_rlp=1.4075705307248088e-05,
             expected_vsbrn=0.41692257126496302,
             expected_vsind=258.97124024420435,
-            expected_vsres=55.488435095110333,
-            expected_vsstt=314.87659791057968,
+            expected_vsres=3.8778020459507383,
+            expected_vsstt=263.26596486142006,
         ),
     ),
 )
@@ -1378,7 +1378,7 @@ def test_vscalc(vscalcparam):
     :type vscalcparam: vscalcparam
     """
 
-    phiint, rlp, vsbrn, vsind, vsres, vsstt = vscalc(
+    phiint, rlp, vsbrn, vsind, vsres, vsstt, rho_cr = vscalc(
         csawth=vscalcparam.csawth,
         eps=vscalcparam.eps,
         facoh=vscalcparam.facoh,