# 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) ﬁve 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 conﬁned by combined electrostatic and magnetic wells. Good agreement is found over a wide range of collisionality, conﬁning 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 ﬂow in the banana regime with zero temperature gradient. (3) The four dimensional (2d2v) version of the code produces the ﬁrst self-consistent simulation results of collisionless damping of geodesic acoustic modes and zonal ﬂow (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 proﬁles inside the magnetic separatrix and parallel ﬂow 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 ﬂuid 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 conﬁguration space via a grid in poloidal magnetic ﬂux (ψ), poloidal angle (θ) and toroidal angle (ζ). The geometry is that of a fully diverted tokamak and so includes boundary conditions for both closed magnetic ﬂux surfaces and open ﬁeld 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-diﬀerencing scheme using a NewtonKrylov iteration advances the system in time [4]. The spatial derivatives are discretized with ﬁnite diﬀerences while a high-order ﬁnite volume method is used in velocity space (E0, µ). A fourth-order upwinding algorithm is used for parallel streaming; and a ﬁfth-order WENO scheme [5] is used for particle cross-ﬁeld 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 ﬁelds 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 eﬀect due to the large electric ﬁeld shearing in the edge and the full ﬁnite Larmour radius eﬀect for short wavelength ﬂuctuations. 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 ﬁve-dimensional (5D) distribution functions are simpliﬁed 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 ﬁeld Φ is split into two parts Φ = Φ0 + δφ: Φ0 is the large amplitude, slowly varying component; δφ is the small amplitude, rapidly varying component. Here E0 is deﬁned 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 ﬂow. 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 ﬁeld and magnetic ﬁeld. Cα is the Coulomb collision operator. The over-bar is used for the gyrocenter variables and denotes the gyroangle averaging. The additional E0 × B ﬂow 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 eﬀect. 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 eﬀect; (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 modiﬁed 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 ﬁrst-order Pad´e approximation to Γ0 with Γ0 − 1 = b/(1 + b) is an excellent ﬁt for 0 ≤ b ≤ 9, and is therefore valid well into the typical ion gyrokinetic regime as used previously

in gyro-kinetic or gyro-ﬂuid 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 ﬁeld 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 deﬁned 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 diﬀerence between ion and electron gyrocenter density! This equation is an extension of typical neoclassical electric ﬁeld 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 ﬁrst-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 ﬁeld-lines). Our present implementation for sheath boundary conditions is for the case of normal intersection of the ﬂux surface with the walls, where the ions are not conﬁned 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 ﬂuid 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 conﬁnement time and τc is the conﬁnement 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 ﬂuid 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 diﬀerent aspects of the 3D, 4D and 5D TEMPEST have been veriﬁed on various known physics problems, such as (1) collisional scattering into a velocity-space loss cone, (2). neoclassical ﬂow and transport, (3). electric ﬁeld 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 conﬁned by combined electrostatic and magnetic wells [9, 10, 11]. Here electrostatic and magnetic ﬁeld are uniform with abrupt conﬁning magnetic ﬁeld and potential barriers at wall. Good agreement is found over a wide range of collisionality, conﬁning 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) conﬁnement time versus density; (b) conﬁnement time versus potential eφ/Te at low collisionality; (c) conﬁnement versus mirror ratio at low collisionality.

3.2 4D Neoclassical tests

If using a shifted Maxwellian distribution that analytically satisﬁes Eq. (1) as an initial condition, then TEMPEST should preserve the solution without any signiﬁcant change (within our ﬁnite diﬀerence 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 ﬁeld Bt = 1.5T , and poloidal

7

TH/P6-23

magnetic ﬁeld Bp = 0.2T . The ion density and temperature proﬁles 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 proﬁles inside the magnetic separatrix; (2) parallel ﬂow stronger than core neoclassical predictions in the SOL as indicated in Fig. 2 (b); (3) a symmetry point is developed for the parallel heat ﬂux on the top of the machine as expected as shown in Fig. 2(c).

1.5×10 24

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) ﬁve 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 conﬁned by combined electrostatic and magnetic wells. Good agreement is found over a wide range of collisionality, conﬁning 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 ﬂow in the banana regime with zero temperature gradient. (3) The four dimensional (2d2v) version of the code produces the ﬁrst self-consistent simulation results of collisionless damping of geodesic acoustic modes and zonal ﬂow (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 proﬁles inside the magnetic separatrix and parallel ﬂow 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 ﬂuid 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 conﬁguration space via a grid in poloidal magnetic ﬂux (ψ), poloidal angle (θ) and toroidal angle (ζ). The geometry is that of a fully diverted tokamak and so includes boundary conditions for both closed magnetic ﬂux surfaces and open ﬁeld 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-diﬀerencing scheme using a NewtonKrylov iteration advances the system in time [4]. The spatial derivatives are discretized with ﬁnite diﬀerences while a high-order ﬁnite volume method is used in velocity space (E0, µ). A fourth-order upwinding algorithm is used for parallel streaming; and a ﬁfth-order WENO scheme [5] is used for particle cross-ﬁeld 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 ﬁelds 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 eﬀect due to the large electric ﬁeld shearing in the edge and the full ﬁnite Larmour radius eﬀect for short wavelength ﬂuctuations. 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 ﬁve-dimensional (5D) distribution functions are simpliﬁed 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 ﬁeld Φ is split into two parts Φ = Φ0 + δφ: Φ0 is the large amplitude, slowly varying component; δφ is the small amplitude, rapidly varying component. Here E0 is deﬁned 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 ﬂow. 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 ﬁeld and magnetic ﬁeld. Cα is the Coulomb collision operator. The over-bar is used for the gyrocenter variables and denotes the gyroangle averaging. The additional E0 × B ﬂow 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 eﬀect. 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 eﬀect; (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 modiﬁed 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 ﬁrst-order Pad´e approximation to Γ0 with Γ0 − 1 = b/(1 + b) is an excellent ﬁt for 0 ≤ b ≤ 9, and is therefore valid well into the typical ion gyrokinetic regime as used previously

in gyro-kinetic or gyro-ﬂuid 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 ﬁeld 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 deﬁned 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 diﬀerence between ion and electron gyrocenter density! This equation is an extension of typical neoclassical electric ﬁeld 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 ﬁrst-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 ﬁeld-lines). Our present implementation for sheath boundary conditions is for the case of normal intersection of the ﬂux surface with the walls, where the ions are not conﬁned 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 ﬂuid 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 conﬁnement time and τc is the conﬁnement 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 ﬂuid 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 diﬀerent aspects of the 3D, 4D and 5D TEMPEST have been veriﬁed on various known physics problems, such as (1) collisional scattering into a velocity-space loss cone, (2). neoclassical ﬂow and transport, (3). electric ﬁeld 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 conﬁned by combined electrostatic and magnetic wells [9, 10, 11]. Here electrostatic and magnetic ﬁeld are uniform with abrupt conﬁning magnetic ﬁeld and potential barriers at wall. Good agreement is found over a wide range of collisionality, conﬁning 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) conﬁnement time versus density; (b) conﬁnement time versus potential eφ/Te at low collisionality; (c) conﬁnement versus mirror ratio at low collisionality.

3.2 4D Neoclassical tests

If using a shifted Maxwellian distribution that analytically satisﬁes Eq. (1) as an initial condition, then TEMPEST should preserve the solution without any signiﬁcant change (within our ﬁnite diﬀerence 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 ﬁeld Bt = 1.5T , and poloidal

7

TH/P6-23

magnetic ﬁeld Bp = 0.2T . The ion density and temperature proﬁles 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 proﬁles inside the magnetic separatrix; (2) parallel ﬂow stronger than core neoclassical predictions in the SOL as indicated in Fig. 2 (b); (3) a symmetry point is developed for the parallel heat ﬂux 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 ﬁnite banana orbit width. Flux surface averaged parallel heat ﬂux q i . The solid lines comes theory and other lines from diﬀerent time; (b) Comparison of simulation results with theory for ﬂux surface averaged parallel ﬂow velocity U i in the banana regime with ν∗i 0.02, ∇Ti = 0 and ﬁnite banana orbit width in X-point divertor geometry. The average is done by integration along the ﬁeld-line from inner plate to outer plate in the SOL. (c) The contours of parallel heat ﬂux 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 ﬁeld drifts, and acceleration. Earlier GAM theory and simulations focused on large aspect ratio and small orbit [12, 13]. Recently Sugama and Watanabe ﬁnd that the damping rate is sensitive to k⊥ρi at large q due to the eﬀect 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 ﬂux surface average. This choice of coeﬃcient for Boltzmann electron model means that there is no cross ﬁeld electron transport. Both radial and poloidal boundary conditions are periodic. The ﬁrst full-f, self-consistent simulation results of collisionless damping of geodesic acoustic modes and zonal ﬂow 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% diﬀerence 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 eﬀect 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 ﬁnite banana orbit eﬀect 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 diﬀerent velocity resolutions; (b) Comparison of simulation results with theory for GAM damping rate with three diﬀerent velocity resolutions in ﬁnite banana orbit regime; (c) GAM damping rate vs q with ﬁnite banana orbit eﬀect (yellow) and without ﬁnite banana orbit eﬀect (black) from Sugama and Watanabe theory [14].

3.4 5D Drift waves and ITG modes

Conventional orderings are diﬀerent 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 ﬁeld solver for shear/zonal ﬂow is diﬀerent from that for turbulence due to: (1) the strong poloidal variation of the electrostatic potential in the divertor X-point geometry originating from diﬀerent boundary conditions in the core, the SOL and private-ﬂux regions; (2) additional terms are promoted by the edge ordering of a large background E × B ﬂow uE vTi. 5D tempest uses ﬁeld-aligned coordinates with 4th-order interpolation and ζ-index shifting for twist-shifted parallel boundary condition in the core inside the separatrix and radial diﬀerence in the usual ﬂux coordinates to minimize the cell distortion in ﬁeld-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 proﬁle and the ﬂat ion temperature proﬁle. For a diﬀerent simulation setup with 10% radial variation of ion temperature proﬁle and the ﬂat ion density proﬁle, 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 diﬀerencing and high-order ﬁnite-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 veriﬁcation 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 ﬁrst-principles theory and simulations and should greatly increase our conﬁdence 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).