Brought to you by:
Paper

Local and global Fokker–Planck neoclassical calculations showing flow and bootstrap current modification in a pedestal

and

Published 8 October 2012 © 2012 IOP Publishing Ltd
, , Citation Matt Landreman and Darin R Ernst 2012 Plasma Phys. Control. Fusion 54 115006DOI 10.1088/0741-3335/54/11/115006

0741-3335/54/11/115006

Abstract

In transport barriers, particularly H-mode edge pedestals, radial scale lengths can become comparable to the ion orbit width, causing neoclassical physics to become radially nonlocal. In this work, the resulting changes to neoclassical flow and current are examined both analytically and numerically. Steep density gradients are considered, with scale lengths comparable to the poloidal ion gyroradius, together with strong radial electric fields sufficient to electrostatically confine the ions. Attention is restricted to relatively weak ion temperature gradients (but permitting arbitrary electron temperature gradients), since in this limit a δf (small departures from a Maxwellian distribution) rather than full-f approach is justified. This assumption is in fact consistent with measured inter-ELM H-Mode edge pedestal density and ion temperature profiles in many present experiments, and is expected to be increasingly valid in future lower collisionality experiments. In the numerical analysis, the distribution function and Rosenbluth potentials are solved for simultaneously, allowing use of the exact field term in the linearized Fokker–Planck–Landau collision operator. In the pedestal, the parallel and poloidal flows are found to deviate strongly from the best available conventional neoclassical prediction, with large poloidal variation of a different form than in the local theory. These predicted effects may be observable experimentally. In the local limit, the Sauter bootstrap current formulae appear accurate at low collisionality, but they can overestimate the bootstrap current in the plateau regime. In the pedestal ordering, ion contributions to the bootstrap and Pfirsch–Schlüter currents are also modified.

Export citation and abstractBibTeXRIS

1. Introduction

Neoclassical effects in a plasma—the flows, fluxes and currents determined by collisions in a toroidal equilibrium in the absence of turbulence—set a minimum level of radial transport [1, 2]. In transport barriers—the pedestal at the edge of an H-mode or internal transport barriers—neoclassical effects are particularly important for several reasons. First, the pressure gradient driven flows and bootstrap current (thought to be determined or at least strongly influenced by neoclassical physics even in the presence of turbulence) become large due to the small radial scale lengths. Second, turbulent radial transport is reduced, so neoclassical radial transport becomes more relevant. Both the flows and bootstrap current will affect the global stability of the transport barrier region. For example, to predict whether given plasma profiles are stable to edge localized modes (ELMs) and to predict the nature of such ELMs [3], accurate calculation of the bootstrap current is essential.

However, conventional neoclassical calculations are not formally valid in the pedestal. The reason is that in conventional neoclassical calculations, the main ion distribution function fi is expanded in an asymptotic series [1, 2] fi = fMi + f1 + ··· with f1/fMi ∼ ρθ/r ≪ 1, where fMi is a Maxwellian, ρθ = (B/Bθ)vi/Ω is the poloidal ion gyroradius, B = |B| is the magnitude of the magnetic field, Bθ is the poloidal magnetic field, r ∼ |∇ ln pi|−1 ∼ |∇ ln Ti|−1 is the scale length of ion pressure pi or ion temperature Ti, is the ion thermal speed, Ω = ZeB/(mic) is the gyrofrequency, Z is the ion charge in units of the proton charge e, mi is the ion mass and c is the speed of light. In the pedestal, r is observed to be comparable to ρθ in present experiments. In particular, the density gradient scale length is generally comparable to ρθ. (For this discussion it does not matter whether r actually scales with ρθ.) The first two terms in the asymptotic series fMi and f1 are then of comparable magnitude, so the asymptotic approach breaks down. In the conventional case, the orbit width (∼ρθ) is thin compared with the equilibrium profiles, so neoclassical effects are radially local: the flows on a given flux surface depend only on the physical quantities and their radial gradients at that surface. However, in a transport barrier where the ion orbit width is not small relative to the equilibrium scales, the ions will sample a range of densities and temperatures during their orbits. Accordingly, ion flows on a given flux surface are influenced by equilibrium parameters from neighboring flux surfaces that lie roughly within a poloidal gyroradius. Thus a radially global (i.e. nonlocal) calculation is required for the ion physics. A nonlocal calculation is unnecessary for electrons since their orbit widths are times smaller than ion orbit widths, but the electron distribution is nonetheless modified due to collisions with the modified ion distribution [4].

In the conventional local theory, a natural scale separation exists between flows within a flux surface, which are first order in the ρθ/r ≪ 1 expansion, and radial transport fluxes, which are second order in this expansion. This scale separation at least partially breaks down in the pedestal, and radial transport fluxes compete with flux-surface flows, even within a purely neoclassical framework. Our work includes this important effect, which strongly impacts the resulting flux-surface flows.

It is harder experimentally to measure the local bootstrap current density than to measure the plasma flow. Since the current is just the difference in ion and electron flows, validation of neoclassical flow calculations would give confidence in bootstrap current predictions. Impurity and main ion flows have been measured and compared with neoclassical predictions in several experiments, with mixed results [512]. Neoclassical theory makes an absolute prediction for the poloidal flow but not the toroidal or parallel flows, since the latter are a function of dΦ0/dψ, and this radial electric field cannot be determined within the lowest-order axisymmetric theory [13]. (Here, 2πψ is the poloidal flux.) The poloidal flows are largest in the steep-gradient transport barrier regions, yet these are precisely the regions in which the theory breaks down. For this reason as well, an improved nonlocal calculation of flows is sought to compare with measurements.

Even in the limit of small collisionality, the form of the collision operator is crucial for determining the neoclassical flows, fluxes, and current. The collision operator rigorously derived from first principles is the Fokker–Planck–Landau operator [14]. In much analytic and numerical work, however, simpler 'model' collision operators are used instead [1, 2, 1520]. Model operators generally yield somewhat different results for all neoclassical quantities [21, 22], so in the following work the exact linearized Fokker–Planck–Landau operator is used. At the same time, it should be remembered that even the exact Fokker–Planck operator is only correct to O(1/ln Λ) where lnΛ is the Coulomb logarithm.

Calculations of neoclassical quantities at realistic aspect ratio and with a realistic treatment of collisions require a numerical treatment. Local neoclassical computations with complete Fokker–Planck–Landau collisions were described by Sauter et al in [2325] and later extended to stellarator geometry in the code NEO described in [26, 27]. More recently, Fokker–Planck–Landau collisions have been implemented in other codes [22], including a second code called NEO [21] (unrelated to [26, 27]), and in [28]. All of these codes are radially local.

In recent years, a number of numerical efforts have been undertaken to compute nonlocal neoclassical effects in transport barriers. Most of these efforts have used the particle-in-cell (PIC) approach [2937]. PIC and continuum codes have differing treatments of collisions and boundary conditions, and face different numerical resolution challenges, so it is good practice to develop both approaches to verify they yield the same physical results. Some investigations of neoclassical effects have been begun in global continuum codes [3840], but these codes use approximate collision models and are ultimately designed for turbulence studies, and very different algorithms have been used than the ones we use here.

In this work, we present a new approach to computing global neoclassical effects. A continuum (Eulerian) framework is used, including the exact linearized Fokker–Planck–Landau collision operator. Our approach includes a general prescription for extending a local neoclassical code to incorporate nonlocal effects in a numerically efficient manner, by making such a local calculation the inner step of an iteration loop. Several terms related to the electric field must first be added to the local code, since in the pedestal these terms cannot be neglected.

Our approach is not completely general, for while we allow the ion density scale length rn to be ∼ρθ, we require the ion temperature scale length be >ρθ. The electron temperature scale length may be either ∼ρθ or >ρθ. The ordering is satisfied in the pedestal on many present tokamaks, including DIII-D, JET, ASDEX-U, NSTX and MAST [4151], when inter-ELM or non-ELMy (without RMP) profiles are carefully examined, though exceptions exist such as I-mode and EDA H-Mode in Alcator C-Mod [52]. Both entropy considerations [53] and data [45, 47] suggest resists becoming as small as ρθ when collisionality is low, while rn and are not similarly constrained. This suggests that future experiments with higher pedestal temperatures and lower collisionalities may increasingly satisfy . The entropy argument is based on a more careful version of the asymptotic analysis above, showing that while would require fi to depart strongly from a Maxwellian flux function, the same need not be true if or rn ∼ ρθ. Collisionless orbits radially average the ion temperature within a poloidal gyroradius, preventing strong ion temperature variation, while the density is not similarly averaged due to electrostatic confinement. The strong temperature gradient case for general collisionality requires use of the full nonlinear Fokker–Planck–Landau collision operator, giving rise to a kinetic equation that is nonlinear in fi. In the weak-T'i case we consider, we will show it is still appropriate to expand about a Maxwellian flux function fMi and use the linearized collision operator. The E × B-drift nonlinearity associated with the poloidal electric field also becomes negligible, making the collisionless part of the kinetic equation linear in δf = fi − fMi. The full-f strong-T'i case will be considered in future work. Any more general full-f nonlinear code must be able to accurately reproduce the weak-T'i limit, and since expansion about a Maxwellian is useful both numerically and analytically, it is worth understanding this limit in detail.

Another simplification in this work is that we assume Bθ ≪ B, separating ρθ from the gyroradius ρ = vi/Ω scale. Without this approximation, the desired ordering r ∼ ρθ would then imply the equilibrium varies on the ρ scale, so a drift-kinetic description would not be possible. This is not a serious limitation and is well satisfied for the edge region, which is characterized by large safety factors.

The primary finding in this work is that the ion flow is significantly altered in magnitude and direction relative to the prediction of local theory, and in particular, the flow's poloidal variation is qualitatively different. The poloidal variation of the flow is effectively determined by the requirement that the total flow be divergence free. In conventional theory, the total flow is approximately given by a sum of parallel, diamagnetic and leading-order E × B components, implying the poloidal flow must vary on a flux surface as Bθ. However, in the pedestal, two other contributions to the flow divergence grow to become leading-order terms: the E × B flow of the poloidally varying part of the density, and the radial variation in particle flux. As a result, the coefficients that multiply the ion temperature gradient in the parallel and poloidal flow are no longer equal, the poloidal flow no longer varies as Bθ, and the flow may change magnitude and sign relative to the local prediction. These effects may be important to consider in any comparison between experimental flow measurements in the pedestal and neoclassical theory [11, 12]. We present the details of one calculation of these effects at experimentally relevant aspect ratio and collisionality, considering a single ion species.

In the following section we review the relevant aspects of local neoclassical theory. At the core of our global solver is a local solver, so in the next section we discuss in detail the local solver used. New comparisons to reduced analytic models are presented. Section 4 then discusses a δf formulation for the global neoclassical problem in a transport barrier with a strong radial density gradient. Changes to the structure of the flow are discussed in section 5, and changes to the Pfirsch–Schlüter and bootstrap currents are calculated in section 6. Even in the δf formulation, the kinetic equation is challenging to solve by direct numerical methods, so section 7 introduces the operator-splitting initial-value-problem approach which reduces the dimension of the numerical problem to solve. Section 8 discusses the need for a sink term in the model and describes the sinks used. Results are presented in section 9, and we conclude in section 10.

2. Definitions and local theory

In the local case, the ion distribution function fi is approximately a Maxwellian with constant density ni(ψ) and temperature Ti(ψ) on each flux surface: fMi = ni[mi/(2πTi)]3/2exp(−miv2/[2Ti]). To next order, fi = fMi − ZeΦ1fMi/Ti + f1 where Φ1(ψ, θ) = Φ − Φ0, Φ is the electrostatic potential, Φ0(ψ) is the flux-surface average of Φ, and it can be shown |Φ1| ≪ |Φ0|. The distribution f1 is found by solving the following drift-kinetic equation:

Here, is the radial magnetic drift, and Ci is the linearized ion–ion Fokker–Planck–Landau collision operator. The derivatives in (1) are performed at fixed and total energy W0 = miv2/2 + ZeΦ0, so

where x = v/vi. We ignore the correction introduced by ion–electron collisions.

It is sometimes convenient to apply the identity vd · ∇ψ = (Iv||/Ω)∇||(v||/B), where I equals the toroidal field Bζ times the major radius R, to rewrite (1) as

Here, F = −(Iv||/Ω)∂fMi/∂ψ, and

Only a Ti gradient can drive g, not gradients in ni or Φ0. This result follows from Ci{v||fMi} = 0, so the dpi/dψ and dΦ0/dψ terms in (2) disappear entirely from Ci{F} and from (3).

Once fi is found, the two moments of greatest interest are the radial heat flux

and the parallel flow

Here, kq and k|| are dimensionless coefficients defined by the right equalities in (5)–(6), epsilon is the inverse aspect ratio, and brackets denote a flux-surface average:

for any quantity A where

dℓθ is the poloidal length element, V' = dV/dψ, and 2πV(ψ) is the volume enclosed by the flux surface. Also, is the ion–ion collision frequency with νi the Braginskii ion collision frequency. The definition for kq in (5) turns out to be convenient as kq then has a finite limit as epsilon → 0 and collisionality → 0.

It can be shown using the following argument that the parallel flow must have the form (6) with k|| constant on a flux surface. First, apply the operation to (3). (Here, σ = sgn(v||).) This operation annihilates the linearized collision operator terms by particle conservation. Pulling ∇|| in front of the velocity integrals, the boundary term from the upper limit of the dμ integral vanishes in the σ sum, leaving B||(∫d3v v||g/B) = 0, and so ∫d3v v||g = niXB where X is constant on a flux surface. Then applying to (4), and noting the last term in (2) vanishes in the v integral, the flow must have the form

Recalling g ∝ dTi/dψ, and normalizing X by convenient constants, then the form (6) results with k|| constant on a flux surface. The form (9) can also be understood from a fluid perspective, as follows. First, the leading-order perpendicular flow is V = cB−2(dΦ0/dψ + [Zeni]−1dpi/dψ)B × ∇ψ. Together with the mass continuity relation and B × ∇ψ = IB − R2B2ζ for toroidal angle ζ (which follows from B = ∇ζ × ∇ψ + Iζ), this implies (9). The constancy of k|| on a flux surface may be used as a test for any numerical scheme.

The constant k|| may also be understood as the magnitude of the poloidal flow where . This result arises because the E × B and diamagnetic perpendicular flows cancel the dpi/dψ and dΦ0/dψ terms in (6) when the poloidal component is formed, leaving only dTi/dψ to drive poloidal flow. The coefficient k|| arises again in the dTi/dψ contribution to the parallel current, as will be shown in section 6:

Here, σneo, , and (defined in section 6) are coefficients determined by the magnetic geometry and electron collisionality, affected by the ions only through the ion charge Z.

The linearized Fokker–Planck–Landau operator for ion–ion collisions, needed to solve (1) or (3), may be written as

where , , is the error function,

is the Lorentz operator, and ξ = v||/v. Also, H and G are the non-Maxwellian corrections to the Rosenbluth potentials, defined by and , with the velocity-space Laplacian . The last three terms in (11) (those following together form the 'field part' of the operator. While historically this part of the operator is often replaced with ad-hoc models, here we retain the exact field terms. The concise form (11) of the field operator is derived in equation (7) of [54], and it is exactly equivalent to the full linearized Fokker–Planck–Landau field operator.

3. Local solver

The basic approach to solving the kinetic equation (1) with the full field operator is to treat H and G as unknown fields along with the distribution function g, and to solve a block linear system for three simultaneous equations: the kinetic equation and the two Poisson equations that define the potentials. Figure 1 illustrates the structure of this linear system. The approach is similar to the innovative method described in [28] but was developed independently. Reference [28] (a radially local code) is focused on the banana regime in which g = g(μ, v) is a function of two phase-space variables, whereas in the analysis here we wish to keep the collisionality general, which means g depends on three or four phase-space variables (in the local and global cases, respectively.)

Figure 1. Refer to the following caption and surrounding text.

Figure 1. Block structure of the linear system for the local solver. A few rows are also reserved for boundary conditions.

Standard image

We may solve either (1) for f1 or (3) for g. The operator and matrix are the same for the two approaches, but the right-hand side vector (the inhomogeneity) is different. The equivalence of the distribution functions obtained by the two approaches is another useful test of convergence. For the second approach, the inhomogeneous term in (3) may be evaluated explicitly:

Deriving this result amounts to evaluating Ci{v||v2fMi}, which is done in e.g. equation (C19) of [55].

We discretize the Rosenbluth potentials by retaining a finite number of Legendre polynomial modes P(ξ). There are several motivations for this choice. First, the Legendre amplitudes of H and G fall off rapidly with ℓ since . Therefore only 2–4 modes are sufficient for convergence, although the code allows for the retention of an arbitrary number of modes. Secondly, the Legendre representation allows a convenient and efficient treatment of the boundary at large v, which can be understood as follows. The distribution function will be within machine precision of zero for v > 6vi, so it is wasteful to store g for this v region. However, H and G scale as powers of v rather than as , so they remain nonnegligible even for v > 6vi. (In fact, for general g, G increases with v.) However, with a Legendre representation , we may exploit the fact that for v > vMax = 4–6vi, the defining equation for H becomes , and so H = Av−(ℓ+1) + Bv. The physical solutions have B = 0, and so the Robin boundary condition v dH/dv + (ℓ + 1)H = 0 may be applied at vMax to ensure H ∝ v−(ℓ+1). In the case of , there are four linearly independent solutions for G. Two are homogeneous solutions to as for H above, and two are particular solutions, which vary as v2 times the homogeneous solutions. Thus where G = Cv−(ℓ+1) + Dv + Ev1−ℓ + Fvℓ+2. The physical solutions have D = F = 0 (see e.g. (45) of [14]), leaving one homogeneous and one particular solution. To accommodate both solutions requires a second order equation as a boundary condition. Writing v2 d2G/dv2 + av dG/dv + b G = 0, and inserting G ∝ v−(ℓ+1), then G ∝ v1−ℓ, yields two equations for (a, b), giving the boundary condition v2 d2G/dv2 + (2ℓ + 1)v dG/dv + (ℓ2 − 1)G = 0.

The other boundary conditions applied are as follows: H = 0 and G = 0 at v = 0 for ℓ > 0, dH/dv = 0 and dG/dv = 0 at v = 0 for ℓ = 0, g = 0 at v = vMax, ∂g/∂ξ = 0 at v = 0, and ∂g/∂v = 0 at (v, ξ) = (0, 0). No boundary conditions are applied to g at ξ = ±1 (i.e. the kinetic equation is applied there with one-sided derivatives.)

While it seems essential to represent the pitch-angle dependence of the potentials using Legendre polynomials, the distribution function itself need not be discretized in the same way, and there are many options available for the other coordinates, so a range of different discretization schemes were investigated. A choice of piecewise Chebyshev spectral colocation and finite-difference methods of various orders were implemented for both the x and ξ grids. The spectral colocation approach is highly accurate for given grid resolution. However, as the matrix is denser in the associated coordinate for these approaches, the solver slows more rapidly compared with finite differencing as the grid resolution increases. Thus, for satisfactory numerical convergence, high-order finite-difference methods are often preferable in practice. For discretization in θ, finite-difference methods of various orders and spectral colocation as well as a sine/cosine modal representation have been implemented. The modal approach is extremely efficient for the simple concentric circular flux-surface model, in which case the matrix is sparse in θ. However, for shaped geometry, the matrix becomes dense in θ for the modal approach, so the colocation approach is both more convenient and similarly accurate. In shaped geometry, despite the accuracy of the spectral approaches, finite-difference differentiation again typically gives satisfactory convergence in less time due to the sparsity of the matrix.

The linear system may be solved using a sparse direct algorithm; a Krylov-space iterative solver may be much faster, but convergence of the algorithm then requires an effective preconditioner. One successful preconditioner is obtained by eliminating the off-diagonal blocks in figure 1 as well as the off-diagonal-in-x terms in the energy scattering operator and boundary conditions. If high-order finite-difference derivatives are used in x, convergence typically also requires that a constant ∼νii be added to the diagonal of the kinetic equation. We find the generalized minimum residual method (GMRES) does not converge consistently, while the stabilized biconjugate gradient and transpose-free quasi-minimal residual methods are more reliable.

Several issues regarding null solutions and symmetry properties of the distribution function are discussed in appendix A.

Figures 2 and 3 show typical results of the local code, plotting the flow and thermal conductivity coefficients k|| and kq as functions of aspect ratio and collisionality. Although the code can use general shaped geometry, for all plots in this paper we use the standard concentric circular flux-surface model B ∝ 1/(1 + epsilon cos θ) and b · ∇θ = constant, to facilitate comparison with previous literature on neoclassical theory. We may then take the definition of the ion collisionality to be ν* = νii/(epsilon3/2vi||θ).

Figure 2. Refer to the following caption and surrounding text.

Figure 2. The flow coefficient k|| defined in (6) for concentric circular flux-surface geometry, as computed by the local code (solid) and the Sauter et al formula (15) (dashed). Squares show the Helander–Sigmar formula (17), an approximate analytical treatment of the ν* → 0 limit, for the same four values of epsilon. The arrow indicates the known analytic ν* → 0, epsilon → 0 limit 1.17.

Standard image
Figure 3. Refer to the following caption and surrounding text.

Figure 3. The thermal conductivity coefficient kq defined in (5) for concentric circular flux-surface geometry. Solid curves are computed by the local code. Dashed lines indicate the Chang–Hinton formula (14) for the same four values of epsilon. Squares show the Taguchi formula (16), an approximate analytical treatment of the ν* → 0 limit, again for the same four epsilon. The arrow indicates the known analytic ν* → 0, epsilon → 0 limit 0.66.

Standard image

Several approximate analytic formulae are also plotted. The Chang–Hinton formula for the heat flux [56]

and the formula of Sauter et al [25]

apply to arbitrary aspect ratio, plasma shaping, and collisionality. (Note α in [25] equals −k|| in our notation.) Taguchi's formula for the heat flux [57]

and a formula for the flow coefficient derived on p 216 of [2]

are applicable at arbitrary aspect ratio and shaping in the limit of small ν*. An expression equivalent to the latter result was also given previously in equation (28) of [7]. (There, and α1 = 1.173. Taking for circular flux surfaces yields the same result as (17).) Here, ft = 1 − fc and

As figure 2 shows, (17) does a reasonable job of predicting the low-collisionality limit of k||. The analytic result k|| = 1.17 obtained using a momentum-conserving pitch-angle scattering model collision operator is indeed the limiting value for ν* → 0 and epsilon → 0 as expected, but epsilon must be <0.01 for this value to be a good approximation. Figure 3 shows Taguchi's formula (16) is extremely accurate. The Chang–Hinton formula (14) is less accurate but it correctly captures the trends with epsilon and ν*.

The local code was also compared with published results from the Fokker–Planck code of [22]; the remarkable agreement of the codes is shown in figure 4.

Figure 4. Refer to the following caption and surrounding text.

Figure 4. The local code agrees to high precision with the Fokker–Planck code of [22] for both the flow (a) and heat flux (b) coefficients. Calculations are shown for epsilon = 0.1. Data reproduced with permission.

Standard image

Another set of transport coefficients arise in the analysis of the electrons. The radial particle and electron heat diffusivity are smaller than the ion heat transport and are always dominated by turbulent transport in practice, so we will not discuss them further. Of greater interest are the electrical conductivity and bootstrap current; these quantities are discussed in section 6.

The distribution function obtained by the local solver has several noteworthy features. In the ν* ≪ 1 limit, analysis shows the g piece of the distribution function should vanish in the trapped region of phase space [1, 2], and a boundary layer exists between the trapped and passing regions [58]. These properties are reproduced in the code, as illustrated in figure 5. The thickness of this boundary layer increases with collisionality. Although the collisionality ν* is typically defined in terms of the thermal speed vi, each value of v in phase space effectively has its own collisionality given by νD/(v epsilon3/2b · ∇θ), with faster particles being less collisional. As shown in the figure, the boundary layer indeed grows narrower with v. In figure 5(c), g at each v is scaled to go to 1 at ξ = 1 for clarity. Also notice in figures 5(a) and (b) that g is nearly constant along particle trajectories, as it should be.

Figure 5. Refer to the following caption and surrounding text.

Figure 5. (a) Particle orbits, i.e. contours of magnetic moment μ, for epsilon = 0.3. (b) The 'collisional response' distribution g at v = vi for ν* = 0.01, showing g ≈ 0 for the trapped region as predicted by banana-regime analytic theory. (c) A slice of g at θ = 0 shows the boundary layer, which is narrower at larger v due to the lower effective collisionality.

Standard image

4. Global kinetic equation

Under what circumstances does the ion distribution remain close to a Maxwellian flux function, i.e. is the approximation fi ≈ fMi(ψ) of section 2 valid? The magnitude of the correction to the flux-function Maxwellian may be estimated from the local theory roughly as f1 ∼ F ∼ (ρθ/r)fMi where r is the scale length of variation in density and/or temperature. In a pedestal, since ρθ/r ∼ 1, then f1 ∼ fMi and the neoclassical expansion breaks down.

However, a more careful estimate reveals there is a regime in which the near-Maxwellian assumption is still valid [59]. To define this regime, first write fMi in an equivalent form:

where

and again W0 = miv2/2 + ZeΦ0 is the leading-order total energy. Then the derivative ∂fMi/∂ψ (at fixed W0) that determines the magnitude of f1 is

(equivalent to (2).) In this form, it is apparent that the magnitude of ∂fMi/∂ψ is determined by rT and rη, the scale lengths of Ti and η, but not directly by rn, the scale length of density. Therefore f1/fMi may be small compared with unity even when rn ∼ ρθ as long as rT and rη are ≫ρθ.

This 'weak-T'i pedestal' regime is the ordering [53, 59] we shall consider for the rest of the analysis: δ ≪ 1 where δ = ρθ/rT, ρθ/rη ∼ δ, and ρθ/rn ∼ 1. This regime is useful for two reasons: the collision operator may be linearized, and, as we will show, the poloidal electric field may be decoupled from the kinetic equation, eliminating the E × B nonlinearity. Therefore the kinetic equation is linear in f1. For rT ∼ ρθ and/or rη ∼ ρθ, a full-f nonlinear kinetic equation must be solved, retaining both the collisional and E × B nonlinearities. Notice rη ≪ 1 implies the ions are electrostatically confined (dΦ0/dψ ≈ −[Zeni]−1dpi/dψ) with (Ze/Ti)dΦ0/dψ ∼ 1/(RBθρθ), so ZeΦ0/Ti ∼ 1. Due to this ordering for the electric field, the vE · ∇fi term in the kinetic equation, neglected in conventional neoclassical calculations, becomes comparable to the v||||fi streaming term. Therefore, although f1 ≪ fMi in the weak-T'i pedestal, conventional neoclassical results must still be modified.

Now, consider the full-f drift-kinetic equation [60]

where S represents any sources/sinks and Cnl is the nonlinear Fokker–Planck–Landau operator. As pointed out by Hazeltine [60], (22) may be derived recursively, and so its validity does not require |v||||fi| ≫ |vd · ∇fi|. Since fi ≈ fMi to leading order, Cnl may be approximated with Ci, the operator linearized about fMi. For the drift velocity vd it will be convenient to use vd = (v||/Ω)∇ × (v||b) (discussed in appendix B) where the gradient acts at fixed μ and W. We make the ansatz Φ1 ∼ δΦ0 and ∂Φ1/∂ψ ∼ δ0/dψ, and we will show shortly that these assumptions are self-consistent. We define f1 by fi = fMi − (ZeΦ1/Ti)fMi + f1, and change the independent variable from W to W0. Neglecting several Φ1 terms that are small in δ,

The contribution from Φ1 to vd · ∇θ is O(δ) smaller than the Φ0 contribution, and (vE · ∇ψ)/(vm · ∇ψ) ∼ ZeΦ1/Ti ∼ δ where vm = vd − vE is the magnetic drift, so Φ1 may be entirely neglected in vd and (23). We therefore approximate vd in (23) with vd0 = vm + vE0 where vE0 = cB−2B × ∇Φ0. To evaluate Φ1 we may use the adiabatic electron density response ne + (eΦ1/Te)ne, where ne(ψ) = Zni(ψ) is the leading-order electron density, with quasineutrality to obtain

Hence, as f1 ∼ δfMi, the ordering ansatz for Φ1 above is self-consistent. As both the collisional and E × B nonlinearities are thus formally negligible, (23) is completely linear.

Just as in the local case, it is convenient to define the collisional response part of the distribution function g using (4). Eliminating f1 in (23) in favor of g, a pair of terms cancels. We also drop the resulting vd0 · ∇[(Iv||/Ω)∂fMi/∂ψ] term because vd0 · ∇(Iv||/B) ∝ ∇ × [(v||/B)B] · ∇(Iv||/B) = 0 exactly and vd0 · ∇(∂fMi/∂ψ) is small in δ. Thus, we obtain

The advantage of this second form is that it makes clear that gradients in ni, Φ0, and/or η cannot affect the g part of the distribution function—only a Ti gradient can drive g. The logic is the same as in the local case: Ci{v||fMi} = 0, so the only gradient surviving in the inhomogeneous drive term (13) is dTi/dψ. This property is obscured in the form (23). While the independence of g from dni/dψ and dΦ0/dψ was well known previously for the local case, it is noteworthy that this property persists in the weak-T'i pedestal case considered here [59].

5. Changes to flow structure

Two noteworthy differences between the local and global analyses are that the parallel flow coefficient k||, as defined in (6), no longer needs to be constant on a flux surface, and it no longer also describes the poloidal flow. To see the first of these points, we may apply the operation ∫d3v to the kinetic equation (25), as detailed in appendix B. The resulting mass conservation equation (ignoring sources) is

Recall from section 2 that in the local case, the v|| term in (26) dominates the others, implying ∫d3v v||g ∝ B. This result was crucial for proving the constancy of k|| on a flux surface, for the dTi/dψ term in the parallel flow is precisely ∫d3v v||g. However, in the global case, (26) indicates that ∫d3v v||g need not vary on a flux surface in proportion to B, so the proof for the constancy of k|| breaks down.

These same results can be derived from a fluid perspective, making no reference to the drift-kinetic equation, starting instead from the fluid mass flow

where . Equivalently,

where ΓE = ∫d3v fivE, is the diamagnetic flow, , , and . The equivalence of (27) and (28) follows from

which may be derived [1] using the more accurate drift . (This drift is identical to our earlier expression to leading order in β ≪ 1.) Notice b· (28) with (4) and V|| = Γ · b/ni gives (6) with k|| = ZeB2〉(cIniB dTi/dψ)−1∫d3v v||v as before.

We now impose mass conservation ∇ · Γ = 0, substituting (4) into (27), applying B × ∇ψ = IB − R2B2ζ, and noting ∫d3vvmfMi = c(ZeB2)−1(dpi/dψ)B × ∇ψ + ∇ × [cpib/(ZeB)]. Cancellations occur to leave T1 + T2 + T3 + T4 = 0 where T1 = ∇ · ∫d3v(v||b + vE0 + vm)g, T2 = −∇ · ∫d3v(vE0 + vm)ZeΦ1fMi/Ti, T3 = ∇ · ∫d3v vE1fMi = ∇ · (cniB−2B × ∇Φ1), and where vE1 = cB−2B × ∇Φ1. As ∇ni = −(Zeni/Ti)∇Φ0 + O(δ) in our ordering, T2 and T3 cancel to leading order in δ. It can be verified that the terms in T4 are O(δ) smaller than the terms in T1, so to leading order, T1 = 0, which is precisely (26), but re-derived from a fluid rather than drift-kinetic perspective. The fluid analysis thereby confirms k|| is no longer constant on each flux surface. Compared with the fluid analysis in the conventional ordering, reviewed following (9), it can be seen that two new contributions to mass conservation become important: E × B convection of the poloidally varying density, and radial variation of the particle flux or (equivalently) diamagnetic flow. Even though |vE0| ≪ vi, E × B convection of the density carried by g matters for mass conservation (26) because the parallel flow only enters multiplied by the small factor Bθ/B. And although , the radial derivative in ∇ · Γ means the next-order correction to in the diamagnetic flow (or equivalently the radial neoclassical flux) must be retained to accurately compute ∇ · Γ.

The poloidal fluid flow Vθ is found by computing Vθ = Γ · eθ/ni, using (27) or (28). Plugging (4) into eθ · (27), several cancelations occur, leaving

So far no terms have been dropped. We now order the terms using the orderings developed in section 4. Using (Ze/Ti)dΦ0/dψ ∼ 1/(RBθρθ) it can be verified that V1 ∼ V2. It can also be verified that each term following V3 is O(δ) smaller than V1 ∼ V2 using ∫d3v g ∼ δni, ZeΦ1/Ti ∼ δ, ∇Φ1 · B × eθ ≈ −IBθ∂Φ1/∂ψ ∼ δIBθ0/dψ, and noting the quantities in square brackets cancel to leading order.

It remains to evaluate V3. The leading-order contribution comes from the radial gradient of the integral of g, since only this derivative has the short scale length ρθ. Thus, we obtain

In the local case, only the v|| term arose in the analogous integral for Vθ. In the pedestal we may define a normalized poloidal flow

so kθ → k|| in the local limit. The property Vθ ∝ dTi/dψ from conventional theory persists in the pedestal, due to (31) and g ∝ dTi/dψ.

6. Electron kinetics and parallel current

The orbit width for electrons is thinner than that of the ions, so direct finite-orbit-width effects for electrons may be neglected. However, the electrons are affected by modifications to the main ion flow. To demonstrate this point, and to show applications of the local Fokker–Planck code to electron quantities, we now analyze the electron kinetics. Since the particle and electron heat transport are essentially always dominated by turbulent transport, we focus here instead on the neoclassical conductivity and bootstrap current. Though the analysis below uses the pedestal ordering, conventional results for the parallel current are exactly recovered in the appropriate limit of the expressions derived here.

Using the gauge of appendix C, the electron kinetic equation may be written

where vme is the electron magnetic drift, w = mev2/2 − eΦ is an independent variable, and Ce is the total electron collision operator. We assume fe ≈ fMe where

Then Ce = Cee + Cei where Cee is equivalent to (11) but with ion quantities replaced by electron quantities, Cei{fe1} ≈ νeiL{fe1} + fMeνeimev||Vi||/Te, , x = v/ve, and . We write fe = fMe exp(eΦ1/Te) + mev||V||fMe/Te + h and solve for h. We also make a change of independent variables in the kinetic equation to w0 = mev2/2 − eΦ0. Using (6), the leading terms in δ and are

where Ωe = −eB/(mec), h0 is the first term in a series h = h0 + h1 + ···, and the inductive term has been taken as higher order. The solution to (35) may be written where p = pe + pi, and hp and are the solutions to

with De = v|||| − Cee − νeiL. Recalling eΦ1/Ti ∼ δ, the O(δ) terms in the kinetic equation are

where vE1 = cB−2B × ∇Φ1, and we have assumed . The solution may be written (dTi/dψ)hΦ/e where ρ0 = vimic/(ZeBav), and . Here, hE, , and hΦ are the solutions of

Note in the local case, ∇||k|| = 0 so the last term in (40) vanishes and . The De operator, which is radially local in that ψ is merely a parameter, may be inverted numerically for the right-hand sides (36)–(37) and (39)–(41) just as described in section 3 for the similar ion operator v|||| − Cii. Then the parallel current is

Now consider the result of applying to (39). This operation annihilates both the right-hand side and the collision operators in De, leaving (∂/∂θ)∫d3v v||hE/B = 0. Therefore the flow carried by hE is ∫d3v v||hE = αEB for some flux function αE. The same logic applies to (37), so for some flux function . Applying ∫d3v to (36), the right-hand side is not annihilated this time, and we instead find ∫d3v v||hp = B−1 + αpB for some flux function αp. Lastly applying ∫d3v to the first equation in (40) and to (41), we obtain and ∫d3v v||hΦ = αΦB − ng/(ZB) where and αΦ are flux functions, ng = Ti(ρ0Ini dTi/dψ)−1∫d3v g is the O(1) normalized density carried by g, and we have invoked (24). Thus, the dTi/dψ term in the parallel current varies poloidally ∝B in the local case where k|| is constant, but not in the global case where k|| varies.

Putting the pieces together, the total parallel current is

where α is another flux function. Multiplying this equation by B, flux-surface averaging, and substituting the result back into (43), we obtain

The dp/dψ term is the standard Pfirsch–Schlüter current, and the 〈j||B〉 term is the ohmic and bootstrap contribution. However, the k|| and ng terms are new in the global case, vanishing in the local case where k|| is constant and |∇ne|ρθ ≪ 1. We may write the ohmic and bootstrap contribution as

where σneo = 〈B∫d3v v||hE〉, , , , and . The term in (45) is new, becoming negligible in the conventional case. For the local case of constant k||, where , it is useful to define so is completely independent of all ion quantities except Z. The definitions of σneo, , and here are consistent with [25]. Interestingly, the new ng terms in (44) and (43) and the new term in (45) are quadratic in the gradients.

Figure 6 shows these coefficients of the bootstrap current and the conductivity as calculated by our code for the local limit k|| = constant, using the circular flux-surface model and Z = 1. The conductivity has been normalized by the parallel Spitzer value. The analytic fits to numerical calculations of the coefficients by Sauter et al are plotted for comparison [25]. The horizontal coordinate in these plots is , which is smaller than the ν*e defined in [25]. We find the Sauter expressions give an excellent fit to the coefficients in the banana regime, though there is some discrepancy at higher collisionality when epsilon > 0.1, the same pattern observed in figure 2. The reason for the discrepancy is unclear, since the fundamental kinetic equations and collision operators we use to generate figures 2 and 6 are identical to those solved by CQLP, the code to which the Sauter expressions are fit. (CQLP uses an adjoint method whereas our results do not, though this difference should not affect the physical results.) We have verified the difference persists when D-shaped Miller equilibrium is used, and the code of [22] produces identical coefficients to ours. As shown in figure 7, the difference between our coefficients and those of [25] can lead us to predict a reduced total bootstrap current density in the pedestal for experimentally relevant plasma parameters when ν*e > 1. This difference is primarily due to our lower , which multiplies the large dp/dψ term. When ν*e < 1, our prediction for the total bootstrap current density becomes indistinguishable from that of [25].

Figure 6. Refer to the following caption and surrounding text.

Figure 6. Parallel current coefficients defined in (45) computed by our local code (dots, connected by solid curves). Dashed curves show the semi-analytic formulae of Sauter et al [25].

Standard image
Figure 7. Refer to the following caption and surrounding text.

Figure 7. (a) Model profiles resembling the DIII-D measurements in figure (8) of [44] and a plausible q profile yield the collisionality profiles in (b). (c) The local code described here predicts identical total bootstrap current density to the formulae of [25] for lower collisionality. (d) We predict somewhat lower current density at the real collisionality, due primarily to the discrepancy in shown in figure 6(b). (e) At higher collisionalities, the differences become significant.

Standard image

As with the flow, the total current vector j remains divergence free in the pedestal ordering:

where is the total diagonal anisotropic stress, including O(p) and O(δp) terms. Equation (46) can be proved from (44) and (26) using (24), fe ≈ fMe + eΦ1fMe/Te, and (29).) Equation (46) indicates the new k|| and ng terms in (44) arise for the same fundamental reason as the conventional Pfirsch–Schlüter current: a parallel current must flow to maintain ∇ · j = 0 given the perpendicular diamagnetic current. In the pedestal, the diamagnetic current associated with the poloidally varying pressure becomes large enough to modify the parallel current on the level of the dTi/dψ terms.

7. Global numerical scheme

It is equally valid and equally numerically challenging to solve either (23) or (25). For the rest of the analysis here we discuss the case of (25) for definiteness.

As we are interested in a narrow radial domain around the pedestal, we assume I, B, and ∇||θ are independent of ψ for simplicity. These approximations are also convenient as they make vm · ∇θ = 0 exactly for our form of the drifts. For simplicity, we also take η and Ti to be constant over the simulation domain. The one place where dTi/dψ must be retained is in the inhomogeneous term, since the drive is ∝dTi/dψ. As the kinetic equation is linear, g may be normalized by dTi/dψ, while every other appearance of Ti is treated as a constant.

Both versions of the global kinetic equation (23) and (25) resemble their local counterparts (1) and (3), but with the additional vd0 · ∇ term in the unknown. Due to the radial derivative in this term, the radial coordinate no longer enters the kinetic equation as a mere parameter, meaning the problem is now four-dimensional: g = g(ψ, θ, μ, W0). In these original variables, the allowed range of each coordinate depends on the other coordinates in a complicated manner. For numerical work it is therefore convenient to change the independent variables from (μ, W0) to (v, ξ) so the coordinate ranges become coordinate-independent. In these variables, the kinetic terms in (23) and (25) become

where

is the drift-kinetic operator implemented in conventional neoclassical codes [16, 22], and

consists of new terms proportional to the radial electric field. In a pedestal, not only is the ∂g/∂ψ term in (47) important, but the terms in KE also become equally important. The aforementioned ordering for dΦ0/dψ implies each term in KE comparable in magnitude to K0. Physically, the latter two terms in (49) are essential for maintaining conservation of μ and total energy as a particle's kinetic energy changes during an orbit. This kinetic energy changes because the electrostatic potential seen by the particle varies over an orbit width.

We choose ni(ψ), which determines Φ0 = (Ze)−1Ti ln(η/ni). On either end of the radial domain, we take ni(ψ) and Φ0(ψ) to be uniform for a distance of several ρθ, as illustrated in figures 8(a) and (b). In this way, the distribution function will approximate the local neoclassical solution at the radial boundaries, so local solutions can be used there as inhomogeneous Dirichlet boundary conditions.

Figure 8. Refer to the following caption and surrounding text.

Figure 8. (a) Equilibrium density profile for the global calculation, normalized to its value at the left boundary. (b) Normalized radial electric field −cI(viB0)−10/dψ. (c) Profile of ν*. (d) Parallel flow coefficient k||, evaluated at the outboard midplane (θ = 0). Black dashed curve indicates the result of the local neoclassical code. The other curves (nearly indistinguishable) demonstrate the convergence of the global code to the various numerical resolution parameters. (e) Radial heat flux, with the same legend as (d). (f) Normalized poloidal flow kθ and k|| evaluated at the outboard and inboard (θ = π) midplanes.

Standard image

As the kinetic equation (25) is linear, it may in principle be solved numerically using a single matrix inversion. Indeed, this is the approach traditionally adopted by local neoclassical codes [1618, 22], including the one described in section 3. However, this approach is already somewhat numerically challenging for the local problem due to the three-dimensional phase space, as the matrix has dimension (Nθ Nv Nξ× (Nθ Nv Nξ), where Nθ, Nv and Nξ are the number of modes or grid points in the respective coordinates. In the nonlocal case, the additional spatial dimension means the matrix size must increase to (Nψ Nθ Nv Nξ× (Nψ Nθ Nv Nξ) for Nψ radial grid points, making the solution much more time- and memory-intensive. Therefore we seek an alternative method.

In the new approach proposed here, a derivative with respect to a fictitious time ∂g/∂t is first added to the left-hand side of (25). For reasonable initial conditions and boundary conditions, g should evolve towards an equilibrium since the equation (25) is dissipative. However, an explicit time advance requires very small time steps for stability due to the many derivatives in the kinetic equation, and an implicit time advance would require the inversion of a matrix just as large as for a direct solution of the original time-independent equation.

An effective solution is to employ the following operator-splitting technique. Consider the following series of two backwards-Euler time steps:

where KL = K0 + KE − Ci is the 'local operator' and KNL = (vm · ∇ψ)∂/∂ψ is the 'nonlocal operator'. When summed together, gt+(1/2) cancels, leaving an equation that is equivalent to first order in Δt to a backwards-Euler time step with the complete operator KNL + KL. However, each of the steps (50)–(51) are much easier than a step with the total operator because the dimensionality is reduced: e.g. ψ is only a parameter in (50), so this step requires the inversion of Nψ matrices, each of size (Nθ Nv Nξ× (Nθ Nv Nξ). Also note that the local operator at each radial grid point need only be LU-factorized once, with the L and U factors reused at each time step for rapid implicit solves. The same is true of the nonlocal operator at each v and ξ.

Several higher order operator-splitting schemes were explored, but none were found to be stable for the equation here.

The procedure outlined here provides a general recipe for extending a conventional neoclassical code into a pedestal code. A conventional neoclassical code inverts an operator K0 − Ci, i.e. many of the terms in KL, so minor modifications would allow such a code to carry out the local part of the time advance. The modifications necessary are adding the electric field terms KE and adding the diagonal associated with the time derivative. The resulting operator is then iterated with the nonlocal operator.

For the results shown here we employ a piecewise-Chebyshev grid in ψ with spectral colocation differentiation. A tiny artificial viscosity is required at the endpoints for numerical stability; the magnitude of this viscosity may be varied by many orders of magnitude with no perceptible change to the results. Inhomogeneous Dirichlet radial boundary conditions are imposed, with the distribution function at these points taken from the local code. For completeness, we have also tried upwinded high-order finite differences for radial differentiation, with the upwinding direction opposite above and below the midplane, corresponding to whether drift trajectories in the region move towards increasing or decreasing ψ. For our sign convention, the magnetic drifts are downward, so the inhomogeneous Dirichlet radial boundary condition must be specified above the midplane at large minor radius and below the midplane at small minor radius. This radial discretization scheme gives equivalent results to the Chebyshev method, but it requires more grid points for convergence, and a numerical instability tends to arise at large times.

8. Need for a sink

In order to reach equilibrium, it is essential to include a heat sink. This requirement may be understood physically as follows. As we take the scale lengths at each radial boundary to be large compared with ρθ, the heat flux into the volume at small minor radius and the heat flux out of the volume at large minor radius are determined by the local neoclassical result (5). These fluxes are different due to the different densities at the two boundaries, and so net heat will constantly leave (or enter) the simulation domain. More rigorously, as shown in appendix B, the 〈∫d3v( · )〉 and moments of the kinetic equation in steady state give

The first equation represents local mass conservation, and the quantity following V' is the particle flux. The particle flux is exactly zero in the local limit, so it vanishes at the radial boundaries, and so in the absence of a source/sink, it must vanish everywhere in the domain. In the second equation, representing global energy conservation, the first term is the difference between the heat into and out of the domain, and the second term represents change in electrostatic energy associated with particle flux. If S = 0, then the latter two of the three terms in (53) vanish, but the first term is nonzero because the heat fluxes at the two radial boundaries are generally unequal. This contradiction proves the kinetic equation has no steady-state solution without a sink S.

In a real pedestal, there will be a divergence of the turbulent fluxes, which would act as a sink term in the long-wavelength (drift-kinetic) equation we simulate here. Determining the phase-space structure of this turbulent sink term from first principles is an extremely challenging task, beyond the scope of this work. We therefore use a variety of ad hoc sink terms, and we find the simulation results are only mildly sensitive to the particular choice of sink.

The standard sink we use is

where γ is a constant. The sum over signs of ξ ensures that S vanishes exactly for an up-down symmetric magnetic field in the local limit, due to the parity of the local solution discussed in appendix A. This sink is quite similar to the one described in [61] for global δf gyrokinetic codes. The constant γ may be varied by several orders of magnitude without major qualitative changes to the results.

Another option we consider for the sink is

where γm and γp are constants, n1 = ∫d3v g, and p1 = ∫d3v (miv2/3)g. The first term in (55) dissipates any mass in g, while the second term dissipates any energy in g.

9. Results

Figures 810 show results of the global calculation for a pedestal with epsilon = 0.3. The simulation domain consists of an annular region in space, i.e. an interval in ψ. The density varies by roughly a factor of 3 from the top of the pedestal to the bottom, with the profile of dimensionless n = ni/ni(r = −) shown in figure 8(a). This density profile implies the electric field profile shown in figure 8(b), which reaches a minimum of ≈ −0.5viBθ/c in the pedestal center. The collisionality ν* ranges from 0.5 to 0.15 over the domain. (We choose this arbitrary range close to one just to emphasize that ν* is not formally large or small in this formulation.) In these plots, the radial coordinate r/ρθ is defined by r/ρθ = ZeB0(micviI)−1ψ where B0 is the toroidal field on axis. The radial location r = 0 is an arbitrary minor radius, (here the middle of the pedestal), not the magnetic axis. The sink used is (54) with γ = 0.1ωt, where ωt = vi||θ is the transit frequency. The simulation is nearly converged by t = 30/ωt, but very small changes in the results continue until t = 200/ωt. We plot results for t = 200/ωt since doubling this duration produces no visible change to the results. By t = 200/ωt, the residual, which we define as a sum over all phase-space grid points of |∂g/∂t|, has been reduced to 0.05% of its initial value.

It is also important to verify that the code has converged with respect to the many other numerical parameters. Figure 8(d) shows the parallel flow coefficient k|| at the outboard midplane for 11 global runs, all with the same physics parameters, but varying each numerical parameter by a factor of two: simulation duration (tmax), time step (dt), artificial radial viscosity, number of poloidal modes (Nθ), number of Legendre polynomials in the Rosenbluth potentials (NL), number of grid points in v, ξ and ψ (Nv, Nξ and Nψ), and domain size in speed (vmax) and radius (ψmax). The changes are barely perceptible, demonstrating very good convergence. For comparison, the profile computed by the local code is also plotted, calculated by solving (3) (i.e. a single linear system solve) at each radial grid point. The local coefficient varies across the pedestal due to the change in collisionality. Resolution parameters were, unless doubled, Nθ = 6, Nψ = 29, Nξ = 25, Nv = 16, rmax = 4ρθ, vmax = 5vi, NL = 2 and dt = 0.01ωt. Running in Matlab on a single Dell Precision laptop with Intel Core i7-2860 2.50 GHz CPU and 16 GB memory, the base case global simulation took roughly 3 hours to reach t = 200ωt, though runs could undoubtedly be greatly expedited if the code were parallelized and rewritten in fortran. Work to this end is underway. The local solver for these parameters took 0.5 seconds per radial grid point.

Figures 8(d)–(f) show the heat flux kqn2 and the flow coefficients k|| and kθ at the outboard and inboard midplanes. (It is radial variation in the heat flux kqn2 and not kq itself that determines local heating, as shown by (5) and (53)). Outside of the pedestal, as expected, these coefficients agree with the local prediction, and k|| and kθ are equal and poloidally invariant. In the pedestal, however, all coefficients are substantially altered from the local prediction; k|| and kθ differ and vary poloidally. The radial heat flux profile is flattened relative to the local prediction.

For the parameters used here, k|| and kθ change sign in the simulation within the pedestal. These coefficients may have either sign in conventional theory, as shown in figure 2, depending on collisionality and magnetic geometry. In both the local and global cases, the integrand v||g that determines k|| ∝ ∫d3v v||g is positive in part of phase space and negative elsewhere, and the balance between these regions determines the overall sign.

Structure with a radial scale comparable to ρθ is observed in the flow coefficients. Ions 'communicate' over distances comparable to the orbit width , and so the effects of the driving electric field well are felt outside of the well itself, with influence decaying on the orbit width scale. The behavior of the flow coefficients on either side of the well need not be monotonic, as the mean flow adjacent to the well arises from a complicated interplay of particles entering from regions of differing collisionality, some particles directly affected by the electric field and some not.

The flow coefficients are also observed to be nonmonotonic functions of r. This behavior is not unreasonable given the radial localization of the electric field well: even the local distribution function is a nonmonotonic function of r, since it is a complicated function of the radially varying collisionality, so in the global case the flow coefficients need not be monotonic in r.

Figure 9 shows the poloidal variation of the flow coefficients near the pedestal top and bottom. As discussed in the appendix, the drift terms in the kinetic equation break the symmetry which the distribution function possesses in the local case, and so the flow coefficients need not be even or odd in θ.

Figure 9. Refer to the following caption and surrounding text.

Figure 9. Poloidal variation of the parallel flow coefficient k|| and normalized poloidal flow kθ at two radial locations straddling the pedestal.

Standard image

The various components of the mass conservation equation (26) were each independently computed from g: ∇ · ∫d3v v||gb, ∇ · ∫d3v vE0g, ∇ · ∫d3v vmg and ∫d3v S. The first three of these integrals summed to nearly zero everywhere in space, with the sink integral negligible in magnitude compared with the others. Thus, the sink has little effect on the mass conservation relation that effectively determines the flows. The vm integral was intermediate in magnitude, leaving a dominant balance between the ∇ · ∫d3v v||gb and ∇ · ∫d3v vE0g terms. The density ∫d3v g has a ∝ cos(θ) behavior in the pedestal, resulting in ∇ · ∫d3v vE0g ∝ (∂/∂θ)∫d3v g ∝ − sin(θ). To balance this term in the mass conservation law, k|| must develop a − cos(θ) structure, which can be seen in figure 9. Outside of the pedestal, the vE0 term becomes negligible due to the reduced |dΦ0/dψ|, so this drive for poloidal variation in k|| is absent.

Figure 10 shows how the results are altered when different choices are made for the sink. For the sink (54), we show results for γ = 0.1ωt (the value used for all other plots) and for γ = ωt. We also show results for the alternative sink (55). For comparison, results are also shown for a run in which no sink was included. For this run the code did not converge in time, due to the constant loss of heat described in section 8, so it was stopped at t = 30ωt (a time before the heat loss becomes excessive, but after the runs with sinks have nearly converged.) The various options yield results that show the same qualitative modification of the coefficients: a well develops in k|| and kθ at the outboard midplane, and the heat flux profile is flattened relative to the local prediction.

Figure 10. Refer to the following caption and surrounding text.

Figure 10. Using the same equilibrium density profile (a) and Er profile (b) as in figure 8, the parallel flow coefficients (c)–(d) and heat flux (e) show some minor dependence on the the choice of sink. However, qualitative features such as the well in k|| and kθ at the outboard midplane are robust.

Standard image

As a further test of the code, we repeated the numerical calculation, treating f1 as the unknown quantity instead of g, in which case the inhomogeneous term in the equation is −vD · ∇fMi instead of Ci(F). Despite the very different phase-space structure of these two source terms, the numerical results from the two approaches agreed, as they should.

10. Discussion

In this work we have demonstrated a method to extend neoclassical calculations to incorporate finite-orbit-width effects in a transport barrier for the case of Bθ ≪ B (which is a good approximation in standard tokamaks) and a relatively weak ion temperature gradient. The method is implemented in a new continuum δf code. Operator splitting is used to improve numerical efficiency, and we have demonstrated that excellent convergence is feasible for experimentally relevant parameters. By construction, the method exactly reproduces conventional (local) results in the appropriate limit of weak radial gradients. The Rosenbluth potentials are solved for along with the distribution function at each step, allowing use of the full linearized Fokker–Planck–Landau collision operator.

A principal finding of this work is that the parallel and poloidal flows may differ significantly from the conventional predictions. While the coefficients of the poloidal flow and dTi/dψ-driven parallel flow are equal in conventional theory, in the pedestal these two coefficients (k|| and kθ) differ. And, while the poloidal variation of the poloidal flow is ∝Bθ in the core, the same is not necessarily true in the pedestal. The poloidal variation of the flow is effectively determined by mass conservation, and in the pedestal, two new terms become important which are normally neglected: E × B convection of the poloidally varying density, and radial variation of the particle flux (which can be related to diamagnetic flow from the correction to the pressure). These effects cause the parallel flow coefficient k|| to take a well-shaped radial profile, with different magnitude and (for some parameters) opposite sign, relative to the conventional local neoclassical result. In addition, the flow coefficients exhibit a strong poloidal variation not previously found. While this poloidal variation resembles cos(θ), it contains other harmonics and is asymmetric about the midplane due to the magnetic drift.

These issues may be important for comparisons of experimental pedestal flows to theory [11, 12]. In general, the flow coefficients may differ in both magnitude and sign relative to local theory, as shown in figure 8. The fluid flow exhibits strong shear, with radial variation on the ρθ scale.

Associated with the modification to the flow, the parallel current is also modified. Due to the additional terms which must be included in the mass conservation equation, the usual division of the parallel current into Pfirsch–Schlüter and ohmic-bootstrap components is modified, as shown in equation (44). In addition, the dTi/dψ contribution to the bootstrap current is altered, as shown in (45). In the δf formulation, the ion temperature scale length cannot be as small as the density scale length in the pedestal, so these modifications to the parallel current are modest. However, similar changes to the current would presumably occur in a full-f calculation when rT ∼ ρθ, giving order-unity changes to the Pfirsch–Schlüter and bootstrap currents in that case. This issue needs to be examined in future studies.

In the development of this work, the local code was also used to test several analytic expressions for conventional neoclassical theory. The flow and heat flux coefficients k|| = kθ = 1.17, kq = 0.66 derived using the momentum-conserving pitch-angle scattering model for collisions are a poor approximation unless epsilon is ≪0.1. Expressions (16)–(17) are a much better approximation at realistic aspect ratio. The Chang–Hinton heat flux captures the trends at finite epsilon and ν* well and gives results correct to within 20%, at least for the circular concentric flux-surface model. The semi-analytic local formulae of Sauter et al [25] for the flow and bootstrap current coefficients were found to be in excellent agreement with our local code when ν* < 0.3, but some disagreement was found for 0.3 < ν* < 100. For pedestal profiles typical of DIII-D, the Sauter bootstrap current formula closely agreed with our code at low collisionality, ν*e < 1, but the Sauter formula can give a bootstrap current more than twice ours when ν*e ⩾ 10. The Sauter formulae are intended to reproduce results from a code based on the same physical model as our conventional local code.

There are many ways in which the global calculations can be extended. First, it would be useful to include impurities, for it is typically the impurity flow that is measured rather than that of the main ions, and the flows of different ion species may be significantly different [7]. Also, the presence of impurities can introduce a direct density gradient dependence [7] to g. Second, the method should be extended to allow strong temperature gradients (rT ∼ ρθ). Doing so will require the full bilinear collision operator and a full-f treatment. However, as the weak-T'i case is less complicated to analyze, due to the linearity of the kinetic equation, thorough understanding of this limit using the present approach is important for benchmarking future more sophisticated full-f codes. Finally, studies of the velocity-space structure responsible for fluxes in turbulence codes may yield more accurate forms of the sink term needed in our approach.

Acknowledgments

The authors wish to thank Daniel Told for suggestions regarding the sink term, Michael Barnes for helpful discussions on operator splitting, and Felix Parra and Peter Catto for many enlightening conversations. We are also grateful to S Kai Wong and Vincent Chan for contributing data for figure 4. This work was supported by the Fusion Energy Postdoctoral Research Program administered by the Oak Ridge Institute for Science and Education.

Appendix A.: Null space and symmetry of the distribution

The local drift-kinetic equations (1) and (3) have two null solutions fMi and v2fMi, meaning that the discretized matrix for the local code should be nearly singular. In practice, the matrix is still sufficiently well conditioned that the linear system may be solved without a problem, yielding a distribution function that contains a small amount of the two null solutions. For many applications this may not be a concern, because these null solutions do not contribute to the heat flux and flow.

For an up–down symmetric tokamak (i.e. if B and ∇||θ are both unchanged under θ → −θ), then a symmetry exists in the local kinetic equations: if f1(θ, ξ, v) is a solution, then so is −f1(−θ, −ξ, v) (and similarly for g). This property can be exploited to simultaneously eliminate the null space from the matrix and to reduce its size [22]. This is carried out by forcing f1 to have the above symmetry by representing it as a sum of two types of modes: those that are even in θ and odd in ξ, and those that are odd in θ and even in ξ. The two null solutions do not possess this symmetry, so they are automatically excluded. Furthermore, the matrix size is reduced without loss of resolution. For example, the ξ grid can be reduced to only cover the interval [0, 1] instead of [−1, 1] if all sin(Mθ) and cos(Mθ) modes are retained. The odd-θ (i.e. sin(Mθ)) modes are forced to be even in ξ by application of the boundary condition ∂f1/∂ξ = 0 at ξ = 0, and the even-θ (i.e. cos(Mθ)) modes are forced to be odd in ξ by application of the boundary condition f1 = 0 at ξ = 0.

Even if the parity of the solution is not enforced automatically by the discretization in this manner, the null solutions can still be excluded by enforcing parity as follows. Given a numerical solution f1(θ, ξ, v) that contains some of the null solutions, the combination [ f1(θ, ξ, v) − f1(−θ, −ξ, v)]/2 can be formed; the result will also satisfy the kinetic equation but have the desired parity.

In the global case, the symmetry of the kinetic equation is broken by the drift terms. However, the local operator has no null space in the initial-value-problem formulation due to the extra contribution on the matrix diagonal from the time derivative.

Appendix B.: Conservation laws for the global drift-kinetic equation

Here we sketch the derivation of general conservation equations, from which (26), (52) and (53) can be obtained. For most of this appendix we do not assume axisymmetry, we retain radial variation of magnetic quantities, and we do not require B · ∇Φ = 0. We do require the electric field to be electrostatic and we assume ∂Φ/∂t and ∂B/∂t can be neglected. The derivation applies both to the full-f and δf contexts, since the necessary integrals of both the bilinear and linearized ion–ion collision operators vanish.

We begin with the ion drift-kinetic equation

where vd = (v||/Ω)∇ × (v||b) and C is either the bilinear or linearized Fokker–Planck–Landau operator. Gradients are all performed at fixed μ and total energy W = miv2/2 + ZeΦ (including the total potential Φ, not just Φ0), so vd includes both the magnetic drift vm and E × B drift vE. This form of vd is convenient because it makes the kinetic equation conservative without cumbersome higher order terms. This vd includes an incorrect O(β) ≪ 1 parallel magnetic drift, but this component of the magnetic drift is typically unimportant compared with parallel streaming motion, and in fact vm · ∇θ is precisely zero in the model magnetic geometry we use in the code. It is convenient to first rewrite

Then ∫d3v is applied to (B1), annihilating C. Notice

where σ = sgn(v||). The divergence in (B2) may be pulled in front of the integrals in (B3), as the contributions from differentiating the integration limits all vanish either due to the σ sum or because v|| = 0 at the lower limit of W. Application of several vector identities to the B × ∇fi term then yields a mass conservation equation:

Flux-surface averaging and neglect of Φ1 then gives (52).

An energy conservation equation may be obtained by observing that the above derivation of the mass conservation equation is essentially unchanged if ∫d3v W is applied to (B1) in place of ∫d3v. Subtracting ZeΦ× (B4) from the result, one obtains

In the special case of axisymmetry and ∇||Φ = 0, flux-surface averaging and integration in ψ then gives (53).

To obtain the momentum conservation equation, it is convenient to specialize to axisymmetry at the start, taking the ∫d3v(Iv||/B) moment of (B1), and using vd · ∇fi = (v||/B)∇ · [fi(mic/Ze)∇ × (v||b)] instead of (B2). The divergence may be brought in front of the W and μ integrals as before. Noting vd · ∇(Iv||/B) = 0 and v||||(Iv||/Ω) = vd · ∇ψ, the result may be written

This result holds in axisymmetry even if ∇||Φ and/or ∇||I are nonzero.

Appendix C.: Convenient gauge

Here we prove that the gauge may always be chosen so

Axisymmetry is not required, and the loop voltage need not be uniform. The utility of (C1) is that the inductive part of E|| has simple spatial variation ∝B.

Suppose we begin in a different gauge, denoted by tildes, in which

We may transform to a new gauge using and for a generator χ. We choose

where the integrand is evaluated at t' and θ' rather than t and θ. We must verify (C3) is single-valued in θ so Φ is single-valued. To this end, notice 〈B · ( )〉 applied to (C2) gives . Therefore χ(θ = 2π) = 0 = χ(0), so χ is indeed periodic. Applying ∇|| to with B· (C2) and (C3) then gives (C1) as desired.

undefined