Proceedings of ITC/ISHW2007

Modeling of Anomalous Transport in ECRH Plasmas at HSX
W. Guttenfelder1, D.T. Anderson1, J.M. Canik1, K.M. Likin1, J. Lore1, J.N. Talmadge1, W. Dorland2, M.
Barnes2
1

University of Wisconsin, Madison, USA2

University of Maryland, USA

The gyrokinetic code GS2 has been used to calculate the linear growth rates of the underlying micro-instabilities
presumably responsible for the anomalous transport in HSX. With ECRH heating (Te >> Ti), the dominant
long-wavelength instability is the trapped electron mode (TEM). To test whether available TEM transport models
can reproduce the transport in HSX, the Weiland ITG/TEM transport model is used, in addition to calculations ofneoclassical transport, to predict plasma profiles in HSX. To approximate the 3D geometry of HSX in the Weilandmodel, the necessary input geometry information is taken from the region in HSX where the fastest growing modes
are spatially localized in the 3D GS2 calculations. Specifically, the local curvature/∇B scale lengths (~R/3) andhelical ripple (εH) are used in place of the axisymmetric values (R and εT = r/R, respectively). With theseapproximations, the TEM linear growth rates predicted by the Weiland model agree quantitatively (within 30%)
with those predicted by the 3D GS2 calculations for HSX experimental parameters. Predicted density andtemperature profiles using the above transport estimates are in reasonable quantitative agreement with a number ofexperimental profiles in the QHS configuration. The predicted confinement times are within ~10% of theexperimental confinement times.
Keywords: turbulent transport model, quasi-linear transport, 3D gyrokinetic stability, trapped electron mode

1. Introduction
It has become routine to optimize stellarators using
neoclassical theory [1]. It is hoped that in the future,stellarators may also be optimized to reduce anomaloustransport, thought to be caused by plasma turbulence.
Indeed, there is evidence in LHD that anomaloustransport is reduced when neoclassical transport isreduced [2]. While significant advances have been madein the predictive capability for turbulent transport intokamaks [3], relatively little work has been performedfor 3D toroidal configurations such as stellarators.
Turbulent transport in tokamaks (caused by driftwaves) has been modeled through the use of quasi-linear
transport estimates that have been scaled to match
non-linear simulations [4-6]. Suppression of turbulent
transport via equilibrium E×B shear has also been
included to successfully model H-mode pedestals andinternal transport barriers [7]. These two features appear
to provide the dominant scaling in predicting turbulent
transport (via drift waves) in tokamak plasma.
Measured turbulence characteristics [8] and energy
confinement time scaling [9] are quite similar intokamaks and stellarators. Furthermore, recent non-linear
simulations [10,11] have demonstrated that predictedauthor’s e-mail: guttenfe@cae.wisc.edu

turbulence characteristics in stellarator geometries display
similarity to those predicted in tokamaks. However,
because of the 3D shaping, stellarator non-linear
simulations require increased resolution to treat thenon-symmetric geometry and various classes of trapped
particles, and are therefore more computationallyexpensive.
Given the similarity in turbulence and energy
confinement between stellarators and tokamaks, a logicalfirst step in performing predictive transport modeling forpresent stellarator experiments is to follow the same
framework as used for tokamak predictions. Following
this reasoning, this paper will present predictions ofdensity and temperature profiles in the HSX stellarator
using a tokamak transport model. To justify the use ofthis model for the 3D geometry in HSX, comparisons will
be made between linear growth rates predicted by thetokamak model, and those calculated using a 3D
gyrokinetic code that uses the HSX equilibrium.

2. 3D Gyrokinetic Microstability Calculations
The initial value gyrokinetic code GS2 [12] has beenpreviously used to calculate linear micro-stability in astellarator configuration [13]. To perform theseProceedings of ITC/ISHW2007
calculations, 3D equilibrium are calculated using VMEC
[14,15], which are then transformed into Boozer
coordinates using TERPSICHORE [16]. GS2 uses aballooning representation, and the necessary geometrycoefficients along a field line are calculated using VVBAL[17]. Since stellarators in general have multiple ripples in|B| along a field line, care must be taken in creating a gridthat captures these additional non-symmetric ripples.
Multiple grids have been created for the quasihelicallysymmetric (QHS) configuration of HSX with varying
number of grid points to determine how well these
additional ripples must be resolved to reach a converged
solution. For the calculations shown in this paper, 400 grid
points (over a range of θ = ±4π) were sufficient. Typically16 energy and ~30 pitch angle grid points were needed forthe calculations to converge. To verify accuracy, a few
cases were benchmarked against the independent linear
gyrokinetic FULL code for the HSX geometry [18]. The
GS2 calculations were performed on the NERSC IBM SP3
supercomputers using up to 256 processors, requiring afew minutes per run.
The linear stability calculations were performed on afield line where the magnetic axis undergoes the largest
excursion to the outboard side, and hence has the strongest
bad curvature drive. Figure 1 (top) shows eigenfunctions
calculated for the QHS configuration for multiple poloidal
wavenumbers (kθρs = 0.3-0.8), at a normalized minor

radius of ρ = 0.86, for experimentally relevant parameters(a/LTi = 0, Te/Ti = 2, a/Ln = a/LTe = 3, υ = 0). For theseconditions, the instability is the collisionless trappedelectron mode (CTEM), with a mode frequency
propagating in the electron diamagnetic drift direction.
Also shown in Fig. 1 are |B| along a field line (middle) andthe term proportional to the curvature drift (bottom). The
eigenfunctions for these unstable TEMs are strongly
localized in the low field, bad curvature region of thedevice, where the trapped particles exist. The structure ofthese instabilities is very similar to that observed intokamaks, due to the quasisymmetric field.

3. 1D Predictive Transport Modeling
To model the anomalous particle and electron heattransport in HSX [19], the Weiland model [4] is utilized.
The Weiland model is a toroidal fluid model for the ion
temperature gradient (ITG) and TEM instabilities. It is alinearized set of equations that provides both linear
stability predictions, and quasi-linear transport estimates ofparticle, and ion and electron heat transport. Thesetransport estimates have been checked against a limited
number of non-linear simulations [4,20] with good results.
The Weiland ITG/TEM model forms the core of theMulti-Mode Model [21] used for numerous tokamak
calculations.
As input, the Weiland model requires density andtemperature gradient scale lengths (Ln, LTe, LTi),temperature ratio (Te/Ti), wavenumber (kθρs), trapped
particle fraction (ft), and a ∇B/κ scale length (LB) toapproximate the toroidal drift terms. In a tokamak, thetrapped particle fraction is simply found using the toroidal
ripple, εT = r/R, f t = 2ε T (1 + ε T ) , and the ∇B/κ scale
length is the major radius, LB = R.
To use the axisymmetric Weiland model for HSX, twoapproximations must be made to account for the local 3D
geometry. These assumptions are based on the 3D
gyrokinetic calculations discussed in Sec. 2. Since theunstable eigenmodes are strongly localized in the low field,
bad curvature region (Fig. 1), the trapped particle fractionand ∇B/κ scale length are taken from this localized regionof HSX. Because of the quasihelical symmetry, the trapped
particle fraction is calculated using the dominant helicalripple (εH = 0.14⋅r/a = 1.4⋅r/R). In this location, where theaxis undergoes an excursion to the outboard side, the local
curvature is roughly three times larger than that of atokamak with the same major radius. Therefore, the ∇B/κscale length is reduced by a factor of three, LB = R/3.
The Weiland model has been used with the above
approximations to calculate the linear stability of HSX
plasmas. Growth rates have also been calculated for thesame input parameters using the 3D GS2 code for the QHSequilibrium. A comparison of growth rates is shown in Fig.
BB

BFig.1 Normalized eigenfunctions (top), |B| (center),and curvature drift (bottom) for GS2
calculations in HSX.

Proceedings of ITC/ISHW2007
2 for a scan over density and electron temperature gradient,with the other parameters the same as in Fig. 1. As seen in
the top two plots, the growth rates from the Weiland model
are close in magnitude to the 3D GS2 calculations. For therange of gradients indicative of experiment (highlightedred lines) the agreement is better than 30%. If the local
HSX geometry approximations are not used, the Weiland
model underpredicts the growth rates by about a factor oftwo (Fig. 2, bottom).
In order to predict density and temperature profiles,
1-D flux surface averaged transport equations for electron
density and temperature are solved numerically (Eqs. 1).
⎛⎞∂
1 ∂∂n2(1a)n+V' ⎜ − D∇ρ + V ( n ) n ∇ρ ⎟ = S(ρ)∂tV' ∂ρ⎜⎝⎟⎠∂ρ∑GS2 - HSX(1b)⎛⎞
1 ∂∂T13 ∂2n T+V ' ⎜ − nχ∇ρ + V ( nT ) nT ∇ρ ⎟⎟ = ∑ P(ρ )∂ρeV ' ∂ρ ⎜⎝2 ∂t⎠
A multi-mode model approach is used, summing transport
from the Weiland model (with geometry approximations),
neoclassical transport (calculated using DKES [22,23]),and a small transport contribution from a resistive
ballooning mode model, as used in [21]. For theneoclassical and Weiland model, the transport components
are represented by both diffusive (D,χ) and convective(V(n), V(nT)) components. The ECRH power deposition
profile is calculated using a ray-tracing code, and the total
absorbed power is determined from the time response of
the diamagnetic flux loop during ECRH turn off. Theparticle source rate profile used is based on 3D neutral gas
calculations [19]. The total magnitude is adjusted tominimize the difference in predicted and measured
densities, but is usually within a factor of two of the neutral
gas calculations that have been scaled to match absolutelycalibrated Hα measurements.
Using the above sources and transport models, Eqs.
1 are integrated to steady state. Figure 3 shows acomparison of predicted and measured electron density
and temperature profiles (using Thomson scattering) for
B=1.0T QHS plasmas for two different injected powers(Pinj = 44 & 100 kW, O1 ECRH). The density profiles
agree very well across the entire minor radius. The
temperature profiles agree outside a normalized minor
radius of ρ = 0.3, but the model underpredicts the coretemperature. For four different injected powers (includingtwo cases not shown here), the rms deviation of thesimulatedandexperimentaldensity
profileWeiland - HSX

experiment

Weiland - TOK

Fig.2 Growth rates (105 s-1) calculated in HSX using
3D GS2 (top), the Weiland model withgeometry approximations (center), and theWeilandmodel
without
geometry
approximations (bottom).
The highlighted
lines represent typical experimental gradients.

simulation

Fig.3 Predicted and measured electron density (top)and temperature (bottom) profiles for the QHS
configuration for Pinj = 44 kW (left) and 100
kW (right).

Proceedings of ITC/ISHW2007(∑ (nsim(ρ j ) − n expj)

2∑ (n ) , following ref. 3) is 9%.
exp 2jThe rms deviation for Te is considerably higher (40%).
The discrepancy of the core Te profiles is due to alarge χe from the Weiland TEM between ρ = 0.2-0.3, anda large neoclassical χe as the magnetic axis is approached.
It is of interest to note that in the region between ρ =0.2-0.3, the radial electric field predicted from theneoclassical ambipolarity constraint varies rapidly due toa change from the ion root to the electron root as the hotcore is approached. In this region, estimated E×B shear
rates [24] are larger than the linear TEM growth rates
predicted from the transport modeling. This may indicatethat E×B shear suppression could be important for
determining the appropriate anomalous transport
contribution.
Although the central Te is not predicted withaccuracy, the simulated energy confinement times agree
within 10% of the experimental confinement times.
Figure 4 shows the kinetic electron energy confinement
times measured experimentally (red) and predicted by thesimulations (blue) for four injected powers (26, 44, 70,100 kW). The predicted confinement times scale like
~P-0.57, similar to empirical scaling laws [9] and to that
expected from gyroBohm transport. Also shown are theenergy confinement times determined from thediamagnetic measurements for similar line-averaged
densities (black). These confinement times have a slightly
weaker scaling with absorbed power, ~P-0.42.

Fig.4 Predicted (blue) and experimental (red) kinetic
energy confinement times. Also shown are
confinement
times
fromdiamagnetic
measurements only (black).

4. Conclusions
The Weiland ITG/TEM anomalous transport modelhas been used to predict transport in the HSX stellarator.
By using approximations to represent the local geometry inHSX, the Weiland model predicts linear growth rates thatagree within 30% of those calculated using the 3D
gyrokinetic code GS2. Although non-linear effects or E×B
shear suppression have not been included, the linear
scaling should be captured reasonably well by thisapproximation. Predicted density profiles and energy
confinement times (using the Weiland model plus
appropriate neoclassical calculations) agree within 10% of
experimentally measured values. Root mean squaredeviations in the electron temperature profile are larger(40%) due to a systematic underprediction of the core Te.

References
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]

H.E. Mynick, Phys. Plasmas 13, 058102 (2006).
H. Yamada et al., Nucl. Fusion 41, 901 (2001).
ITER Physics Basis, Ch. 2 Nucl. Fusion 39, 2175 (1999).
H. Nordman et al., Nucl. Fusion 30, 983 (1990).
M. Kotschenreuther et al., Phys. Plasmas 2, 2381 (1995).
R.E. Waltz et al., Phys. Plasmas 4, 2482 (1997).
J.E. Kinsey et al., Phys. Plasmas 12, 052503 (2005).
A.J. Wooton et al., Plasma Phys. Contr. Fusion 34, 2023
(1992).
[9] U. Stroth et al., Nucl. Fusion 36, 1063 (1996).
[10] P. Xanthopoulos et al., Phys. Rev. Lett. 99, 035002 (2007).
[11] T.-H. Watanabe et al., Nucl. Fusion 47, 1383 (2007).
[12] M. Kotschenreuther et al,. Comp. Phys. Comm. 88, 128(1995).
[13] E.A. Belli et al., Bull. Am. Phys. Soc. 46, No. 8, 232(2001).
[14] S.P. Hirshman & D.K. Lee, Comp. Phys. Comm. 39, 161(1986).
[15] S.P. Hirshman et al., J. Comp. Phys. 87, 396 (1990).
[16] D.V. Anderson et al., Int. J. Supercomp. Appl. 4, 34(1990).
[17] W.A. Cooper, Plasma Phys. Contr. Fusion 34, 1011
(1992).
[18] G. Rewoldt et al., Phys. Plasmas 12, 102512 (2005).
[19] J.M. Canik et al., Phys. Rev. Lett. 98, 085002 (2007).
[20] A. Dimits et al., Phys. Plasmas 7, 969 (2000).
[21] G. Bateman et al., Phys. Plasmas 5, 1793 (1998).
[22] S.P. Hirshman et al., Phys. Fluids 29, 2951 (1986).
[23] W.I. van Rij & S.P. Hirshman, Phys. Fluids B 1, 563(1989).
[24] T.S. Hahm, Phys. Plasmas 4, 4074 (1997).