Proceedings of ITC/ISHW2007

Recent progress in NEO-2 - A code for neoclassical transportcomputations based on field line tracing∗
Winfried KERNBICHLER1) , Sergei V. KASILOV1,2) , Georg O. LEITOLD1) , Viktor V. NEMOV1,2) andKlaus ALLMAIER1)
1)

Association EURATOM-ÖAW, Institut für Theoretische Physik - Computational Physics, TU Graz, Petersgasse 16,A–8010 Graz, Austria
2)Institute of Plasma Physics, National Science Center “Kharkov Institute of Physics and Technology”, AkademicheskayaStr. 1, 61108 Kharkov, Ukraine
NEO-2 is a code for the computation of neoclassical transport coefficients and current drive efficiency
in toroidal devices which is based on the field line integration technique. The possibility to use the completelinearized collision integral is realized in this code. In this report the results of comparison of the full matrixof transport coefficients in a tokamak with analytical models are presented. Effects of simplifications of thelinearized collision model (e.g., reduction to a Lorentz model) are studied in order to provide a comparison withvarious momentum correction techniques used for the computation of transport coefficients in stellarators.
Keywords: neoclassical transport, linearized collison operator
Accurate computations of transport coefficients, bootstrap current and the generalized Spitzer function in tokamaks and stellarators is an important problem for stellarator optimization, generation of neoclassical data bases, andmodeling of current drive. Based on the field line integration technique [1], the code NEO-2 has been developed forthis purpose [2, 3]. This code solves the linearized drift kinetic equation in regimes where the effect of electric fieldon transport and bootstrap coefficients is negligible. Recently this code has been upgraded for computations of
the full transport matrix and the possibility of treatmentof magnetic fields in Boozer coordinates has been realized
in addition to magnetic fields in real space coordinates. In
the following, a comparison of NEO-2 results with resultsof analytical theory for tokamak are presented. Transportcoefficients are computed here only for the electron component assuming the ions to be immobile.
Approximating the energy dependence of the perturbation of the distribution function by the expansion over afinite number of test functions,MX
fm (ψ, s, λ)ϕm (v/vT ) , (1)δ f (ψ, s, v, λ) ≈ f0 (ψ, v)
m=0ϕm (x)= π3/4s2Γ(m + 1) (3/2) 2
L(x ) ,Γ(m + 5/2) m

where f0 (ψ, v) is the Maxwellian and L(3/2)
are associatedmLaguerre polynomials (Sonine polynomials) of the order∗ This work, supported by the European Communities under the contract of Association between EURATOM and the Austrian Academy ofSciences, was carried out within the framework of the European Fusion
Development Agreement. The views and opinions expressed herein do
not necessarily reflect those of the European Commission. Additional
funding is provided by the Austrian Science Foundation, FWF, under contract number P16797-N08.
author’s e-mail: winfried.kernbichler@TUGraz.at

3/2, the linearized drift kinetic equation (DKE) is transformed to a set of coupled two dimensional differentialequations,σ)M (X∂ fmσ1
νmm′ L fmσ′ + Kmm′ fmσ′ + Dmm′ fmσ′−κ∂s|λ|m′ =0ρ (3) σσ(2)(2)
A1 a(1)=
m + A2 am qG + A3 am qEh|∇ψ|iwith the pitch-angle scattering operator
L fmσ′==∂
1 ∂(1 − λ2 ) fmσ′ (ψ, s, λ)2|λ| ∂λ∂λ!
∂ η|λ| ∂ σf ′ (ψ, s, η) ,2∂η B̂ ∂η m(3)and the integral part of the linearized collision operator
LKmm′fmσ′ =1 X ℓ
I ′ Pℓ (λ)|λ| ℓ=0 mm

Z 1−1′dλ′ Pℓ (λ′ ) fmσ′ (ψ, s, λ′ ) , (4)where Pℓ are Legendre polynomials. Here, ψ is a flux
surface label, s is the distance counted along the m.f.l.,
λ = vk /v is pitch, σ is the sign of vk , η = (1 − λ2 )/ B̂
is the normalized perpendicular invariant (magnetic moment), B̂ = B/B0 is the magnetic field module normalized
to a reference magnetic field B0 which is equal to (0, 0) harmonic of magnetic field expansion in Boozer coordinates,
B0 = hB3 i/hB2 i, κ = 2/lc with lc = vT τee being the mean
1/2
free path, ρ = vT /ωc , vT = (2T e /me ) , ωc = eB0 /me c and√τee = 3m2e v3T / 16 πne e4 Λ are electron Larmor radius,
thermal velocity, cyclotron frequency and collision time,respectively, and h. . .i means the average over the volume
between two neighboring magnetic surfaces (flux surface(i)average). The quantities νmm′ , Imm′ , Dmm′ and am
are matrix elements independent of plasma parameters, whereas

Proceedings of ITC/ISHW2007
the quantities Ai are the thermodynamic forces defined as,!
1 ∂n
3 ∂Te ∂ΦA1 =−+h|∇ψ|i,n ∂ψ 2T ∂ψ T ∂ψA2 =1 ∂Th|∇ψ|i,T ∂ψA3 =e hEk B̂i.
T h B̂2 i(5)The source terms qGσ and qσE with drives by gradients and byparallel electric field, respectively, are defined as follows,

qGσ=V̂G=!
∂ |λ|
V̂G ,
qσE = σ B̂,∂η B̂!
1 4− η |∇ψ|kG ,
3 B̂

(6)
HH
γ33
=

where kG is the geodesic curvature.
Equations (2) are discretized over the η variable on anadaptive non-equidistant grid and the resulting set of coupled ordinary differential equations is solved with help ofintegration along the field lines. In particular, all necessary
flux surface averages are computed during this integration
as follows,Z L
Z L
hαidsds= lim
α /β.
(7)
L→∞
hβi0 B0 B
As a result of computation, one obtains the matrix of transport coefficients Leij defined through the thermodynamicforces A j and thermodynamic fluxes Ii via
Ii = −

3X

Leij A j ,dreff = dV/S , V is a volume limited by a flux surface and S
is a flux surface area. In order to obtain these matrices fordifferent definitions of plasma radius, e.g. for the radiusdefined via the toroidal flux, rψ = (2ψ/B00 )1/2 , coefficients
α1 and α2 should be multiplied with drψ /dreff .
For the large aspect ratio tokamak with circular flux
surfaces coefficients γi j computed by NEO-2 are compared to the analytical results of Refs. [4, 5, 6]. In particular, dimensionless transport coefficients for the HinHH
ton and Hazeltine model [4], γiHH
=j , are given by γi j−1/2
2 −3/2HH
Ki j q ǫt Zeff for i, j = 1, 2, γi j = Ki j qǫt /2 for i = 1, 2and j = 3 or i = 3 and j = 1, 2,(8)j=1where I1 = Γe , I2 = Qe /T e and I3 = h jk B̂i/e. Here Γeand Qe are average particle and heat flux densities defined
as total particle and heat fluxes divided by the flux surface
area and jk is parallel electron current density.
Symmetric matrix Leij is conveniently expressed
through dimensionless coefficients γi j which depend only
on the device geometry, the mean free path lc and effective
charge number Z as follows, Li j = ne γi j βi β j /τee , whereβ1 = β2 = ρ and β3 = lc . For this purpose solution to (2)is formally presented as a superposition of thermodynamicforces,ρ  ˆσ,(1)fmσ =A1 fm + A2 fˆmσ,(2) + A3 fˆmσ,(3) .
(9)h|∇ψ|i
Then γi j are determined by normalized solutions for single
drive problems, fˆmσ,(i) , as follows+* Z 1/B̂
αi α j X X (i)σ,( j) −σˆ(10)dη fm qib B̂
γi j =lc m σ=±1 m0

where qσ1 = qσ2 = qGσ , qσ3 = qσE , α1 = α2 = lc h|∇ψ|i−1 ,(i)α3 = 1 and numerical coefficients bm
are, again, independent of problem parameters. With help of these coefficientsthe quantity λbb [3] is expressed as λbb = 2γ31 .
It should be noted that matrices Leij and γi j correspondto the effective radius reff used as a radial variable whereK33 ǫt1/2 − 1
2Zeff 0.29 + 0.46 (1.08 + Zeff )−1(11)and matrix Ki j is defined by Eqs. (6.125) and (6.126) ofRef. [4]. Here q is safety factor and ǫt = r/R is the inverse aspect ratio. The results of the comparison are presented in Figs. 1 through 6 for ι = 1/q = 0.362 andǫt = 0.075. The results of NEO-2 are computed with associated Laguerre polynomials up to forth order and Legendre polynomials up to third order. For all coefficients thedependence on collisionality as well as Zeff is well reproduced. The main differences come from the finite toroidicity. Whereas NEO-2 does not assume smallness of themagnetic field modulation, theoretical approximations are
based on the expansion over ǫt . It should be noted that inNEO-2 all nine transport coefficients are computed independently and Onsager symmetry of these coefficients isused for the control of the computation accuracy which
improves with both, grid resolution and (mainly) number of Laguerre polynomials in modeling energy dependence. For the present computation violation of symmetries γ13 = γ31 and γ23 = γ32 is around 1% and violation ofsymmetry γ12 = γ21 is around 10%.
Beside the full linearized collision operator, two different model operators are used in Figs. 7 through 9,
namely, the mono-energetic collision model and monoenergetic collision model with momentum recovery. Lastℓtwo models are obtained by putting in (2) Dmm′ = Imm′ = 0
or only Dmm′ = 0, respectively. In particular, the monoenergetic model here corresponds to the most common
mono-energetic approach where transport coefficients are
given by the convolution over energy of the results for
the Lorentz model. It can be seen that mono-energetic
model overestimates particle diffusion coefficient γ11 whilemono-energetic model with momentum recovery underestimates this coefficient as compared to the full linearized
collision model. The bootstrap coefficient γ31 , in turn,is underestimated by the mono-energetic model while themono-energetic model with momentum recovery overestimates it. Finally, conductivity coefficient is well reproduced by the mono-energetic model while the monoenergetic model with momentum recovery significantly
overestimates it. Differences between all three models areProceedings of ITC/ISHW2007
441010

F: Zeff=1

F: Zeff=1

HH
AS
F: Z =2

HH
AS
F: Z =2eff

eff

HH
AS
F: Zeff=4310HH
AS
F: Zeff=4310HHAS2γ12γ11HHAS101110100102100−810−610−4−2101001010210−810Lc/lc−610−4−21010021010Lc/lc

Fig. 1 Results of NEO-2 with full linearized collision operator(F) and analytical models of Ref. [4] (HH) and Refs. [5,6] (AS) for the dimensionless diffusion coefficient γ11 at
three values of the effective charge Zeff .

Fig. 2 The same as in Fig.1 for γ12 .

410F: Zeff=1
HHASF: Z =2eff

HHAS3

γ2210210110−810−610−4−21010021010Lc/lc

Fig. 3 The same as in Fig.1 for γ22 .

12F: Z =1eff

HHAS
F: Zeff=210

HHAS
F: Z =4eff

HHAS

831

[1] V. V. Nemov, S. V. Kasilov, W. Kernbichler, and M.F. Heyn,Phys. Plasmas 6, 4622 (1999).
[2] W. Kernbichler, S. V. Kasilov, G. O. Leitold, V. V. Nemov,and K. Allmaier, Computation of neoclassical transport instellarators with finite collisionality, Proc. of the 15th Int.
Stellarator Workshop, Madrid, 3-7 October 2005, P2-15(2005).
[3] W. Kernbichler, S. V. Kasilov, G. O. Leitold, V. V. Nemov,and K. Allmaier, 33rd EPS Conference on Plasma Physics,Rome, 19–23 June 2006, ECA Vol. 30I P-2.189 (2006).
[4] F. L. Hinton and R. D. Hazeltine, Rev. Mod. Phys 48, 239(1976)[5] O. Sauter, C. Angioni, and Y. R. Lin-Liu, Phys. Plasmas 6,2834 (1999).
[6] C. Angioni and O. Sauter, Phys. Plasmas 7, 1224 (2000).
[7] K .Allmaier et al., ICNTS - Benchmarking of Bootstrap Current Coefficients, this conference (2007).
[8] K. Allmaier et al., Computations of neoclassical transport instellarators using a δ f method with reduced variance, thisconference (2007).
[9] C. D. Beidler et al., ICNTS - Benchmarking of Momentum
Correction Techniques, this conference (2007).

HHAS
F: Zeff=4γ

naturally reduced with higher Zeff . Currently NEO-2 hasbeen and is being used for the benchmarking of various
methods for the computation of mono-energetic transportcoefficients and bootstrap coefficient [7,8] as well as momentum correction techniques [9].

6420−810−610−4−21010L /lc cFig. 4 The same as in Fig.1 for γ31 .

010210Proceedings of ITC/ISHW2007
2515

F: Z =1

F: Z =1HH
AS
F: Zeff=2

M: Z =1eff

effeff

L: Z =1eff

F: Z =2eff

HHAS
F: Z =420

M: Z =2eff

L: Z =2eff

eff

HHAS

F: Zeff=410

M: Z =4eff

15γγ3132

L: Zeff=410550−8−610−410−2101000−61021010−510−4−110Fig. 8 The same as in Fig. 7 for the bootstrap coefficient γ3110.7−210c cFig. 5 The same as in Fig.1 for γ32 .

0.810L /lc c0.9−310

L /l

3.5

F: Zeff=1

F: Z =1HHAS
F: Zeff=2

M: Zeff=1eff

3L: Z =1eff

F: Zeff=2HH
AS
F: Zeff=4

M: Z =2eff

2.5HH
AS

L: Z =2eff

F: Zeff=4
M: Z =4eff

2L: Zeff=4

γ33

γ33

0.60.51.50.410.30.5

0.2

0.1−8−610−410−21001021010c c410

F: Zeff=1
M: Zeff=1
L: Z =1eff

F: Z =2eff

M: Z =2eff

L: Z =2eff

F: Z =4310eff

M: Z =4effL: Z =4γ11eff

2101−510−4−310−510−4−31010−21010−210−110Lc/lc

Fig. 7 Dimensionless particle diffusion coefficient γ11 for thefull linearized collision operator (F), mono-energetic approach (L) and mono-energetic approach with momentum recovery (M).

−110Fig. 9 The same as in Fig. 7 for the conductivity coefficient −γ33Fig. 6 The same as in Fig.1 for −γ33 .

10−6100−610Lc/lc

L /l