Mixed Bohm/gyro-Bohm Transport Model


Glenn Bateman, Alexei Pankin, and Arnold Kritz
Lehigh University Physics Department
16 Memorial Drive East, Bethlehem PA 18015
bateman@fusion.physics.lehigh.edu
pankin@fusion.physics.lehigh.edu
kritz@fusion.physics.lehigh.edu

The mixed_Bohm_gyro_Bohm module is used to compute anomalous plasma transport coefficients using the Mixed-Bohm/gyro-Bohm transport model. The Mixed Bohm/gyro-Bohm anomalous transport model (also called the JET or JETTO transport model) contains Bohm and gyro-Bohm contributions. The Bohm contribution, which is linear in the gyro-radius, is a non-local transport model, in which the transport throughout the plasma depends on a finite difference approximation to the electron temperature gradient at the edge of the plasma. The gyro-Bohm contribution is a local transport model, which is proportional to the square of the gyro-radius. The version of the Mixed Bohm/gyro-Bohm model in this module is described in detail in Nuclear Fusion 38 (1998) 1013. The magnetic and flow shear stabilization are described in Plasma Physics and Controlled Fusion 44 (2002) A495.

There are two subroutines in this module:

The mixed_model subroutine computes arrays of effective diffusivities using arrays representing the profiles from the center to the edge of the plasma.

The mixed subroutine computes scalar effective diffusivities given the plasma properties at a given point, and also given the finite difference approximation to the electron temperature gradient at the edge of the plasma.

The version of the Mixed-Bohm/gyro-Bohm model that is used in this subroutine is described in Nuclear Fusion 38 (1998) 1013 by M. Erba, T. Aniel, V. Basiuk, A. Becoulet, and X. Litaudon. Both the electron and ion thermal diffusivities consist of two terms. One term has Bohm scaling

$$\chi^{\rm Bohm} \equiv \rho_s c_s q^2 \frac{ a ( d p_e / d r ) }{ p_e } \Delta_{T_e},$$ (1)

while the other term has gyro-Bohm scaling

$$\chi^{\rm gyro-Bohm} \equiv \frac{ \rho_s^2 c_s }{ a } \frac{ a ( d T_e / d r ) }{ T_e }.$$ (2)

The notation is described in Table . In the Bohm diffusivity expression, $\Delta_{T_e}$ is a finite difference approximation to the normalized temperature electron temperature difference at the plasma edge

$$\Delta_{T_e} \equiv \frac{T_e(r/a=0.8) - T_e(r/a=1)}{T_e(r/a=1)}.$$ (3)

The resulting anomalous ion and electron thermal diffusivities are constructed from the sum of these Bohm and gyro-Bohm terms, with empirically determined coefficients [3]

$$\chi_i^{\rm JET} = 1.6 \times 10^{-4} \chi^{\rm Bohm} + 1.75 \times 10^{-2} \chi^{\rm gyro-Bohm}$$ (4)

$$\chi_e^{\rm JET} = 8 \times 10^{-5} \chi^{\rm Bohm} + 3.5 \times 10^{-2} \chi^{\rm gyro-Bohm}$$ (5)

and the hydrogenic charged particle diffusivity is given by

$$D^{\rm JET} \propto \frac{ \chi_i \chi_e }{ \chi_i + \chi_e }$$ (6)

It should be noted that the Mixed-Bohm/gyro-Bohm model is not a local transport model -- it cannot be evaluated with just the local plasma parameters. The Bohm contribution to the core transport $\chi^{\rm Bohm}$ is proportional to a finite difference approximation to the edge electron temperature gradient through the term $\Delta_{T_e}$.

Table: Notation

Variable	Units	Meaning
$a$	m	minor radius (half-width) of plasma
$B_T$	Tesla	vacuum toroidal magnetic field at
		major radius $R$ along flux surface
$c_s$	m/s	$[k_b T_e / m_i]^{1/2}$ speed of sound
$D$	m$^2$/s	effective charged particle diffusivity
		charged particle flux divided by density gradient
$e$	C	electron charge
$I_p$	MA	toroidal plasma current
$k_b$		conversion from keV to Joules
$m_i$	kg	average ion mass
$n_e$	m$^{-3}$	electron density
$q$		magnetic q-value
$r$	m	minor radius (half-width) of each flux surface
$R$	m	major radius to geometric
		center of flux surface
$T_e$	keV	electron temperature
$T_i$	keV	ion temperature
$Z_{\rm eff}$		$\sum_s n_s Z_s^2 / n_e$ summed over each species
$\chi$	m$^2$/s	effective thermal diffusivity
		heat flux divided by density time temperature gradient
$\delta$		plasma triangularity
$\kappa$		plasma elongation
$\nu_{*}$		collision frequency divided by bounce frequency
$\rho_s$	m	gyroradius [ $c_s m_i / ( e B_T )$]
$\rho_*$		normalized gyroradius ($\rho_s / a$)

c@mixed_model.tex
c--------1---------2---------3---------4---------5---------6---------7-c
c
      module mixed_Bohm_gyro_Bohm
      private
      public mixed_model
c
c..physical constants
c
      real, parameter  :: zckb = 1.60210e-16 ! energy conversion factor [Joule/keV]
      real, parameter  :: zcme = 9.1091e-31  ! electron mass [kg]
      real, parameter  :: zcmp = 1.67252e-27 ! proton mass [kg]
      real, parameter  :: zce  = 1.60210e-19 ! electron charge [Coulomb]
c
      contains
c
      subroutine mixed_model (
     &   rminor,  rmajor,  tekev,   tikev,   q
     & , btor,    aimass,  charge,  wexbs
     & , grdte,   grdne,   shear
     & , t_e_kev_edge, rminor_edge, npoints
     & , chi_i_mix,  themix,   thdmix
     & , thigb,   thegb,    thibohm, thebohm
     & , ierr
     & , lflowshear)
c
c
c    All the following 1-D arrays are assumed to be defined on flux
c    surfaces called zone boundaries where the transport fluxes are
c    to be computed.  The number of flux surfaces is given by npoints
c    (see below).  For example, if you want to compute the transport
c    on only one flux surface, set npoints = 1.
c
c  Input arrays:
c  -------------
c
c  rminor(jz)   = minor radius (half-width) of zone boundary [m]
c  rmajor(jz)   = major radius to geometric center of zone bndry [m]
c
c  tekev(jz)    = T_e (electron temperature) [keV]
c  tikev(jz)    = T_i (temperature of thermal ions) [keV]
c  q(jz)        = magnetic q-value
c  btor(jz)     = ( R B_tor ) / rmajor(jz)  [tesla]
c
c  aimass(jz)   = mean atomic mass of main thermal ions [AMU]
c               = ( sum_i n_i M_i ) / ( sum_i n_i ) where
c                 sum_i = sum over all ions, each with mass M_i
c
c  charge(jz)   = charge number of the main thermal ions
c                 = 1.0 for hydrogenic ions
c
c  wexbs(jz)    = ExB shearing rate in [rad/s]
c
c    All of the following normalized gradients are at zone boundaries.
c    r = half-width, R = major radius to center of flux surface
c
c  grdte(jz) = -R ( d T_e / d r ) / T_e
c  grdne(jz) = -R ( d n_e / d r ) / n_e
c  shear(jz) =  r ( d q   / d r ) / q    magnetic shear
c
c    The following variables are scalars:
c
c  t_e_kev_edge = electron temperature at the edge of the plasma [keV]
c  rminor_edge  = minor radius at the edge of the plasma [m]
c
c  npoints = number of values of jz in all of the above arrays [integer]
c
c  Control variable (input):
c  -------------------------
c
c  lflowshear = 0 for no magnetic and flow shear stabilization
c             = 1 to use magnetic and flow shear stabilzation by
c                 [T.J. Tala et al Plasma Phys. Controlled Fusion 44
c                 (2002) A495]
c
c  Output:
c  -------
c
c    The following effective diffusivities are given in MKS units m^2/sec
c
c  chi_i_mix(jz) = total ion thermal diffusivity from the MIXED model
c  themix(jz) = total electron thermal diffusivity from the MIXED model
c  thdmix(jz) = total hydrogenic ion diffusivity from the MIXED model
c
c    The following contributions to the effective diffusivities are
c  for diagnostic purposes:
c
c  thigb(jz) = gyro-Bohm contribution to the ion thermal diffusivity
c  thegb(jz) = gyro-Bohm contribution to the electron thermal diffusivity
c
c  thibohm(jz) = Bohm contribution to the ion thermal diffusivity
c  thebohm(jz) = Bohm contribution to the electron thermal diffusivity
c
c  ierr    = returning with value .ne. 0 indicates error
c
c***********************************************************************
c
c-----------------------------------------------------------------------
c
c  Compile this routine and routines that it calls with a compiler
c  option, such as -r8, to convert real to double precision when used on
c  workstations.
c
c-----------------------------------------------------------------------
c
c  External dependencies:
c
c  Call tree: mixed_model calls the following routines
c
c  mixed             - Computes diffusivity from the MIXED model
c
c-----------------------------------------------------------------------

      implicit none
c
c-----------------------------------------------------------------------
c..input variables
c
      integer, intent(in)  :: npoints ! number of radial points
c
      real, intent(in) ::
     &   rminor(*),  rmajor(*)
     & , tekev(*),   tikev(*),    q(*),       btor(*)
     & , wexbs(*),   aimass(*),   charge(*)
c
      real, intent(in) ::                    ! gradients
     &   grdne(*)
     & , grdte(*), shear(*)
c
      real, intent(in) ::  t_e_kev_edge, rminor_edge
c
c  optional input switch for flow-shear stabilization option
c
      integer, intent(in), optional :: lflowshear
c
c-----------------------------------------------------------------------
c..output variables
c
      real, intent(out) ::
     &   thigb(*),   thegb(*),    thibohm(*), thebohm(*)
c
      real, intent(out) ::
     &   chi_i_mix(*),   themix(*),  thdmix(*)
c
      integer, intent(out) :: ierr
c
c-----------------------------------------------------------------------
c..local variables
c
      integer  :: jz, j1, j2, jm
c
c  npoints1 = value of jz corresponding to 80 0f the normalized radius
c
      integer  :: npoints1, llflow

      integer  :: lswitch5
c
      real     :: zte_p8, zte_edge, zi
c
c..local variables connected to the mixed module
c

      real :: zq,    zsound,     zgyrfi,    zrhos, zra
     & , zchii, zchie, zdhyd
     & , zrmaj, zaimass, zcharge,   zbtor, zrmin
     & , zgte,  zti,  zte, zgne
     & , zrlpe, zshear, zgradte
     & , zchbe, zchbi, zchgbe, zchgbi
c
c.. variables for exb model
c
      real  zwexb
c
c  zwexb    = local copy of ExB shearing rate
c
c..initialize arrays
c
c
      thigb(1:npoints)  = 0.
      thegb(1:npoints)  = 0.
      thibohm(1:npoints)= 0.
      thebohm(1:npoints)= 0.
      chi_i_mix(1:npoints) = 0.
      themix(1:npoints) = 0.
      thdmix(1:npoints) = 0.
c
c..initialize switches
c
      if (present(lflowshear)) then ! flow shear correction
        llflow = lflowshear
      else
        llflow = 0
      endif
c
c-----------------------------------------------------------------------
c
c..physical constants
c
c     define the jz value at 0.8 0f normalized radius
c
      npoints1 = int(npoints*0.8)
c
c.. start the main do-loop over the radial index "jz"..........
c
c
      do 300 jz = 1, npoints
c
c  compute scalar quantities necessary for mixed module
c
        zshear   = shear(jz)
        zte      = tekev(jz)
        zti      = tikev(jz)
        zaimass  = aimass(jz)
        zcharge  = charge(jz)
        zrmin    = rminor(jz)
        zrmaj    = rmajor(jz)
        zgte     = grdte(jz)
        zgne     = grdne(jz)
        zgradte   = abs( ( zgte / zrmaj ) * zte)
        zbtor    = btor(jz)
        zq       = q(jz)
        zra      = rminor(jz) / rminor(npoints)
        zrlpe    = (zgte + zgne)
        zte_p8   = tekev(npoints1)
        zte_edge = tekev(npoints)
c
        zwexb = wexbs(jz)
c
c
      call mixed(
     &  zbtor,          zgradte,        zq,
     &  zra,            zrlpe,
     &  zrmaj,          zshear,         zte,            zte_p8,
     &  zte_edge,       zti,            zwexb,          zaimass,
     &  zcharge,        llflow,
     &  zchie,          zchii,          zdhyd,
     &  zchbe,          zchbi,          zchgbe,         zchgbi,
     &  ierr)
c
c
c  If ierr not equal to 0 an error has occured
c
        if (ierr .ne. 0) return
c
c  compute effective diffusivites for diagnostic purposes only
c
         thdmix(jz)  = zdhyd
         themix(jz)  = zchie
         chi_i_mix(jz)  = zchii
c
c  put value of gbohm and bohm terms into kb and rb terms for
c  diagnostic purposes
c
         thebohm(jz) = zchbe
         thibohm(jz) = zchbi
c
         thegb(jz)   = zchgbe
         thigb(jz)   = zchgbi
c
c..end of mixed model
c
c
 300  continue
c
c
c   end of the main do-loop over the radial index, "jz"----------
c
      return
      end subroutine mixed_model
c--------1---------2---------3---------4---------5---------6---------7-c
c

Subroutine for Computing Particle and Energy Fluxes
Using the JET Mixed Bohm/gyro-Bohm
Transport Model Version 1.2: 22 August 2003
Implemented by M. Erba, G. Bateman, A. H. Kritz, T. Onjun, and A. Pankin
Lehigh University

For questions about this routine, please contact:
Arnold Kritz, Lehigh: kritz@plasma.physics.lehigh.edu
Glenn Bateman, Lehigh: glenn@plasma.physics.lehigh.edu
Thawatchai Onjun, Lehigh: onjun@fusion.physics.lehigh.edu
Alexei Pankin, Lehigh: pankin@fusion.physics.lehigh.edu

c@mixed.tex
c--------1---------2---------3---------4---------5---------6---------7-c
c
      subroutine mixed(
     &  btor,           gradte,         q_safety,
     &  ra,             rlpe,
     &  rmaj,           shear,          tekev,          te_p8,
     &  te_edge,        tikev,          wexb,           aimass,
     &  zi,             lflow,
     &  chi_e,          chi_i,          d_hyd,
     &  chi_e_bohm,     chi_i_bohm,
     &  chi_e_gyro_bohm, chi_i_gyro_bohm,
     &  ierr)
c
c Revision History
c ----------------
c      date          description
c
c   25-Aug-2003      New version of flow shear sabilization correction
c                    Converted to Fortran-90 module
c   27-Jul-2000      Major rewrite by Bateman
c   12-Jun-2000      gyro-Bohm term restored
c   06-May-1999      Revised as module by Thawatchai Onjun
c   24-Feb-1999      First Version by Matteo Erba
c
c Mixed Transport Model (by M.Erba, V.V.Parail, A.Taroni)
c
c Inputs:
c    btor:      Toroidal magnetic field strength in Tesla
c                 at geometric center along magnetic flux surface
c    gradte:    Local electron temperature gradient [keV/m]
c    q_safety:  Local value of q, the safety factor
c    ra:        Normalized minor radius r/a
c    rlpe:      a/L_pe, where a is the minor radius of the plasma,
c                    and L_pe is the local electron pressure scale length
c    rmaj:      Major radius [m]
c    shear:     Local magnetic Shear
c    tekev:     Local electron temperature [keV]
c    tikev:     Ion temperature [keV]
c    te_p8:     Electron temperature at r/a = 0.8
c    te_edge:   Electron temperature at the edge [r/a = 1]
c    wexb:      EXB Rotation
c                    "Effects of {ExB} velocity shear and magnetic shear
c                    on turbulence and transport in magnetic confinement
c                    devices", Phys. of Plasmas, 4, 1499 (1997).
c    zi:        Ion charge
c
c Outputs:
c    chi_e:      Total electron thermal diffusivity [m^2/sec]
c    chi_i:      Total ion thermal diffusivity [m^2/sec]
c    d_hyd:      Hydrogenic ion particle diffusivity [m^2/sec]
c
c    chi_e_bohm:      Bohm contribution to electron thermal diffusivity
c                         [m^2/sec]
c    chi_i_bohm:      Bohm contribution to ion thermal diffusivity [m^2/sec]
c    chi_e_gyro_bohm: gyro-Bohm contribution to electron thermal diffusivity
c                 [m^2/sec]
c    chi_i_gyro_bohm:    gyro-Bohm contribution to ion thermal diffusivity
c                 [m^2/sec]
c    ierr:      Status code returned; 0 = OK, .ne.0 indicates error
c
c Coefficients set internally:
c
c    alfa_be:   Bohm contribution to electron thermal diffusivity
c    alfa_bi:   Bohm contribution to ion thermal diffusivity
c    alfa_gbe:  gyro-Bohm contribution to electron thermal diffusivity
c    alfa_gbi:  gyro-Bohm contribution to ion thermal diffusivity
c
c    coef1:     First coefficient for empirical hydrogen diffusivity
c    coef2:     Second coefficient for empirical hydrogen diffusivity
c
c Other internal variables:
c
c    gamma:     The characteristic growthrate for the ITG type of
c                    electrostatic turbulence
c    func:      Function for EXB and magnetic shear stabilization
c
      IMPLICIT NONE
c
c Declare variables
c
c..Input variables
c
      real, intent(in) ::
     &  btor,           ra,             rmaj,           rlpe,
     &  shear,          tekev,          te_p8,          te_edge,
     &  tikev,          gradte,         aimass,
     &  wexb,           zi,             q_safety
c
      integer, intent(in)  :: lflow
c
c..Output variables
c
      real, intent(out)    ::
     &  chi_e,               chi_i,               d_hyd,
     &  chi_e_gyro_bohm,     chi_i_gyro_bohm,
     &  chi_e_bohm,          chi_i_bohm
c
c
      integer, intent(out) :: ierr
c
c..Local variables
c
      REAL
     &  alfa_be,        alfa_bi,        alfa_gbe,       alfa_gbi,
     &  alfa_d,         chi0,           coef1,          coef2,
     &  delta_edge,
     &  em_i,           func,           gamma,
     &  omega_ce,       omega_ci,       rho,
     &  v_sound,        vte_sq,         vti,
     &  zepsilon
c
c
c..initialize diffusivities
c
      chi_e = 0.0
      chi_i = 0.0
      d_hyd = 0.0
c
      chi_e_bohm = 0.0
      chi_i_bohm = 0.0
      chi_e_gyro_bohm = 0.0
      chi_i_gyro_bohm = 0.0
c
c check input for validity
c
      zepsilon = 1.e-10
c
      ierr = 0
      if ( tekev .lt. zepsilon ) then
         ierr=1
         return
      elseif ( tikev .lt. zepsilon ) then
         ierr=2
         return
      elseif ( te_p8 .lt. zepsilon ) then
         ierr=3
         return
      elseif ( te_edge .lt. zepsilon ) then
         ierr=4
         return
      elseif ( rmaj  .lt. zepsilon ) then
         ierr=5
         return
      endif