Proceedings of ITC/ISHW2007

Neoclassical viscosities in NCSX and QPS
with few toroidal periods and low aspect ratios
S.Nishimura, D.R.Mikkelsena, D.A.Spongb, S.P.Hirshmanb,
L.P.Kua, H.E.Mynicka, M.C.Zarnstorff aaNational Institute for Fusion Science, 322-6 Oroshi-cho, Toki 509-5292, Japan
Princeton Plasma Physics Laboratory, Princeton, New Jersey 08543-0451, USAb
Oak Ridge National Laboratory, Oak Ridge, Tennessee 37831-6169, USA

Previously reported benchmarking examples for the analytical formulas of neoclassical viscosities were madeimplicitly assuming applications in a future integrated simulation system for the LHD (Large Helical Device).
Therefore the toroidal period numbers assumed there were mainly N=10. In this kind of calculation, however, an
implicit (or sometimes explicit) assumption of /N<<1 is sometimes included. This assumption is included not onlyin simplified bounce averaged drift kinetic equations for ripple diffusions, but also in the equation before theaveraging for non-bounce-averaged effects determining neoclassical parallel viscosity and the banana-plateau
diffusions. To clarify the applicability of the analytical methods even for configurations with extremely low
toroidal period numbers (required for low aspect ratios), we show here recent benchmarking examples in NCSX(National Compact Stellarator Experiment) with N=3 and QPS (Quasi-poloidal Stellarator) with N=2.
Keywords: neoclassical transport, neoclassical viscosity, moment equation approach, drift kinetic equation,
non-symmetric toroidal plasmas, low-aspect ratio advanced stellarators

1. Inroduction
The moment equation approach for neoclassicaltransport [1,2] in non-symmetric toroidal plasmas had
been developed mainly for neoclassical parallel flows andthe associated parallel viscosity [3-6]. Even though it was
shown in Ref.[7] that a consistent frame work including
not only the flows but also radial diffusions (in otherwords, not only the parallel viscosity, but also poloidal
and toroidal viscosities) can be constructed in this line ofmoment approach, methods to calculate all of required
viscosity coefficients in general collisionality regimes ingeneral toroidal configurations had not been shown.
Motivated by design activities of advanced stellarators, amethod to obtain the full neoclassical viscosity
coefficients was developed in Ref.[8]. It was shown therethat three mono-energetic viscosity coefficients
M*(parallel viscosity against flows), N*(driving force forbootstrap currents), and L*(radial diffusion) are required
to describe the full neoclassical characteristics of general
non-symmetric toroidal configurations. Since existing
numerical methods such as variational methods andMonte Carlo methods for the drift kinetic equation
described in the 3-D phase space (poloidal angle ,
toroidal angle , pitch angle ) could be applicable [8],applications of the new theory were done for various
types of advanced helical/stellarator configurations [9,10].
However, this step of the development of the moment
approach was still in the “basic frame work”. Evenauthor’s e-mail: nishimura.shin@lhd.nifs.ac.jp

though there were many alternative methods for the L* inthe collisionless limit (aforementioned ripple diffusions),
the other viscosity coefficients could be obtained only by
the Drift Kinetic Equation Solver (DKES) code
[8,9,10,11]. There are many demands for faster and easierestimation methods for the neoclassical quantities, in
integrated simulation systems using iterative calculationsof the equilibrium and the transport [12], in theconfiguration optimizations [9], and in experimental
studies investigating dependences on configurations.
Since the viscosities (or resulting neoclassical transportcoefficients) are direct consequences of guiding center
drift motions, evaluating them is an important part ofunderstanding the characteristics of the designed
magnetic configurations [9,10,13,14]. Even in tokamak
experiments, neoclassical toroidal viscosity effects due to
breaking the axisymmetry have recently been studied [15],and thus the framework of the moment approach for
non-symmetric configurations and the methods tocalculate the viscosity coefficients are now required forall studies of toroidal plasma confinements. In theories ofaxisymmetric tokamaks, simple analytical methods basedon asymptotic expansions of the drift kinetic equations
and connections of results of them are commonly used
[1,16]. Since this approach using analytical methods will
be required also in a future integrated simulation system
for the LHD(Large Helical Device)[12], we had
previously carried out derivations and benchmark tests ofProceedings of ITC/ISHW2007
the analytical formulas for the three mono-energetic
coefficients [17]. The previous benchmarking examples
were made implicitly assuming applications in the LHD[a helical heliotron with major and minor radii of
R0=3.9m, a=0.6m, and magnetic field strength of B03T].
Therefore the toroidal period numbers assumed there(Boozer)
were mainly N=10, and the assumed Bmn[B=Bmn cos(m nN )] spectra didn’t include n0,1. One
open
issue
remaining
there
is
the
ripple-trapped/untrapped boundary layer. Even though
there are many alternative calculating methods for thebounce-averaged motion of the ripple-trapped particles,
this boundary layer causes coupling effects between the
bounce-averaged motion of ripple-trapped particles andthe non-bounce-averaged motion of untrapped particles(collisional detrapping/entrapping). In Ref.[17], aprevious boundary layer theory by Shaing and Callen for
rippled tokamaks [18] was applied with an extension tomulti-helicity stellarators. One effect of this coupling,which was investigated also in Ref.[17], is a difference of2(BS)
N* (or G  B N*/M*) in the 1/ regime (Es/v 0)
regime (Es/v0,
from the collisionless detrapping/v0) value given by Refs.[3,4]. However, we had not
1/2
shown any benchmarking examples for the 1/
diffusion, which is another important effect discussed in
Ref.[18]. Even for the N* in the 1/ regime, thenumerical examples for more general cases had not beenshown. To investigate these effects in configurations
including Bmn of n0,1 and with extremely low toroidal
period numbers (required for low aspect ratios) giving
larger /N, we show in this paper recent calculation
examples in NCSX (National Compact Stellarator
Experiment)[10,13] and QPS (Quasi-poloidal Stellarator)[9,10,14].

2. Calculations in NCSXThe NCSX is a quasi-axisymmetric (QA) toroidal
system with N=3, R0=1.4m, a=0.32m, and B02T. Herewe consider the calculation on a flux surface withnormalized
toroidalfluxof( / edge)1/2=0.51(corresponding to r0.165m) in a standard configuration(NCSX-m50) with a finite beta of =4% and a finite
toroidal current of Ip=178 kA as an example. The
notations for the flux surface coordinates (mainly the
Boozer coordinates) in Refs.[8,17] are followed and thus
the radial derivatives of the poloidal and toroidal
magnetic fluxes are  ' =0.1178T•m, ' =0.2513T•m, andcovariant poloidal and toroidal components of themagnetic field are B =0.0036T•m, B =2.3210T•m,(Boozer)
in arespectively, on this flux surface. The Bmn
range of 0m16 and n 11 are used. As described in
Ref.[18], the boundary layer structure is determined by aPAS
drift kinetic equation (V//Ca )fa1=0, where V// PASisb• (μ=const) is the linearized orbit propagator and Ca
the pitch-angle-scattering operator with the collision
frequency of  Da [8]. Since V// fa10 for the non-trivial
solution of this equation, the existing bounce- or
ripple-averaging methods assuming V//fa1=0 is not
appropriate for this analysis. Therefore we should use thebounce- or ripple-averaging methods to obtain fa1/μ in
the ripple-trapped pitch-angle range, which gives theboundary condition for the boundary layer analysis [18],
together with the analytical solution for the boundary
layer as complimentary methods. For this analyticalsolution, a model expression of the magnetic fieldisstrength
B/B0=1+T( )+H( )cos{L N ( )}required.
We
use
here
 H ( )={Bmax ( )  Bmin ( )} / (2B0 ) for each poloidal
angle to define H( ). This is a truncation of the Fourier
series used by Todoroki [20] who expanded not theamplitude but the phase of B/B01T( ).
is
given
by
Analogously,
1+T( )
1 +  T ( )={Bmax ( ) + Bmin ( )} / (2B0 ) . The residual
ripple-well structure is distorted, or sometimes eliminated± /2, by finite rotational transform per toroidal
atperiod (’/ ’)/N in cases with small H. For this kind of
situations, the effective ripple well depth eff and effective
ripple well length correction * were introduced in thetheory for rippled tokamaks [18]. We use also thistechnique with an extension to helical/stellarator
configurations. The error of a well-known Shaing-Hokin
formula [21] for the 1/ ripple diffusions in H0 limits(for e.g., H0.01), which was pointed out in Ref.[8], isstrongly reduced by introducing this method. By using
1/2
diffusionthese notations, an expression for the 1/
coefficient in Ref.[18], which is a contribution offorripple-trapped
pitch-angle
range
0 212
 {wμB0(1+Teff)}/(2μB0 eff), can be extended to aform
includingmoregeneralnon-symmetricconfigurations as
L*(1/2) = 2.92

2 v 
 2   Da 

1/2

2 3/4 V ' ( ')2  4 2 

1/2

  2sin 1 *
d
3/4
 eff
B0

  2
 N '  L  '  '  /









1/2

2
2
2
 
    T  1   * 2 T H + (1   * 2 )  H  





 
9



(1)
Here, for the aforementioned distribution function fa1/μ2
in 0 1 as the boundary condition, an analyticalsolution by Shaing and Hokin [21] is used to consider
analytically the dependence on the magnetic
3/4 1/2
N,
configurations. From this form of L*(1/2)   eff
we can understand that this component of the diffusions
3/2 0
N
can dominate over the 1/ diffusion of L*(1)   eff

Proceedings of ITC/ISHW2007[21] only in configurations with small ripple amplitude
eff and with small toroidal period numbers N, andtherefore it will appear in QA configurations rather than
the rippled tokamaks considered in Ref.[18]. In fact,
previous numerical results in CHS-qa [23] with N=2
1/2
showed a clear 1/ dependence of L* in a wide rangeof collisionality (/v).
NCSX-M50 s=0.26
20

15

(BS)

G

Analytical w/oboundary layer correction

105

G(BS)(analytical)
G(BS)(DKES)(Es/v=0)
G(BS)(DKES)(Es/v=1.0e-4T)
G(BS)(DKES)(Es/v=3.0e-4T)
G(BS)(DKES)(Es/v=1.0e-3T)
G(BS)(DKES)(Es/v=3.0e-3T)

0 -7
-6
-5
-4
-3
-2
-1
10
10
10
10
10
10
10
1
/v [m ]101

NCSX-M50
s=0.26

100

0

10

1

10

L* (DKES)(Es/v=0)
1/nu (NEO)
NEO+analitical 1/sqrt(nu)L* (analytical)(banana-plateau)
L*(DKES)(Es/v=1.0e-4T)
L*(DKES)(Es/v=3.0e-4T)
L*(DKES)(Es/v=1.0e-3T)
L*(DKES)(Es/v=3.0e-3T)

3. Calculations in QPS
The QPS adopts a quasi-poloidal configuration
reducing the radial drift of the trapped particles [19], with
N=2, R0=1m, a=0.3m, and B01T. The parameters of theflux surface, where the calculation examples are made,

10-1
L*
1 -2

-2

[m T ] 10

NEO+
analytical
1/sqrt(nu)

10-3
10-4

analytical(ban-plat)

dependences of the DKES results on the EB drift
parameter Es/v. The N* in NCSX is insensitive to Es/v
3even in the range of Es/v310 T since the 1/ diffusionof the ripple-trapped particles accompanying theboundary layer correction N*(boundary) in Eq.(14) in
Ref.[17] is strongly reduced in this configuration. In spiteof this reduction of L*(1) and accompanied N*(boundary),
we can see another boundary layer effect in Fig.1(a). The
N* given by the DKES transiently becomes larger at
/v~103m1 compared with the analytical formula. Thistransient increase is peculiar to QA configurations where
the 1/1/2 component becomes comparable or dominates
over the 1/ component in the radial diffusion, and thus
was found also in CHS-qa [23]. Although this effect in
1/2
the 1/ regime cannot be calculated by a method in
Ref.[17] assuming a collisionless limit of the 1/ regime(This previous formula gives too small values for the QA
configurations and thus is not included in Fig.1(a)), thetransient increase, which is about 30% at most, will not
be so important in the energy integrated coefficients.

NEO

10-5
10-7 10-6 10-5 10-4 10-3 10-2 10-1
1
/v [m ]

12
10

0

1

QPS s=0.238
G(BS) E00000
G(BS) E00010
G(BS) E00030
G(BS) E00100
G(BS) E00300
G(BS)(analytical)
G(BS) (1/nu)

8
6

G

(BS)

Fig.1
Mono-energetic viscosity coefficients in the
NCSX given by the analytical methods (solid curves) andby the numerical method in the 3-D phase space (DKES)(open
symbols).
(a)thegeometricalfactor2(BS)
G B N*/M*, (b) components of the diagonal
diffusion L*.

42

 regime

0
(BS)

-2
10-6

G
10

-5

-4

10

2

=<B >N*/M*
10

-3

10-2
1

/v [m ]

102

10-1

1 2

100

101

QPS s=0.238
L* E00000
L* E00010
L* E00030
L* E00100
L* E00300
L*
1/sqrt(nu)1/nu (S-H)

101

L* [m T ]Figure 1(a) shows the mono-energetic viscosity
coefficients N* in the NCSX obtained by the analyticalformulas [17] and those by the DKES [11] with theconversion formulas in Ref.[8]. The mono-energetic
coefficient L* (Es/v0) is analytically given by sum of
three components: (1) L*(1) given by appropriate
bounce-averaging codes with field line integral methods,(2) L*(1/2) given by Eq.(1), and (3) contributions of
non-bounce-averaged drifts given by Eq.(16) in Ref.[17](L*(banana-plateau)). We used here the NEO code [22] for
the L*(1) in the NCSX, and Figure 1(b) shows thesecomponents L*(1), L*(1)+L*(1/2), L*(banana-plateau), andthe DKES results. The sum L*(1)+L*(1/2) approximately
predicts a deviation of the DKES from a pure 1/
scaling given by the bounce-averaging codes at
3 1
/v<10 m . In these figures, we show also the1/ regime

10

10

100
10-1
10-2
10-3 -6
10analytical cal. for ban-plat component
without the boundary layer correction
10

-5

10-4

10-3

10-2
1

10-1

100

101

/v [m ]Fig.2 Mono-energetic viscosity coefficients in the QPS
given by the analytical methods (solid curves) and by the
DKES (lines with symbols). (a) the geometrical factor
2(BS)G B N*/M*, (b) components of the diagonal
diffusion L*.

Proceedings of ITC/ISHW2007
1/2

are ( /edge) =0.49 (corresponding to r0.14m),  '=0.0275T•m,  ' =0.1423T•m, B =0, and B =1.1403T•m.
(Boozer)ranges are 0m20 and n 20. Results in
The Bmnthe QPS are shown in Fig.2. In Fig.2(b), the 1/ diffusion
coefficient L*(1/) given by the Shaing-Hokin formula[21] including the minor modifications of B expression inSec.2 is shown to confirm a validity of following
discussions on the boundary layer correction. Even for
3/2
L*(1)   eff, the Shaing-Hokin theory still retains an
accuracy of factor 2. Therefore we can investigate
characteristics of the boundary layer correction on theparallel viscosity N*(boundary) with a weaker dependenceon eff and H by the analytical method. As confirmed in
Ref.[17], we have to interpret a previous ”1/ regime”
formula for the parallel viscosity derived by Shaing,etal.[3,4] and the N* connection formulae including it (redsolid curve in Fig.2(a)) as expressions for strong Es/v
1/2
limit (i.e., the  regime or the  regime discussed later).
The correct 1/ regime (Es/v0) value is given by adding
a boundary layer correction term N*(boundary) which wasneglected in Refs.[3,4]. Although the calculation ofN*(boundary) in the QPS (H/ 1.5T/ andM=±2.4rad) requires some minor modifications forEq.(14) in Ref.[17], they will be reported in a separated
article. In Fig.2(a), we showed the 1/ regime asymptotic(sym)(asym)+N*+N*(boundary) [17].
value of N* given by N*
It approximately predicts the numerical result for a weak4radial electric ranges of Es/v<110 T by the DKES.

4. Concluding remarks
The mono-energetic neoclassical viscosity coefficientsare investigated in two low aspect stellarator
configurations with contrasting design concepts. For M*,
N*(asym), (Xa, G Xa(asym)) defined in Ref.[17] due to purenon-bounce-averaged motions [3,4], the validity of theanalytically approximated formulas [8,17] has been
confirmed even in the NCSX and in the QPS. In the two
configurations, there are contrasting effects of theripple-trapped/untrapped boundary layer at  21 causing
coupling effects between the bounce-averaged motion ofripple-trapped particles and the non-bounce-averaged
motion of untrapped particles. The 1/1/2 ripple diffusion
L*(1/2) in the QA configurations is peculiar to theconfigurations with small ripples. In the ripple-trapped
pitch angle range 0 21 in these configurations, the1/1/2 component of the perturbed distribution function is
not negligible compared with the small 1/ component.
However, their effects on the ripple-untrapped
pitch-angle range  2>1 is not important. In contrast tothis, the boundary layer affects on this range  2>1 in
configurations without quasi-axisymmetry and makeother boundary layer corrections on the viscosity
coefficients; N*(boundary) appearing in the 1/ regime(Es/v0) and also L*(boundary) near the collisionality regime
boundary between 1/ and plateau regimes. Even though
the integration constant in 0 21 is negligible comparedwith a large 1/ component, boundary layer effects as adriving force of (eff)1/2 for the flows in  2>1 is not
negligible for the 0 component of the distributionfunctionin
these
configurations
without
quasi-axisymmetry.

Acknowledgements
A main part of this work was done during one author(S.N.)'s visit to the Princeton Plasma Physics Laboratory(PPPL) and the Oak Ridge National Laboratory (ORNL)
sponsored by the Graduate University for Advanced
Studies (GUAS). A NIFS coordinated research program
NIFS07KNXN103 also supports this work.

References[1] S.P.Hirshman and D.J.Sigmar, Nucl.Fusion 21, 1079 (1981)[2] K.C.Shaing and J.D.Callen, Phys.Fluids 26, 3315 (1983)[3] K.C.Shaing, E.C.Crume,Jr.,et al., Phys.Fluids B 1,148(1989)[4] K.C.Shaing, B.A.Carreras et al., Phys.Fluids B1, 1663(1989)[5] N.Nakajima and M.Okamoto, J.Phys.Soc.Jpn.61, 833 (1992)[6] N.Nakajima, M.Okamoto, and M.Fujiwara.
J.Plasma Fusion Res.68, 46 (1992)[7] H.Sugama and W.Horton, Phys.Plasmas 3, 304 (1996)[8] H.Sugama and S.Nishimura, Phys.Plasmas 9, 4637 (2002)[9] D.J.Strickler, S.P.Hirshman, D.A.Spong, et al.,Fusion Sci. Technol.45,15(2004)[10] D.A.Spong, Phys.Plasmas 12, 056114 (2005)[11] W.I.van Rij and S.P.Hirshman, Phys.Fluids B 1,563 (1989)[12] Y.Nakamura, M.Yokoyama, N.Nakajima, et al.,Fusion Sci.Technol. 50, 457 (2006)[13] D.R.Mikkelsen, H.Maassberg, M.C.Zarnstorff, et al.,Fusion Sci.Technol. 51,166 (2007)[14] D.A.Spong, S.P.Hirshman, J.F.Lyon, et al.,Nucl.Fusion 45,918(2005)[15] W.Zhu, S.A.Sabbagh, R.E.Bell, et al.,Phys.Rev.Lett.96, 22002 (2006)[16] W.A.Houlberg, K.C.Shaing, S.P.Hirshman, and
M.C.Zarnstorff, Phys.Plasmas 4, 3230 (1997)[17] S.Nishimura, H.Sugama, and Y.Nakamura,
Fusion Sci.Technol.51,61(2007)[18] K.C.Shaing and J.D.Callen, Phys.Fluids 25, 1012 (1982)[19] H.E.Mynick, Phys.Fluids 26, 1008 (1983)[20] J.Todoroki, J.Phys.Soc.Jpn 8, 2758 (1990)[21] K.C.Shaing and S.A.Hokin, Phys.Fluids 26, 2136 (1983)[22] V.V.Nemov, S.V.Kasilov, W.Kernbichler, and M.F.Heyn,
Phys.Plasmas 6, 4622 (1999)[23] S.Okamura, K.Matsuoka, S.Nishimura, et al.,Nucl.Fusion 44, 575 (2004)