Continuum Edge Gyrokinetic Theory and Simulations
Transcript Of Continuum Edge Gyrokinetic Theory and Simulations
1
TH/P6-23
Continuum Edge Gyrokinetic Theory and Simulations1
X.Q. Xu 1), K. Bodi 2), J. Candy 3), B. I. Cohen 1), R. H. Cohen 1), P. Colella 4), M. R. Dorr 1), J. A. Hittinger 1), G. D. Kerbel 1), S. Krasheninnikov 2), W. M. Nevins 1), H. Qin 5), T. D. Rognlien 1), P. B. Snyder 3), M. V. Umansky 1), Z. Xiong 1)
1) Lawrence Livermore National Laboratory, Livermore, CA 94551 USA 2) University of California, San Diego, La Jolla, CA 92093 USA 3) General Atomics, San Diego, CA 92186 USA 4) Lawrence Berkeley National Laboratory, Berkeley, CA 94720 USA 5) Princeton Plasma Physics Laboratory, Princeton, NJ 08543 USA
e-mail contact of main author: [email protected]
Abstract. The following results are presented from the development and application of TEMPEST, a fully nonlinear (full-f) five dimensional (3d2v) gyrokinetic continuum edgeplasma code. (1) As a test of the interaction of collisions and parallel streaming, TEMPEST is compared with published analytic and numerical results for endloss of particles confined by combined electrostatic and magnetic wells. Good agreement is found over a wide range of collisionality, confining potential, and mirror ratio; and the required velocity space resolution is modest. (2) In a large-aspect-ratio circular geometry, excellent agreement is found for a neoclassical equilibrium with parallel ion flow in the banana regime with zero temperature gradient. (3) The four dimensional (2d2v) version of the code produces the first self-consistent simulation results of collisionless damping of geodesic acoustic modes and zonal flow (Ronsenbluth-Hinton residual) with Boltzmann electrons using a full-f code. In divertor geometry, it is found that the endloss of particles and energy induces pedestal-like density and temperature profiles inside the magnetic separatrix and parallel flow stronger than the core neoclassical predictions in the SOL. (4) Our 5D gyrokinetic formulation yields a set of nonlinear electrostatic gyrokinetic equations that are for both neoclassical and turbulence simulations.
1. Introduction
Understanding the structure of edge transport barrier in high-performance (H-mode) discharges requires a kinetic description of the plasmas because the radial width of the pedestal observed in experiments is comparable to the radial width of individual ion orbits (leading to large distortion of the local distribution function from a Maxwellian), while the ion and electron mean-free-paths are long compared to the connection length in the hot plasma at the top of the edge pedestal (violating the assumptions underlying a collisional fluid mode). Gyrokinetic formulation (2v) [1] is a reasonable approximation for edge plasmas because much pedestal physics is typically low frequency phenomenon compared to ion gyro-frequency. But existing gyrokinetic theories and codes do not apply to edge plasmas because they cannot treat fully nonlinear electromagnetic perturbations with multi-scale-length structures in space-time for full divertor geometry.
1This work was performed under the auspices of the U.S. Department of Energy by University of California Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48 at LLNL, Grant No. DE-FG0204ER54739 at UCSD, Grants DE-FG03-95ER54309 at general Atomics, and DE-AC02-76CHO3073 at PPPL.
2
TH/P6-23
We report on the development and application of TEMPEST, a fully nonlinear (full-f) gyrokinetic code, to simulate H-mode edge plasmas. This 5-dimensional (ψ, θ, ζ, E0, µ) continuum code represents velocity space via a grid in equilibrium energy (E0) and magnetic moment (µ) variables, and configuration space via a grid in poloidal magnetic flux (ψ), poloidal angle (θ) and toroidal angle (ζ). The geometry is that of a fully diverted tokamak and so includes boundary conditions for both closed magnetic flux surfaces and open field lines. A set of gyrokinetic equations [2, 3] are discretized on a computational grid for both circular and divertor geometry with a magnetic separatrix. The equations are solved via a Method-of-Lines approach and an implicit backward-differencing scheme using a NewtonKrylov iteration advances the system in time [4]. The spatial derivatives are discretized with finite differences while a high-order finite volume method is used in velocity space (E0, µ). A fourth-order upwinding algorithm is used for parallel streaming; and a fifth-order WENO scheme [5] is used for particle cross-field drifts. Boundary conditions at conducting material surfaces are implemented on the plasma side of the sheath. The code includes kinetic or Boltzmann electrons. A nonlinear Fokker-Planck collision operator (CQL) from STELLA [6] has been extracted and integrated into the gyrokinetic simulation code using the same implicit Newton-Krylov solver. A new Fokker-Planck collision operator in (E0, µ) space is under development for improved accuracy and conservation properties. The gyrokinetic Poisson equation is solved using GMRES with a multi-grid preconditioner.
2. Basic Gyrokinetic Equation
A set of generalized gyrokinetic Vlasov-Maxwell equations in the gyrocenter coordinate system valid for edge-plasma conditions has been derived by the Lie transform perturbation method to the Poincar´e-Cartan-Einstein 1-form and the pullback transformation for the distribution function [2]. This formalism allows inclusion of nonlinear large-amplitude, time-dependent background electromagnetic fields in addition to small-amplitude, shortwavelength electromagnetic perturbations. As an example, the pullback transformation in the gyrokinetic Poisson equation is explicitly expressed in terms of moments of the gyrocenter distribution function, thus describing the important gyro-orbit squeezing effect due to the large electric field shearing in the edge and the full finite Larmour radius effect for short wavelength fluctuations. The familiar polarization drift density in the gyrocenter Poisson equation is replaced by a more general expression.
2.1 Fully Nonlinear Ion Gyrokinetic Equations
The ion gyro-kinetic equations presently implemented in TEMPEST for the time-dependent five-dimensional (5D) distribution functions are simplified from our recent new formulation [2] and earlier Hahm’s work [3]. Due to the large potential and its multiple spatial-time scale nature in the edge pedestal, in order to accurately simulate particle parallel streaming, the field Φ is split into two parts Φ = Φ0 + δφ: Φ0 is the large amplitude, slowly varying component; δφ is the small amplitude, rapidly varying component. Here E0 is defined as the total energy including Φ0, but not δφ. Then E0 is a constant of motion if δφ ∼ 0 and a coordinate aligned with the direction of phase-space flow. The kinetic equation for the gyrocenter distribution function Fα(x¯, µ¯, E¯0, t) in gyrocenter coordinates (x¯ = x − ρα, ρα = b × v/Ωcα),
3
TH/P6-23
“total energy” E¯0, and magnetic moment µ¯, has the form:
∂Fα + v¯d · ∂Fα + (v¯ α + vBanos)b · ∂Fα (1)
∂t
∂x¯⊥
∂x¯
+ q ∂ Φ0 + µ¯ ∂B − B v¯ q ∂ δφ − v · (q∇¯ δφ ) ∂Fα = C(F , F ),
∂t
∂t B∗ ∂s
d0
∂E0
αα
v¯d = qcBb∗ × q∇¯ Φ + µ¯∇¯ B + v¯2 MqBα∗c (∇¯ ×b). (2)
v¯d0 = qcBb∗ × q∇¯ Φ0 + µ¯∇¯ B + v¯2 MqBα∗c (∇¯ ×b). (3)
2
v¯ = ±
(E − µ¯B − q Φ ), v
=
µc (b
·
∇¯ ×b),
(4)
Mα 0
0
Banos q
B∗ ≡ B 1 + b · (v ∇¯ ×b) , Ω = qB , µ = Mαv⊥2 ,
(5)
α
Ωcα
cα Mαc
2B
δφ = Φ − Φ0 .
(6)
Here Zαe, Mα are the electric charge and mass of electrons (α = e), ions (α = i). µ is the guiding center magnetic moment. The left-hand side of Eq. (1) describes the particle motion in the electric field and magnetic field. Cα is the Coulomb collision operator. The over-bar is used for the gyrocenter variables and denotes the gyroangle averaging. The additional E0 × B flow terms due to the large amplitude and the slow variation Φ0 from the complete formulation [2] will be added.
2.2 Fully Nonlinear Gyrokinetic Poisson equation
The complete gyrokinetic Poisson equation has been derived [2], including orbit squeezing by large Er shearing and full FLR effect. To make it numerically attractive, two additional approximations are made here: (1). the spatial variation of the transverse µ¯ moments Mn(x¯) calculated from Fα(x¯, µ¯, E¯0, t) is much slower than that of potential in evaluation of the full FLR effect; (2). the total transverse distribution function is Maxwellian.
2.2.1 Fully Nonlinear Gyro-kinetic Poisson equation in the arbitrary wavelength regime
In the arbitrary wavelength regime, the self-consistent electrostatic potential is computed from the gyro-kinetic Poisson equation:
1 0 = −4πe α ZαNα(x, t) − ne(x, t) − α λ2Dα [Γ0(b) − 1] Φ, (7)
where Γ0(b) = I0(b)e−b, b = ρ2α∇2⊥/2, I0(b) is the usual zeroth-order modified Bessel function.
The ion gyroradius is ρα = 2T⊥α/Mα/Ωα, the ion gyrofrequency is Ωα = ZαeB/Mαc, and the ion Debye length is λ2Dα = T⊥α/4πNαZα2e2. Although Eq. (7) is similar to the usual gyrokinetic Poisson equation [3], there is an important distinction. Our gyro-kinetic Poisson equation is fully nonlinear and the gyrocenter center density Nα and perpendicular ion pressure p⊥α are calculated from the gyrocenter distribution function Fα(x¯, µ¯, E¯0, t).
2π N (x, t) ≡
B∗dv¯ dµ¯F ,
2π n (x, t) ≡
B∗dv dµf ,
(8)
α
Mα
α
e
me
e
4
TH/P6-23
p = πB dv dµ¯(v2 F ), T = p⊥α
(9)
⊥α
⊥α
⊥α Nα(x, t)
Here the dot product between the density and potential as well as the Debye shielding have been dropped for simplicity.
The first-order Pad´e approximation to Γ0 with Γ0 − 1 = b/(1 + b) is an excellent fit for 0 ≤ b ≤ 9, and is therefore valid well into the typical ion gyrokinetic regime as used previously
in gyro-kinetic or gyro-fluid simulations [7, 8]. For single ion species, substituting a simple functional transformation Φ = φL + [T⊥α/(NαZα2e)] [ZαNα(x, t) − ne(x, t)] and the Pad´e approximation into Eq. (7) yields
ρ2α ∇2 φ = − Tα 1 + ρ2α ∇2 ln T⊥α [Z N (x, t) − n (x, t)] .
(10)
2 ⊥L
Nα Zα2 e
2 ⊥ Nα
αα
e
where φL is obtained by the gyrokinetic Poisson solver.
2.2.2 Fully Nonlinear Gyro-kinetic Poisson equation in the long wavelength regime
In the long wavelength limit k⊥ρα 1, the self-consistent electric field is typically computed from the gyrokinetic Poisson equation for multiple species
ρ2α
2
α 2λ2Dα ∇⊥Φ +
ρ2α ∇⊥ ln Nα α 2λ2Dα
· ∇⊥Φ + ∇2Φ
= −4πe
Z N (x, t) − n (x, t) − ρ2α 1 ∇2 p . (11)
αα
e
α
α 2λ2Dα NαZαe ⊥ ⊥α
The ion gyroradius is ρα = 2T⊥α/Mα/Ωα, the ion gyrofrequency is Ωα = ZαeB/Mαc, and the ion Debye length is λ2Dα = T⊥α/4πnαZα2e2.
There are two important distinctions between Eq. (11) and the usual gyro-kinetic Poisson equation [3]. Our gyro-kinetic Poisson equation is fully nonlinear and the gyrocenter center density Nα and perpendicular ion pressure p⊥α are calculated from the gyrocenter distribution function Fα(x¯, µ¯, E¯0, t), as defined in Eqs. (8)-(9). The last term of Eq. (11) is the diamagnetic density from the long wavelength expansion of the gyroaveraged gyrocenter density Nα(x, t). Although the diamagnetic density is small compared to the ion gyrocenter density, it is of the same order as both the polarization density in high-beta plasmas and as the difference between ion and electron gyrocenter density! This equation is an extension of typical neoclassical electric field calculation including poloidal variation.
2.3 Boundary conditions
2.3.1 Radial boundary conditions
The radial Robin boundary conditions are used for Fα and potential δφ at the inner core surface ψ = ψc and the outer wall surface ψ = ψw. This is a generalization of the Dirichlet and Neumann boundary conditions. Since the gyrokinetic equation has just a first-order radial advection term, there is only one physical boundary condition, depending on the direction of convection, No boundary condition should be imposed for particles convecting out of the simulation domain and therefore a extrapolation is used at that boundary.
5
TH/P6-23
2.3.2 Poloidal boundary conditions
The boundary conditions in the θ direction for Fα and for Φ are the sheath boundary conditions at the divertor plates and periodic in the “core” (closed field-lines). Our present implementation for sheath boundary conditions is for the case of normal intersection of the flux surface with the walls, where the ions are not confined for the perfectly absorbing wall and the current through the sheath is zero with no biasing; and there is an energetic group of impinging electrons that can escape the sheath potential and reach the wall with the energy E0 + eδφsh − µB > 0. Here Φsh = Φ0sh + δφsh is the sheath potential.
i. Sheath boundary conditions for potential
If the gyrokinetic ion and fluid electron model are used, the sheath potential is determined by the quasi-neutrality condition:
Φ = Te ln
4Γi,sh
,Γ
2πB ∞
=
dE
(E0−qΦ0sh)/B
dµ v
F σ.(12)
sh
e
i,sh
ne,shζ 8Te,sh/πme
Mα2 qΦ0sh
0 0
|v | i
The σ = ± represents the plus and minus sheet of parallel velocity with Fiσ = 0 for only incoming sheet. Here it is assumed that impinging electrons have a Maxwellian distribution.
The factor ζ ≡ 1/(1 + τp/τe) includes the correction of electron long mean-free paths physics. τp is long mean-free path confinement time and τc is the confinement time for the collisional sheath limited case. ζ ≡ 1 if the electrons are in the short mean-free paths regime.
If both electron and ion are kinetic, the sheath potential is determined:
Γ
2πB ∞
=
dE
(E0−qΦ0sh)/B
dµ v
F σ,
i,sh
Mα2 qΦ0sh
0 0
|v | i
Γ
2πB ∞
=
dE
(E0−eδφsh)/B
dµ v
F σ,
e,sh
m2e eδφsh
0 0
|v | e
Γi,sh = Γe,sh.
(13) (14)
ii. Sheath boundary conditions for distribution functions If the gyrokinetic ion and fluid electron model is used, the ion distribution function is:
Fα(ψ, θ, E0, µ) = Fα(ψ, θ, E0, µ), v ≥ 0
(15)
0,
v ≤0
A convention regarding the sign of the parallel velocity is that it is positive when there is a positive projection on the θ axis. Here positive θ axis is pointing to the plate/wall.
If both electron and ion are kinetic, the electron distribution function is:
fe(ψ, θ, E0, µ, σ = 1) = fe(ψ, θ, E0, µ, σ = 1),
(16)
fe(ψ, θ, E0, µ, σ = −1) = fe(ψ, θ, E0, µ, σ = 1), |v | ≤ vSH
(17)
0,
|v | ≥ vSH
Here vsh = 2eΦsh/me.
6
TH/P6-23
3. TEMPEST simulation Results
TEMPEST has been developed in the modern framework for either circular or divertor geometry. With importance of benchmark in mind, both full-f or δf options are available. TEMPEST is runnable as (1) 3D for the SOL endloss with F (θ, E0, µ), (2) 4D for transport F (ψ, θ, E0, µ), and (3) 5D for turbulence F (ψ, θ, ζ, E0, µ). The different aspects of the 3D, 4D and 5D TEMPEST have been verified on various known physics problems, such as (1) collisional scattering into a velocity-space loss cone, (2). neoclassical flow and transport, (3). electric field generation and geodesic acoustic mode damping, (4). drift waves and ion temperature gradient (ITG) modes.
3.1 3D Pastukhov collisional endloss
As a test of collisional velocity-space transport and parallel streaming, 3D TEMPEST (1d2v) simulation results are compared with published analytic and numerical results as shown in Fig. 1 for the endloss of particles confined by combined electrostatic and magnetic wells [9, 10, 11]. Here electrostatic and magnetic field are uniform with abrupt confining magnetic field and potential barriers at wall. Good agreement is found over a wide range of collisionality, confining potential, and mirror ratio; and the required velocity-space resolution is modest. In these simulations, the linearized CQL package is used.
a)
b)
c)
FIG. 1: Collisional endloss (“Pastukohov”) test cases: (a) confinement time versus density; (b) confinement time versus potential eφ/Te at low collisionality; (c) confinement versus mirror ratio at low collisionality.
3.2 4D Neoclassical tests
If using a shifted Maxwellian distribution that analytically satisfies Eq. (1) as an initial condition, then TEMPEST should preserve the solution without any significant change (within our finite difference truncation accuracy) after running some time steps. We tested such a case using the following simulation parameters: inverse aspect ratio = a/R0 = 0.03, the major radius R0 = 17.1 meter, toroidal magnetic field Bt = 1.5T , and poloidal
7
TH/P6-23
magnetic field Bp = 0.2T . The ion density and temperature profiles used are ni(ψ) = Nix exp(− ln(Nix/Nio)ψ/Lψ), Ti(ψ) = Tix exp(− ln(Tix/Tio)ψ/Lψ). with Tix = 3keV, Tio = 0.95Tix, Nix = 1 × 1020m−3, and Nio = 0.95Nix. The mesh resolutions are nψ = 30, nθ = 50, nE = 60, and nµ = 30. In the simulations Φ is set to be zero for simplicity. As shown in Fig. 2a), the simulation results remain in good agreement with theoretical prediction even after 10000 time steps ( 50 thermal ion transit time). The solid line in plot of q i comes from theoretical prediction for a shifted Maxwellian distribution, q α = 2.5NαU αTα, and U α = − ΩIα MTαα ∂ l∂nψNα .
In divertor geometry, for a given particle and heat sources on the inner core boundary surface, because the parallel transport is much larger than neoclassical radial transport, it is found that the endloss of particles and energy (discussed in section 2.3) induces (1) pedestal-like density and temperature profiles inside the magnetic separatrix; (2) parallel flow stronger than core neoclassical predictions in the SOL as indicated in Fig. 2 (b); (3) a symmetry point is developed for the parallel heat flux on the top of the machine as expected as shown in Fig. 2(c).
1.5×10 24
(eV/m2/s)
1.0×10 24
0.5×1024
10000 8000 6000 4000
(m/s)
theory
time 3 time 2 time 1 time 0
q (R,Z)(eV/m2/s) ||
2.5 2.0 Z(m) 1.5 1.0
separatrix
0 0.0
0.1 0.2 0.3 0.4 Radial Position (m)
a)
2000
0 0.0 0.02 0.04 0.06 0.08 Radial position (m)
b)
0.5
1.0 1.2 1.4 1.6 1.8 2.0 2.2 R(m)
c)
FIG. 2: (a) Comparison between simulation results with theory for a collisionless case with ∇Ti = 0 and zero finite banana orbit width. Flux surface averaged parallel heat flux q i . The solid lines comes theory and other lines from different time; (b) Comparison of simulation results with theory for flux surface averaged parallel flow velocity U i in the banana regime with ν∗i 0.02, ∇Ti = 0 and finite banana orbit width in X-point divertor geometry. The average is done by integration along the field-line from inner plate to outer plate in the SOL. (c) The contours of parallel heat flux q (R, Z) in the divertor geometry.
3.3 4D Geodesic-Acoustic Modes
The Geodesic-Acoustic Mode (GAM) is an asymmetric mode, which involves parallel ion dynamics, cross field drifts, and acceleration. Earlier GAM theory and simulations focused on large aspect ratio and small orbit [12, 13]. Recently Sugama and Watanabe find that the damping rate is sensitive to k⊥ρi at large q due to the effect of large banana orbit [14]. In our 4D GAM simulations, the charge is radially separated by an initial sinusoidal perturbation on
8
TH/P6-23
ion density. The electron model is Boltzmann ne = ni(ψ, θ, t = 0) exp(eφ/Te)/ exp(eφ/Te) , where represents the flux surface average. This choice of coefficient for Boltzmann electron model means that there is no cross field electron transport. Both radial and poloidal boundary conditions are periodic. The first full-f, self-consistent simulation results of collisionless damping of geodesic acoustic modes and zonal flow are plotted in Fig. 3. Good agreement is shown between theory [14] and simulations for the frequency of GAMs in Fig. 3a) and damping rate in Fig. 3b). The 30% difference between theory and simulation is attributed to the fact that the theory uses an asympototic 1/q2 expansion for large q, while in our simulation q=2.2 is moderate. The large effect of the orbit size on the GAM damping rate is illustrated in Fig. 3c). For the same parameters, the damping rate is almost zero if the finite banana orbit effect is ignored.
φ(t) / φ(t=0) -γ / (vt/qR)
1.0
Velocity mesh
(nE,nµ) = (25,50) (nE,nµ) = (45,50)
γTEMPEST = 0.072(vtd/R)
0.5
(nE,nµ) = (91,101)
Numerical residual
GAM damping rate 0.15
Sugama & Watanabe 06 0.10 kρp = 0.1
γSugama-Watanabe = 0.052(vtd/R) 0
Spatial mesh
(nψ,nθ) = (30,50)
-0.5
a/R0 = 0.02
q = 2.2
0
25
50
75
Time (R0/vtd)
a)
0.05
kρp = 0
TEMPEST q
0
1 2 3 4 5
b)
q
c)
FIG. 3: (a) Time evolution of the zonal-GAM potential shows GAM oscillation, collisionless damping, and residual for a large-aspect-ratio circular geometry with q = 2.2 and = 0.02 with three different velocity resolutions; (b) Comparison of simulation results with theory for GAM damping rate with three different velocity resolutions in finite banana orbit regime; (c) GAM damping rate vs q with finite banana orbit effect (yellow) and without finite banana orbit effect (black) from Sugama and Watanabe theory [14].
3.4 5D Drift waves and ITG modes
Conventional orderings are different for 4D neoclassical transport and 5D (ψ, θ, ζ, E0, µ) turbulence where ζ as the binormal coordinate. Our 5D gyrokinetic formulation yields a set of nonlinear electrostatic gyrokinetic equations that are valid for both neoclassical and turbulence simulations. In particular, the field solver for shear/zonal flow is different from that for turbulence due to: (1) the strong poloidal variation of the electrostatic potential in the divertor X-point geometry originating from different boundary conditions in the core, the SOL and private-flux regions; (2) additional terms are promoted by the edge ordering of a large background E × B flow uE vTi. 5D tempest uses field-aligned coordinates with 4th-order interpolation and ζ-index shifting for twist-shifted parallel boundary condition in the core inside the separatrix and radial difference in the usual flux coordinates to minimize the cell distortion in field-aligned coordinates due to magnetic shear. 5D TEMPEST shows
9
TH/P6-23
a good agreement for drift wave frequency between theory and simulations with 10% radial variation of ion density profile and the flat ion temperature profile. For a different simulation setup with 10% radial variation of ion temperature profile and the flat ion density profile, 5D Tempest also shows ITG mode. Benchmarks with theory and other codes are under progress.
4. Summary and conclusions
Recently developed full-f, 5D continuum edge-plasma code TEMPEST utilizes high-order spatial differencing and high-order finite-volume scheme for velocity space to accurately simulate particle convection and Coulomb collisions. TEMEST runs both in full divertor geometry for the important edge kinetic physics and in a circular geometry for benchmarking with analytical theories and other existing core gyrokinetic turbulence codes. TEMPEST demonstrates expected physics results in 3D, 4D and 5D verification tests. The further improvement and development of TEMPEST will yield a valuable predictive model for the edge pedestal. This work is focused on a fundamental understanding of relevant physics from first-principles theory and simulations and should greatly increase our confidence in predictions of ITER edge-plasma performance.
References
[1] Brizard A J and Hahm T S, submitted to Rev. Mod. Phys., (2006) [2] Qin H, Cohen R H, Nevins W M, and Xu X Q, Contrib. Plasma Phys., 46, (2006) 477 [3] Hahm T S, 1996 Phys. Plasmas 3, 4658 [4] Byrne S D and Hindmarsh A C, 1999 Int. J. High Perf. Comput. Appl., 13, 354 [5] Jiang G S and Shu C W, 1996 J. Comput. Phys. 126, 202 [6] Kupfer K, Harvey R W, Sauter O, et al., 1996 Phys. Plasmas 3, 3644 [7] Dorland W and Hammett G W, Phys. Fluids B 5, 812 (1993) [8] Dimits A, IAEA technical committee meeting 92 [9] Rognlien, T D and Cutler, Nuclear Fusion 20, 1003(1980) [10] Cohen, R H, Nuclear Fusion 19, 1295(1979) [11] Najmabadi,et al., Nuclear Fusion 24, 75(1984) [12] Rosenbluth, M N and Hinton F L,Phys. Rev. Lett., 80 (1998) 724 [13] Lebedev V B, Yushmanov P N, and Diamond P H, et al. Phys. Plasmas 3 (1996) 3023 [14] Sugama H and Watanable T H, submitted to J. Phys. Plasmas, (2006).
TH/P6-23
Continuum Edge Gyrokinetic Theory and Simulations1
X.Q. Xu 1), K. Bodi 2), J. Candy 3), B. I. Cohen 1), R. H. Cohen 1), P. Colella 4), M. R. Dorr 1), J. A. Hittinger 1), G. D. Kerbel 1), S. Krasheninnikov 2), W. M. Nevins 1), H. Qin 5), T. D. Rognlien 1), P. B. Snyder 3), M. V. Umansky 1), Z. Xiong 1)
1) Lawrence Livermore National Laboratory, Livermore, CA 94551 USA 2) University of California, San Diego, La Jolla, CA 92093 USA 3) General Atomics, San Diego, CA 92186 USA 4) Lawrence Berkeley National Laboratory, Berkeley, CA 94720 USA 5) Princeton Plasma Physics Laboratory, Princeton, NJ 08543 USA
e-mail contact of main author: [email protected]
Abstract. The following results are presented from the development and application of TEMPEST, a fully nonlinear (full-f) five dimensional (3d2v) gyrokinetic continuum edgeplasma code. (1) As a test of the interaction of collisions and parallel streaming, TEMPEST is compared with published analytic and numerical results for endloss of particles confined by combined electrostatic and magnetic wells. Good agreement is found over a wide range of collisionality, confining potential, and mirror ratio; and the required velocity space resolution is modest. (2) In a large-aspect-ratio circular geometry, excellent agreement is found for a neoclassical equilibrium with parallel ion flow in the banana regime with zero temperature gradient. (3) The four dimensional (2d2v) version of the code produces the first self-consistent simulation results of collisionless damping of geodesic acoustic modes and zonal flow (Ronsenbluth-Hinton residual) with Boltzmann electrons using a full-f code. In divertor geometry, it is found that the endloss of particles and energy induces pedestal-like density and temperature profiles inside the magnetic separatrix and parallel flow stronger than the core neoclassical predictions in the SOL. (4) Our 5D gyrokinetic formulation yields a set of nonlinear electrostatic gyrokinetic equations that are for both neoclassical and turbulence simulations.
1. Introduction
Understanding the structure of edge transport barrier in high-performance (H-mode) discharges requires a kinetic description of the plasmas because the radial width of the pedestal observed in experiments is comparable to the radial width of individual ion orbits (leading to large distortion of the local distribution function from a Maxwellian), while the ion and electron mean-free-paths are long compared to the connection length in the hot plasma at the top of the edge pedestal (violating the assumptions underlying a collisional fluid mode). Gyrokinetic formulation (2v) [1] is a reasonable approximation for edge plasmas because much pedestal physics is typically low frequency phenomenon compared to ion gyro-frequency. But existing gyrokinetic theories and codes do not apply to edge plasmas because they cannot treat fully nonlinear electromagnetic perturbations with multi-scale-length structures in space-time for full divertor geometry.
1This work was performed under the auspices of the U.S. Department of Energy by University of California Lawrence Livermore National Laboratory under contract No. W-7405-Eng-48 at LLNL, Grant No. DE-FG0204ER54739 at UCSD, Grants DE-FG03-95ER54309 at general Atomics, and DE-AC02-76CHO3073 at PPPL.
2
TH/P6-23
We report on the development and application of TEMPEST, a fully nonlinear (full-f) gyrokinetic code, to simulate H-mode edge plasmas. This 5-dimensional (ψ, θ, ζ, E0, µ) continuum code represents velocity space via a grid in equilibrium energy (E0) and magnetic moment (µ) variables, and configuration space via a grid in poloidal magnetic flux (ψ), poloidal angle (θ) and toroidal angle (ζ). The geometry is that of a fully diverted tokamak and so includes boundary conditions for both closed magnetic flux surfaces and open field lines. A set of gyrokinetic equations [2, 3] are discretized on a computational grid for both circular and divertor geometry with a magnetic separatrix. The equations are solved via a Method-of-Lines approach and an implicit backward-differencing scheme using a NewtonKrylov iteration advances the system in time [4]. The spatial derivatives are discretized with finite differences while a high-order finite volume method is used in velocity space (E0, µ). A fourth-order upwinding algorithm is used for parallel streaming; and a fifth-order WENO scheme [5] is used for particle cross-field drifts. Boundary conditions at conducting material surfaces are implemented on the plasma side of the sheath. The code includes kinetic or Boltzmann electrons. A nonlinear Fokker-Planck collision operator (CQL) from STELLA [6] has been extracted and integrated into the gyrokinetic simulation code using the same implicit Newton-Krylov solver. A new Fokker-Planck collision operator in (E0, µ) space is under development for improved accuracy and conservation properties. The gyrokinetic Poisson equation is solved using GMRES with a multi-grid preconditioner.
2. Basic Gyrokinetic Equation
A set of generalized gyrokinetic Vlasov-Maxwell equations in the gyrocenter coordinate system valid for edge-plasma conditions has been derived by the Lie transform perturbation method to the Poincar´e-Cartan-Einstein 1-form and the pullback transformation for the distribution function [2]. This formalism allows inclusion of nonlinear large-amplitude, time-dependent background electromagnetic fields in addition to small-amplitude, shortwavelength electromagnetic perturbations. As an example, the pullback transformation in the gyrokinetic Poisson equation is explicitly expressed in terms of moments of the gyrocenter distribution function, thus describing the important gyro-orbit squeezing effect due to the large electric field shearing in the edge and the full finite Larmour radius effect for short wavelength fluctuations. The familiar polarization drift density in the gyrocenter Poisson equation is replaced by a more general expression.
2.1 Fully Nonlinear Ion Gyrokinetic Equations
The ion gyro-kinetic equations presently implemented in TEMPEST for the time-dependent five-dimensional (5D) distribution functions are simplified from our recent new formulation [2] and earlier Hahm’s work [3]. Due to the large potential and its multiple spatial-time scale nature in the edge pedestal, in order to accurately simulate particle parallel streaming, the field Φ is split into two parts Φ = Φ0 + δφ: Φ0 is the large amplitude, slowly varying component; δφ is the small amplitude, rapidly varying component. Here E0 is defined as the total energy including Φ0, but not δφ. Then E0 is a constant of motion if δφ ∼ 0 and a coordinate aligned with the direction of phase-space flow. The kinetic equation for the gyrocenter distribution function Fα(x¯, µ¯, E¯0, t) in gyrocenter coordinates (x¯ = x − ρα, ρα = b × v/Ωcα),
3
TH/P6-23
“total energy” E¯0, and magnetic moment µ¯, has the form:
∂Fα + v¯d · ∂Fα + (v¯ α + vBanos)b · ∂Fα (1)
∂t
∂x¯⊥
∂x¯
+ q ∂ Φ0 + µ¯ ∂B − B v¯ q ∂ δφ − v · (q∇¯ δφ ) ∂Fα = C(F , F ),
∂t
∂t B∗ ∂s
d0
∂E0
αα
v¯d = qcBb∗ × q∇¯ Φ + µ¯∇¯ B + v¯2 MqBα∗c (∇¯ ×b). (2)
v¯d0 = qcBb∗ × q∇¯ Φ0 + µ¯∇¯ B + v¯2 MqBα∗c (∇¯ ×b). (3)
2
v¯ = ±
(E − µ¯B − q Φ ), v
=
µc (b
·
∇¯ ×b),
(4)
Mα 0
0
Banos q
B∗ ≡ B 1 + b · (v ∇¯ ×b) , Ω = qB , µ = Mαv⊥2 ,
(5)
α
Ωcα
cα Mαc
2B
δφ = Φ − Φ0 .
(6)
Here Zαe, Mα are the electric charge and mass of electrons (α = e), ions (α = i). µ is the guiding center magnetic moment. The left-hand side of Eq. (1) describes the particle motion in the electric field and magnetic field. Cα is the Coulomb collision operator. The over-bar is used for the gyrocenter variables and denotes the gyroangle averaging. The additional E0 × B flow terms due to the large amplitude and the slow variation Φ0 from the complete formulation [2] will be added.
2.2 Fully Nonlinear Gyrokinetic Poisson equation
The complete gyrokinetic Poisson equation has been derived [2], including orbit squeezing by large Er shearing and full FLR effect. To make it numerically attractive, two additional approximations are made here: (1). the spatial variation of the transverse µ¯ moments Mn(x¯) calculated from Fα(x¯, µ¯, E¯0, t) is much slower than that of potential in evaluation of the full FLR effect; (2). the total transverse distribution function is Maxwellian.
2.2.1 Fully Nonlinear Gyro-kinetic Poisson equation in the arbitrary wavelength regime
In the arbitrary wavelength regime, the self-consistent electrostatic potential is computed from the gyro-kinetic Poisson equation:
1 0 = −4πe α ZαNα(x, t) − ne(x, t) − α λ2Dα [Γ0(b) − 1] Φ, (7)
where Γ0(b) = I0(b)e−b, b = ρ2α∇2⊥/2, I0(b) is the usual zeroth-order modified Bessel function.
The ion gyroradius is ρα = 2T⊥α/Mα/Ωα, the ion gyrofrequency is Ωα = ZαeB/Mαc, and the ion Debye length is λ2Dα = T⊥α/4πNαZα2e2. Although Eq. (7) is similar to the usual gyrokinetic Poisson equation [3], there is an important distinction. Our gyro-kinetic Poisson equation is fully nonlinear and the gyrocenter center density Nα and perpendicular ion pressure p⊥α are calculated from the gyrocenter distribution function Fα(x¯, µ¯, E¯0, t).
2π N (x, t) ≡
B∗dv¯ dµ¯F ,
2π n (x, t) ≡
B∗dv dµf ,
(8)
α
Mα
α
e
me
e
4
TH/P6-23
p = πB dv dµ¯(v2 F ), T = p⊥α
(9)
⊥α
⊥α
⊥α Nα(x, t)
Here the dot product between the density and potential as well as the Debye shielding have been dropped for simplicity.
The first-order Pad´e approximation to Γ0 with Γ0 − 1 = b/(1 + b) is an excellent fit for 0 ≤ b ≤ 9, and is therefore valid well into the typical ion gyrokinetic regime as used previously
in gyro-kinetic or gyro-fluid simulations [7, 8]. For single ion species, substituting a simple functional transformation Φ = φL + [T⊥α/(NαZα2e)] [ZαNα(x, t) − ne(x, t)] and the Pad´e approximation into Eq. (7) yields
ρ2α ∇2 φ = − Tα 1 + ρ2α ∇2 ln T⊥α [Z N (x, t) − n (x, t)] .
(10)
2 ⊥L
Nα Zα2 e
2 ⊥ Nα
αα
e
where φL is obtained by the gyrokinetic Poisson solver.
2.2.2 Fully Nonlinear Gyro-kinetic Poisson equation in the long wavelength regime
In the long wavelength limit k⊥ρα 1, the self-consistent electric field is typically computed from the gyrokinetic Poisson equation for multiple species
ρ2α
2
α 2λ2Dα ∇⊥Φ +
ρ2α ∇⊥ ln Nα α 2λ2Dα
· ∇⊥Φ + ∇2Φ
= −4πe
Z N (x, t) − n (x, t) − ρ2α 1 ∇2 p . (11)
αα
e
α
α 2λ2Dα NαZαe ⊥ ⊥α
The ion gyroradius is ρα = 2T⊥α/Mα/Ωα, the ion gyrofrequency is Ωα = ZαeB/Mαc, and the ion Debye length is λ2Dα = T⊥α/4πnαZα2e2.
There are two important distinctions between Eq. (11) and the usual gyro-kinetic Poisson equation [3]. Our gyro-kinetic Poisson equation is fully nonlinear and the gyrocenter center density Nα and perpendicular ion pressure p⊥α are calculated from the gyrocenter distribution function Fα(x¯, µ¯, E¯0, t), as defined in Eqs. (8)-(9). The last term of Eq. (11) is the diamagnetic density from the long wavelength expansion of the gyroaveraged gyrocenter density Nα(x, t). Although the diamagnetic density is small compared to the ion gyrocenter density, it is of the same order as both the polarization density in high-beta plasmas and as the difference between ion and electron gyrocenter density! This equation is an extension of typical neoclassical electric field calculation including poloidal variation.
2.3 Boundary conditions
2.3.1 Radial boundary conditions
The radial Robin boundary conditions are used for Fα and potential δφ at the inner core surface ψ = ψc and the outer wall surface ψ = ψw. This is a generalization of the Dirichlet and Neumann boundary conditions. Since the gyrokinetic equation has just a first-order radial advection term, there is only one physical boundary condition, depending on the direction of convection, No boundary condition should be imposed for particles convecting out of the simulation domain and therefore a extrapolation is used at that boundary.
5
TH/P6-23
2.3.2 Poloidal boundary conditions
The boundary conditions in the θ direction for Fα and for Φ are the sheath boundary conditions at the divertor plates and periodic in the “core” (closed field-lines). Our present implementation for sheath boundary conditions is for the case of normal intersection of the flux surface with the walls, where the ions are not confined for the perfectly absorbing wall and the current through the sheath is zero with no biasing; and there is an energetic group of impinging electrons that can escape the sheath potential and reach the wall with the energy E0 + eδφsh − µB > 0. Here Φsh = Φ0sh + δφsh is the sheath potential.
i. Sheath boundary conditions for potential
If the gyrokinetic ion and fluid electron model are used, the sheath potential is determined by the quasi-neutrality condition:
Φ = Te ln
4Γi,sh
,Γ
2πB ∞
=
dE
(E0−qΦ0sh)/B
dµ v
F σ.(12)
sh
e
i,sh
ne,shζ 8Te,sh/πme
Mα2 qΦ0sh
0 0
|v | i
The σ = ± represents the plus and minus sheet of parallel velocity with Fiσ = 0 for only incoming sheet. Here it is assumed that impinging electrons have a Maxwellian distribution.
The factor ζ ≡ 1/(1 + τp/τe) includes the correction of electron long mean-free paths physics. τp is long mean-free path confinement time and τc is the confinement time for the collisional sheath limited case. ζ ≡ 1 if the electrons are in the short mean-free paths regime.
If both electron and ion are kinetic, the sheath potential is determined:
Γ
2πB ∞
=
dE
(E0−qΦ0sh)/B
dµ v
F σ,
i,sh
Mα2 qΦ0sh
0 0
|v | i
Γ
2πB ∞
=
dE
(E0−eδφsh)/B
dµ v
F σ,
e,sh
m2e eδφsh
0 0
|v | e
Γi,sh = Γe,sh.
(13) (14)
ii. Sheath boundary conditions for distribution functions If the gyrokinetic ion and fluid electron model is used, the ion distribution function is:
Fα(ψ, θ, E0, µ) = Fα(ψ, θ, E0, µ), v ≥ 0
(15)
0,
v ≤0
A convention regarding the sign of the parallel velocity is that it is positive when there is a positive projection on the θ axis. Here positive θ axis is pointing to the plate/wall.
If both electron and ion are kinetic, the electron distribution function is:
fe(ψ, θ, E0, µ, σ = 1) = fe(ψ, θ, E0, µ, σ = 1),
(16)
fe(ψ, θ, E0, µ, σ = −1) = fe(ψ, θ, E0, µ, σ = 1), |v | ≤ vSH
(17)
0,
|v | ≥ vSH
Here vsh = 2eΦsh/me.
6
TH/P6-23
3. TEMPEST simulation Results
TEMPEST has been developed in the modern framework for either circular or divertor geometry. With importance of benchmark in mind, both full-f or δf options are available. TEMPEST is runnable as (1) 3D for the SOL endloss with F (θ, E0, µ), (2) 4D for transport F (ψ, θ, E0, µ), and (3) 5D for turbulence F (ψ, θ, ζ, E0, µ). The different aspects of the 3D, 4D and 5D TEMPEST have been verified on various known physics problems, such as (1) collisional scattering into a velocity-space loss cone, (2). neoclassical flow and transport, (3). electric field generation and geodesic acoustic mode damping, (4). drift waves and ion temperature gradient (ITG) modes.
3.1 3D Pastukhov collisional endloss
As a test of collisional velocity-space transport and parallel streaming, 3D TEMPEST (1d2v) simulation results are compared with published analytic and numerical results as shown in Fig. 1 for the endloss of particles confined by combined electrostatic and magnetic wells [9, 10, 11]. Here electrostatic and magnetic field are uniform with abrupt confining magnetic field and potential barriers at wall. Good agreement is found over a wide range of collisionality, confining potential, and mirror ratio; and the required velocity-space resolution is modest. In these simulations, the linearized CQL package is used.
a)
b)
c)
FIG. 1: Collisional endloss (“Pastukohov”) test cases: (a) confinement time versus density; (b) confinement time versus potential eφ/Te at low collisionality; (c) confinement versus mirror ratio at low collisionality.
3.2 4D Neoclassical tests
If using a shifted Maxwellian distribution that analytically satisfies Eq. (1) as an initial condition, then TEMPEST should preserve the solution without any significant change (within our finite difference truncation accuracy) after running some time steps. We tested such a case using the following simulation parameters: inverse aspect ratio = a/R0 = 0.03, the major radius R0 = 17.1 meter, toroidal magnetic field Bt = 1.5T , and poloidal
7
TH/P6-23
magnetic field Bp = 0.2T . The ion density and temperature profiles used are ni(ψ) = Nix exp(− ln(Nix/Nio)ψ/Lψ), Ti(ψ) = Tix exp(− ln(Tix/Tio)ψ/Lψ). with Tix = 3keV, Tio = 0.95Tix, Nix = 1 × 1020m−3, and Nio = 0.95Nix. The mesh resolutions are nψ = 30, nθ = 50, nE = 60, and nµ = 30. In the simulations Φ is set to be zero for simplicity. As shown in Fig. 2a), the simulation results remain in good agreement with theoretical prediction even after 10000 time steps ( 50 thermal ion transit time). The solid line in plot of q i comes from theoretical prediction for a shifted Maxwellian distribution, q α = 2.5NαU αTα, and U α = − ΩIα MTαα ∂ l∂nψNα .
In divertor geometry, for a given particle and heat sources on the inner core boundary surface, because the parallel transport is much larger than neoclassical radial transport, it is found that the endloss of particles and energy (discussed in section 2.3) induces (1) pedestal-like density and temperature profiles inside the magnetic separatrix; (2) parallel flow stronger than core neoclassical predictions in the SOL as indicated in Fig. 2 (b); (3) a symmetry point is developed for the parallel heat flux on the top of the machine as expected as shown in Fig. 2(c).
1.5×10 24
1.0×10 24
0.5×1024
10000 8000 6000 4000
theory
time 3 time 2 time 1 time 0
q (R,Z)(eV/m2/s) ||
2.5 2.0 Z(m) 1.5 1.0
separatrix
0 0.0
0.1 0.2 0.3 0.4 Radial Position (m)
a)
2000
0 0.0 0.02 0.04 0.06 0.08 Radial position (m)
b)
0.5
1.0 1.2 1.4 1.6 1.8 2.0 2.2 R(m)
c)
FIG. 2: (a) Comparison between simulation results with theory for a collisionless case with ∇Ti = 0 and zero finite banana orbit width. Flux surface averaged parallel heat flux q i . The solid lines comes theory and other lines from different time; (b) Comparison of simulation results with theory for flux surface averaged parallel flow velocity U i in the banana regime with ν∗i 0.02, ∇Ti = 0 and finite banana orbit width in X-point divertor geometry. The average is done by integration along the field-line from inner plate to outer plate in the SOL. (c) The contours of parallel heat flux q (R, Z) in the divertor geometry.
3.3 4D Geodesic-Acoustic Modes
The Geodesic-Acoustic Mode (GAM) is an asymmetric mode, which involves parallel ion dynamics, cross field drifts, and acceleration. Earlier GAM theory and simulations focused on large aspect ratio and small orbit [12, 13]. Recently Sugama and Watanabe find that the damping rate is sensitive to k⊥ρi at large q due to the effect of large banana orbit [14]. In our 4D GAM simulations, the charge is radially separated by an initial sinusoidal perturbation on
8
TH/P6-23
ion density. The electron model is Boltzmann ne = ni(ψ, θ, t = 0) exp(eφ/Te)/ exp(eφ/Te) , where represents the flux surface average. This choice of coefficient for Boltzmann electron model means that there is no cross field electron transport. Both radial and poloidal boundary conditions are periodic. The first full-f, self-consistent simulation results of collisionless damping of geodesic acoustic modes and zonal flow are plotted in Fig. 3. Good agreement is shown between theory [14] and simulations for the frequency of GAMs in Fig. 3a) and damping rate in Fig. 3b). The 30% difference between theory and simulation is attributed to the fact that the theory uses an asympototic 1/q2 expansion for large q, while in our simulation q=2.2 is moderate. The large effect of the orbit size on the GAM damping rate is illustrated in Fig. 3c). For the same parameters, the damping rate is almost zero if the finite banana orbit effect is ignored.
φ(t) / φ(t=0) -γ / (vt/qR)
1.0
Velocity mesh
(nE,nµ) = (25,50) (nE,nµ) = (45,50)
γTEMPEST = 0.072(vtd/R)
0.5
(nE,nµ) = (91,101)
Numerical residual
GAM damping rate 0.15
Sugama & Watanabe 06 0.10 kρp = 0.1
γSugama-Watanabe = 0.052(vtd/R) 0
Spatial mesh
(nψ,nθ) = (30,50)
-0.5
a/R0 = 0.02
q = 2.2
0
25
50
75
Time (R0/vtd)
a)
0.05
kρp = 0
TEMPEST q
0
1 2 3 4 5
b)
q
c)
FIG. 3: (a) Time evolution of the zonal-GAM potential shows GAM oscillation, collisionless damping, and residual for a large-aspect-ratio circular geometry with q = 2.2 and = 0.02 with three different velocity resolutions; (b) Comparison of simulation results with theory for GAM damping rate with three different velocity resolutions in finite banana orbit regime; (c) GAM damping rate vs q with finite banana orbit effect (yellow) and without finite banana orbit effect (black) from Sugama and Watanabe theory [14].
3.4 5D Drift waves and ITG modes
Conventional orderings are different for 4D neoclassical transport and 5D (ψ, θ, ζ, E0, µ) turbulence where ζ as the binormal coordinate. Our 5D gyrokinetic formulation yields a set of nonlinear electrostatic gyrokinetic equations that are valid for both neoclassical and turbulence simulations. In particular, the field solver for shear/zonal flow is different from that for turbulence due to: (1) the strong poloidal variation of the electrostatic potential in the divertor X-point geometry originating from different boundary conditions in the core, the SOL and private-flux regions; (2) additional terms are promoted by the edge ordering of a large background E × B flow uE vTi. 5D tempest uses field-aligned coordinates with 4th-order interpolation and ζ-index shifting for twist-shifted parallel boundary condition in the core inside the separatrix and radial difference in the usual flux coordinates to minimize the cell distortion in field-aligned coordinates due to magnetic shear. 5D TEMPEST shows
9
TH/P6-23
a good agreement for drift wave frequency between theory and simulations with 10% radial variation of ion density profile and the flat ion temperature profile. For a different simulation setup with 10% radial variation of ion temperature profile and the flat ion density profile, 5D Tempest also shows ITG mode. Benchmarks with theory and other codes are under progress.
4. Summary and conclusions
Recently developed full-f, 5D continuum edge-plasma code TEMPEST utilizes high-order spatial differencing and high-order finite-volume scheme for velocity space to accurately simulate particle convection and Coulomb collisions. TEMEST runs both in full divertor geometry for the important edge kinetic physics and in a circular geometry for benchmarking with analytical theories and other existing core gyrokinetic turbulence codes. TEMPEST demonstrates expected physics results in 3D, 4D and 5D verification tests. The further improvement and development of TEMPEST will yield a valuable predictive model for the edge pedestal. This work is focused on a fundamental understanding of relevant physics from first-principles theory and simulations and should greatly increase our confidence in predictions of ITER edge-plasma performance.
References
[1] Brizard A J and Hahm T S, submitted to Rev. Mod. Phys., (2006) [2] Qin H, Cohen R H, Nevins W M, and Xu X Q, Contrib. Plasma Phys., 46, (2006) 477 [3] Hahm T S, 1996 Phys. Plasmas 3, 4658 [4] Byrne S D and Hindmarsh A C, 1999 Int. J. High Perf. Comput. Appl., 13, 354 [5] Jiang G S and Shu C W, 1996 J. Comput. Phys. 126, 202 [6] Kupfer K, Harvey R W, Sauter O, et al., 1996 Phys. Plasmas 3, 3644 [7] Dorland W and Hammett G W, Phys. Fluids B 5, 812 (1993) [8] Dimits A, IAEA technical committee meeting 92 [9] Rognlien, T D and Cutler, Nuclear Fusion 20, 1003(1980) [10] Cohen, R H, Nuclear Fusion 19, 1295(1979) [11] Najmabadi,et al., Nuclear Fusion 24, 75(1984) [12] Rosenbluth, M N and Hinton F L,Phys. Rev. Lett., 80 (1998) 724 [13] Lebedev V B, Yushmanov P N, and Diamond P H, et al. Phys. Plasmas 3 (1996) 3023 [14] Sugama H and Watanable T H, submitted to J. Phys. Plasmas, (2006).