diff --git a/process/physics.py b/process/physics.py index ac895a08cc..f72a1a8a21 100644 --- a/process/physics.py +++ b/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, diff --git a/process/pulse.py b/process/pulse.py index 561b2ce487..37737c8710 100755 --- a/process/pulse.py +++ b/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, ) diff --git a/source/fortran/physics_variables.f90 b/source/fortran/physics_variables.f90 index 995df178f9..120302b253 100644 --- a/source/fortran/physics_variables.f90 +++ b/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 diff --git a/source/fortran/times_variables.f90 b/source/fortran/times_variables.f90 index d08c094f9d..6933135af3 100644 --- a/source/fortran/times_variables.f90 +++ b/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 ', & diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 125569741a..f4a227da60 100644 --- a/tests/unit/test_physics.py +++ b/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,