2994 new eccd current drive model to account for zeff and te with a cutoff function

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
+[^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. 

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(

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"

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

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,18 +1158,21 @@ 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')
        case ('xi_ebw')
 	  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, &