diff --git a/documentation/proc-pages/eng-models/heating-and-current-drive.md b/documentation/proc-pages/eng-models/heating-and-current-drive.md index 26effb3a43..cea6e5c553 100644 --- a/documentation/proc-pages/eng-models/heating-and-current-drive.md +++ b/documentation/proc-pages/eng-models/heating-and-current-drive.md @@ -17,7 +17,8 @@ The fraction of the required plasma current to be produced by non-inductive mean - `iefrf` = 9: Oscillating Field current drive (RFPs only - OBSOLETE-REMOVED), - `iefrf` = 10: ECRH user input gamma, - `iefrf` = 11: ECRH "HARE" model [^3], -- `iefrf` = 12: EBW user scaling input. Scaling (S. Freethy) +- `iefrf` = 12: EBW user scaling input[^4], +- `iefrf` = 13: ECRH O-mode cutoff with Zeff and Te [^4], (Note that, at present, the neutral beam models do not include the effect of an edge transport barrier (pedestal) in the plasma profile.) @@ -63,4 +64,6 @@ Switch `ignite` can be used to denote whether the plasma is ignited, i.e. fully [^2]: T. C. Hender, M. K. Bevir, M. Cox, R. J. Hastie, P. J. Knight, C. N. Lashmore-Davies, B. Lloyd, G. P. Maddison, A. W. Morris, M. R. O'Brien, M.F. Turner abd H. R. Wilson, *"Physics Assessment for the European Reactor Study"*, AEA Fusion Report AEA FUS 172 (1992) -[^3]: E. Poli, M. Müller, H. Zohm, M. Kovari, *"Fast evaluation of the current driven by electron cyclotron waves for reactor studies"*, Physics of Plasmas 1 December 2018; 25 (12): 122501 \ No newline at end of file +[^3]: E. Poli, M. Müller, H. Zohm, M. Kovari, *"Fast evaluation of the current driven by electron cyclotron waves for reactor studies"*, Physics of Plasmas 1 December 2018; 25 (12): 122501 + +[^4]: Laqua, H & Maassberg, H & Marushchenko, Nikolai & Volpe, Francesco & Weller, A & Kasparek, W. (2003). *"Electron-Bernstein-Wave Current Drive in an Overdense Plasma at the Wendelstein 7-AS Stellarator"*, Physical review letters. 90. 075003. 10.1103/PhysRevLett.90.075003. \ No newline at end of file diff --git a/process/current_drive.py b/process/current_drive.py index 34958c0be4..ff897f43bf 100644 --- a/process/current_drive.py +++ b/process/current_drive.py @@ -194,8 +194,68 @@ def cudriv(self, output: bool): ) effcdfix = effcdfix * density_factor - effrfssfix = effrfssfix * density_factor + elif current_drive_variables.iefrffix == 13: + # ECCD model for O-mode cut-off with added Te and Zeff dependance + # Scaling author: Simon Freethy + # Ref : PROCESS issue #2994 + + fc = ( + 1 + / (2 * np.pi) + * constants.echarge + * physics_variables.bt + / constants.emass + ) + fp = ( + 1 + / (2 * np.pi) + * np.sqrt( + ( + (physics_variables.dene / 1.0e19) + * constants.echarge**2 + / (constants.emass * constants.epsilon0) + ) + ) + ) + + xi_CD = 0.18e0 # This is tuned to the results of a GRAY study + xi_CD = xi_CD * ( + 4.8e0 / (2 + physics_variables.zeff) + ) # Zeff correction + effrfssfix = ( + xi_CD + * physics_variables.te + / ( + 3.27e0 + * physics_variables.rmajor + * (physics_variables.dene / 1.0e19) + ) + ) + + # O-mode case + if current_drive_variables.wave_mode == 0: + f_cutoff = fp + + # X-mode case + elif current_drive_variables.wave_mode == 1: + f_cutoff = 0.5 * ( + fc + + np.sqrt( + current_drive_variables.harnum * fc**2 + 4 * fp**2 + ) + ) + + # Plasma coupling only occurs if the plasma cut-off is below the cyclotron harmonic + a = 0.1 # This controls how sharply the transition is reached + cutoff_factor = 0.5 * ( + 1 + + np.tanh( + (2 / (a)) + * ((current_drive_variables.harnum * fc - f_cutoff) / fp - a) + ) + ) + effcdfix = effrfssfix * cutoff_factor elif current_drive_variables.iefrffix != 0: raise RuntimeError( f"Current drive switch is invalid: {current_drive_variables.iefrffix = }" @@ -226,7 +286,7 @@ def cudriv(self, output: bool): * 1.0e6 ) faccdfix = auxiliary_cdfix / physics_variables.plascur - elif current_drive_variables.iefrffix in [3, 7, 10, 11, 12]: + elif current_drive_variables.iefrffix in [3, 7, 10, 12, 13]: # Injected power pinjemwfix = current_drive_variables.pinjfixmw @@ -403,11 +463,100 @@ def cudriv(self, output: bool): dene20 * physics_variables.rmajor ) current_drive_variables.effcd = effrfss - # EBWs can only couple to plasma if cyclotron harmonic is above plasma density cut-off; # this behaviour is captured in the following function (ref issue #1262): # current_drive_variables.harnum = cyclotron harmonic number (fundamental used as default) # contant 'a' controls sharpness of transition + a = 0.1e0 + + fc = ( + 1.0e0 + / (2.0e0 * np.pi) + * current_drive_variables.harnum + * constants.echarge + * physics_variables.bt + / constants.emass + ) + fp = ( + 1.0e0 + / (2.0e0 * np.pi) + * np.sqrt( + physics_variables.dene + * constants.echarge**2 + / (constants.emass * constants.epsilon0) + ) + ) + + density_factor = 0.5e0 * ( + 1.0e0 + np.tanh((2.0e0 / a) * ((fp - fc) / fp - a)) + ) + + current_drive_variables.effcd = ( + current_drive_variables.effcd * density_factor + ) + effrfss = effrfss * density_factor + + elif current_drive_variables.iefrf == 13: + # ECCD model for O-mode cut-off with added Te and Zeff dependance + # Scaling author: Simon Freethy + # Ref : PROCESS issue #2994 + + fc = ( + 1 + / (2 * np.pi) + * constants.echarge + * physics_variables.bt + / constants.emass + ) + fp = ( + 1 + / (2 * np.pi) + * np.sqrt( + ( + (physics_variables.dene / 1.0e19) + * constants.echarge**2 + / (constants.emass * constants.epsilon0) + ) + ) + ) + + xi_CD = 0.18e0 # This is tuned to the results of a GRAY study + xi_CD = xi_CD * ( + 4.8e0 / (2 + physics_variables.zeff) + ) # Zeff correction + effrfss = ( + xi_CD + * physics_variables.te + / ( + 3.27e0 + * physics_variables.rmajor + * (physics_variables.dene / 1.0e19) + ) + ) + + # O-mode case + if current_drive_variables.wave_mode == 0: + f_cutoff = fp + + # X-mode case + elif current_drive_variables.wave_mode == 1: + f_cutoff = 0.5 * ( + fc + + np.sqrt( + current_drive_variables.harnum * fc**2 + 4 * fp**2 + ) + ) + + # Plasma coupling only occurs if the plasma cut-off is below the cyclotron harmonic + a = 0.1 # This controls how sharply the transition is reached + cutoff_factor = 0.5 * ( + 1 + + np.tanh( + (2 / (a)) + * ((current_drive_variables.harnum * fc - f_cutoff) / fp - a) + ) + ) + current_drive_variables.effcd = effrfss * cutoff_factor else: raise RuntimeError( f"Current drive switch is invalid: {current_drive_variables.iefrf = }" @@ -442,7 +591,7 @@ def cudriv(self, output: bool): gamrf = effrfss * (dene20 * physics_variables.rmajor) current_drive_variables.gamcd = gamrf # ECCD - elif current_drive_variables.iefrf in [3, 7, 10, 11, 12]: + elif current_drive_variables.iefrf in [3, 7, 10, 12, 13]: # Injected power (set to close to close the Steady-state current equilibrium) current_drive_variables.echpwr = ( 1.0e-6 @@ -588,6 +737,11 @@ def cudriv(self, output: bool): ) elif current_drive_variables.iefrf == 12: po.ocmmnt(self.outfile, "EBW current drive") + elif current_drive_variables.iefrf == 13: + po.ocmmnt( + self.outfile, + "Electron Cyclotron Current Drive (O-mode cutoff with Zeff & Te)", + ) po.ovarin( self.outfile, @@ -610,6 +764,11 @@ def cudriv(self, output: bool): ) elif current_drive_variables.iefrffix == 12: po.ocmmnt(self.outfile, "EBW current drive") + elif current_drive_variables.iefrffix == 13: + po.ocmmnt( + self.outfile, + "Electron Cyclotron Current Drive (O-mode cutoff with Zeff & Te)", + ) po.ovarin( self.outfile, @@ -721,13 +880,27 @@ def cudriv(self, output: bool): "(gamma_ecrh)", current_drive_variables.gamma_ecrh, ) - elif current_drive_variables.iefrf == 12: + if current_drive_variables.iefrf == 12: po.ovarre( self.outfile, "EBW plasma heating efficiency", "(xi_ebw)", current_drive_variables.xi_ebw, ) + if current_drive_variables.iefrf in [12, 13]: + po.ovarre( + self.outfile, + "EC harmonic number", + "(harnum)", + current_drive_variables.harnum, + ) + if current_drive_variables.iefrf == 13: + po.ovarin( + self.outfile, + "EC cutoff wave mode switch", + "(wave_mode)", + current_drive_variables.wave_mode, + ) if current_drive_variables.iefrffix != 0: po.ovarre( diff --git a/process/io/plot_proc.py b/process/io/plot_proc.py index 06d07af343..9610947753 100755 --- a/process/io/plot_proc.py +++ b/process/io/plot_proc.py @@ -2283,7 +2283,7 @@ def plot_current_drive_info(axis, mfile_data, scan): if (iefrf == 5) or (iefrf == 8): nbi = True axis.text(-0.05, 1, "Neutral Beam Current Drive:", ha="left", va="center") - if (iefrf == 3) or (iefrf == 7) or (iefrf == 10) or (iefrf == 11): + if (iefrf == 3) or (iefrf == 7) or (iefrf == 10) or (iefrf == 11) or (iefrf == 13): ecrh = True axis.text(-0.05, 1, "Electron Cyclotron Current Drive:", ha="left", va="center") if iefrf == 12: @@ -2301,7 +2301,13 @@ def plot_current_drive_info(axis, mfile_data, scan): if (iefrffix == 5) or (iefrffix == 8): secondary_heating = "NBI" - if (iefrffix == 3) or (iefrffix == 7) or (iefrffix == 10) or (iefrffix == 11): + if ( + (iefrffix == 3) + or (iefrffix == 7) + or (iefrffix == 10) + or (iefrffix == 11) + or (iefrffix == 13) + ): secondary_heating = "ECH" if iefrffix == 12: secondary_heating = "EBW" diff --git a/source/fortran/current_drive_variables.f90 b/source/fortran/current_drive_variables.f90 index adbe82f375..4d30d4e6c6 100644 --- a/source/fortran/current_drive_variables.f90 +++ b/source/fortran/current_drive_variables.f90 @@ -68,8 +68,14 @@ module current_drive_variables !! current drive efficiency (A/W) real(dp) :: harnum - !! cyclotron harmonic frequency number, used in EBW cut-off + !! cyclotron harmonic frequency number, used in cut-off function + integer :: wave_mode + !! Switch for ECRH wave mode : + !! + !! - =0 O-mode + !! - =1 X-mode + real(dp) :: enbeam !! neutral beam energy (keV) (`iteration variable 19`) @@ -246,7 +252,8 @@ subroutine init_current_drive_variables echpwr = 0.0D0 echwpow = 0.0D0 effcd = 0.0D0 - harnum = 1.0D0 + harnum = 2.0 + wave_mode = 0 enbeam = 1.0D3 etacd = 0.0D0 etacdfix = 0.0D0 diff --git a/source/fortran/input.f90 b/source/fortran/input.f90 index 22b23caf11..78501cb1d8 100644 --- a/source/fortran/input.f90 +++ b/source/fortran/input.f90 @@ -260,7 +260,7 @@ subroutine parse_input_file(in_file,out_file,show_changes) use current_drive_variables, only: pinjfixmw, etaech, pinjalw, etanbi, & ftritbm, gamma_ecrh, pheat, rho_ecrh, beamwd, enbeam, pheatfix, bscfmax, & forbitloss, nbshield, tbeamin, feffcd, iefrf, iefrffix, irfcd, cboot, & - etalh, frbeam, harnum, xi_ebw + etalh, frbeam, harnum, xi_ebw, wave_mode use divertor_variables, only: fdfs, anginc, divdens, divclfr, c4div, & c5div, ksic, fififi, flux_exp, divplt, delld, c2div, beta_div, betao, divdum, tdiv, c6div, & omegan, prn1, fgamp, frrp, xpertin, c1div, betai, bpsout, xparain, fdiva, & @@ -1158,7 +1158,10 @@ subroutine parse_input_file(in_file,out_file,show_changes) 'User input ECRH gamma_CD') case ('harnum') call parse_real_variable('harnum', harnum, 1.0D0, 10.0D0, & - 'cyclotron harmonic frequency number') + 'Cyclotron harmonic frequency number') + case ('wave_mode') + call parse_int_variable('wave_mode', wave_mode, 0, 1, & + 'Cyclotron wave mode switch') case ('rho_ecrh') call parse_real_variable('rho_ecrh', rho_ecrh, 0.0D0, 1.0D0, & 'normalised minor radius at which electron cyclotron current drive is maximum') @@ -1166,10 +1169,10 @@ subroutine parse_input_file(in_file,out_file,show_changes) call parse_real_variable('xi_ebw', xi_ebw, 0.0D0, 1.0D0, & 'User input EBW scaling for Plasma Heating') case ('iefrf') - call parse_int_variable('iefrf', iefrf, 1, 12, & + call parse_int_variable('iefrf', iefrf, 1, 13, & 'Switch for curr drive efficiency model') case ('iefrffix') - call parse_int_variable('iefrffix', iefrffix, 0, 12, & + call parse_int_variable('iefrffix', iefrffix, 0, 13, & 'Switch for 2nd curr drive efficiency model') case ('irfcd') call parse_int_variable('irfcd', irfcd, 0, 1, &