The influence of sonic toroidal rotation on gyrokinetic stability and transport is studied, with important implications for heavy impurity dynamics. When centrifugal drifts and electrostatic trapping corrections are included, significant modifications to the calculated transport of heavy impurities are observed. These high-rotation corrections add to the standard Coriolis drift and toroidal rotation shear drive which are normally included in gyrokinetics. Yet, because of their complexity, centrifugal and electrostatic trapping terms (quadratic in the main ion Mach number) are not generally included in gyrokinetic codes. In this work, we explore the implications of using reduced descriptions of the rotational physics. For heavy impurities such as tungsten, cross terms due to the centrifugal force can dominate the rotation dynamics, and neglecting them is shown to lead to large errors in the impurity particle flux.

External torque driven by neutral beam injection in tokamak plasmas can drive plasma flows comparable to the ion sound speed. The presence of this toroidal plasma flow can have a profound effect on the intensity of drift-wave turbulence and corresponding levels of radial transport. In particular, the radial variation (or shear) in the E × B velocity suppresses turbulence by acting to twist plasma eddies.1,2 In contrast, toroidal flow shear provides a Kelvin-Helmholtz-type linear drive that can increase turbulent fluctuation levels.3 It is common in gyrokinetics to consider the weak rotation limit,4,5 retaining only the E × B flow, Coriolis drift, and toroidal rotation shear. In the gyrokinetic equations, these effects are represented by terms linear in the Mach number, and they provide the dominant symmetry-breaking mechanism for momentum transport. However, correct treatment of the sonic rotation regime requires the additional consideration of centrifugal effects.6–9 Sonic toroidal flow, on the order of the thermal speed, produces a centrifugal force that pushes ions toroidally outward, causing them to redistribute non-uniformly around a flux surface. As a result of quasi-neutrality, a poloidally varying electrostatic potential is generated to balance the density asymmetry, pulling the electrons with the ions.10 The centrifugal force induces a centrifugal drift, while the poloidally varying electrostatic potential produces modifications to the E × B drift and to the particle trapping. Because of their complexity, these new sonic terms (quadratic in the main ion Mach number) are generally ignored in gyrokinetic codes. In this work, we will show that they are nevertheless essential for the study of transport of heavy impurities, such as tungsten from the divertor wall, as the influence of the centrifugal force is amplified by their greater mass.

The first gyrokinetic code to correctly implement sonic rotation was GKW.11 With GKW, it has been shown that the centrifugal force leads to an enhanced mirror force, which increases the trapped electron drive and reduces the residual zonal flows, yielding a small increase in the ion heat diffusivity for ion-temperature-gradient (ITG)-dominated turbulence.12 However, centrifugal effects can be large for heavy impurities, such as tungsten, even when they are negligible for the bulk plasma.13 Studies of nonlinear turbulent transport are more limited and have focused on experimental scenarios.14 

In the present work, we will describe the CGYRO15 implementation of sonic rotation, including all quadratic Mach number terms. A major focus of the present work will be to understand the implications of neglecting centrifugal effects by comparing results from the complete sonic formulation to the simpler weak rotation approximation. We remark that the CGYRO implementation complements the rotation formulation in the drift kinetic code NEO.16–18 NEO implements sonic toroidal rotation using a physically consistent and mathematically equivalent ordering hierarchy.8,19 That is, the equations solved in CGYRO complement those solved in NEO insofar as together they represent the complete first-order deviation of the plasma from a local Maxwellian.

The rest of this paper is organized as follows. In Sec. II, the basic simulation model for full sonic rotation is described. Simulation results exploring the impact of rotation and rotation shear on gyrokinetic linear stability are shown in Sec. III. Benchmarks with GKW are also presented for verification purposes. In Sec. IV, the impact of rotation on heavy impurity transport, such as tungsten and carbon, is explored with both quasilinear and nonlinear simulations. Finally, a brief summary is given in Sec. V.

We adopt the non-orthogonal field-aligned coordinate system (ψ,θ,α), together with the Clebsch representation for the magnetic field20 

B=α×ψsuchthatB·α=B·ψ=0.
(1)

The angle α is written in terms of the toroidal angle φ as αφ+ν(ψ,θ), where the function ν(ψ,θ) is defined with respect to the current function I(ψ) by Eq. (9) in Ref. 21. In Eq. (1), ψ is the poloidal flux divided by 2π, and θ simultaneously refers to an angle in the poloidal plane (at fixed φ) and a parameterization of distance along a field line (at fixed α). In these coordinates, the Jacobian determinant and the flux-surface average are, respectively,

Jψ1ψ×θ·αandfdθdφJψfdθdφJψ.
(2)

We consider solution of the gyrokinetic equation using the sonic rotation formalism of Sugama,8 which allows for flow velocities on the order of the ion thermal speed. The formalism is analogous to the gyrokinetic equations in the comoving frame of a toroidally rotating plasma developed by Peeters et al.9 on which the GKW code is based. The distribution function and fields are expanded in powers of ρ*=ρi/a, the ratio of the ion gyroradius to the system size. When the rotation is sonic, the electrostatic potential must be one order larger in ρ* than for the diamagnetic ordering usually assumed in kinetic theory, i.e.,

Φ=Φ1+Φ0+Φ1+.
(3)

It has previously been shown by Hinton and Wong19 that in an axisymmetric system (considered here), Φ1 is a flux function, such that the zeroth-order flow velocity

U0=ω0R2φ
(4)

is purely toroidal and species-independent. Here, ω0(ψ)cd(Φ1)/dψ is the angular rotation frequency.

The O(1) equilibrium equation guarantees that the zeroth-order distribution function retains the Maxwellian form, but now with respect to the rotating velocity frame

f0a=na(ψ,θ)(2πvta2)3/2eua2.
(5)

Above, ua=v/(2vta) is the normalized speed in the rotating velocity frame with respect to the thermal speed vta=Ta/ma, such that v2=v2+v2. While the equilibrium temperature Ta=Ta(ψ) in Eq. (5) is a flux function, the equilibrium density na is not, varying as na(ψ,θ)=Na(ψ)eλa, where λa is the dimensionless effective potential energy, defined as

λa(ψ,θ)zaeTaΦ̃0Ma2R22R02.
(6)

Here, Φ̃0=Φ0Φ0 is the poloidally varying part of the equilibrium potential, and we define the species Mach number as

Ma=ω0R0vta,
(7)

where R0 is the midplane major radius of the flux surface. While Φ0 [like Ta(ψ) and Na(ψ)] is undetermined, Φ̃0 is determined by the O(1) Poisson equation

azaena(ψ,θ)=0.
(8)

This equation is equivalent to the usual quasi-neutrality relation. In general, for an arbitrary number of species, solution of Eq. (8) for Φ̃0 requires solving a nonlinear algebraic equation. In CGYRO, we solve this numerically using Newton's method (we remark that the same equation must also be solved in NEO). For the simple case of a pure plasma, Eq. (8) can be solved analytically, yielding

(1+ziTeTi)eTeΦ̃0(θ)=(Mi2Me2)[R(θ)2R2]2R02.
(9)

The poloidally varying potential results from the density asymmetry created by the non-uniform distribution of ions around a flux surface due to the centrifugal force and pulls the electrons with the ions. Note that it is quadratic in the main ion Mach number. Regarding the input density, we find it convenient for contact with experiment to use the values at the outboard midplane. Note that, without loss of generality, we can rewrite the density at a given flux surface ψ as

na(θ)=na(θ0)exp{zaeTaΦ*(θ)+Ma2[R(θ)2R(θ0)22R02]},
(10)

where Φ*(θ)Φ0(θ)Φ0(θ0) and we take θ0=0. Thus, given na(θ0), Eq. (8) can be solved for Φ*(θ). In the diamagnetic limit (Ma = 0), Φ*=0.

The magnetostatic equilibrium is also altered by the rotation. Momentum balance and Ampère's law determine the basic set of magnetostatic equilibrium equations. We can write these as

1cJ×B=(p)eff,
(11)
×B=4πcJ,
(12)

where the effective equilibrium pressure gradient is defined as

(p)effp+anamaU0·U0.
(13)

Here, p=anaTa is the total equilibrium kinetic plasma pressure. Note that this is poloidally varying due to the poloidal variation of the equilibrium density. While approximations such as trace impurities and adiabatic electrons are often made regarding the total pressure, in CGYRO, the exact total pressure is used, and thus, the magnetostatic equilibrium is fully consistent with the equilibrium density and potential used in the kinetic equations. The effective pressure gradient modifies not only the gyrokinetic drift velocity (as will be detailed shortly) but also the local Grad-Shafranov equilibrium calculation.21 These effects may be important in the KBM-dominated plasma edge.

The gyrokinetic equation yields the O(ρ*) fluctuating distribution in terms of the distribution function of gyrocenters, ha(R). We write the equation in (ξ,ua) velocity-space coordinates which are defined with respect to the rotating velocity frame described in Sec. II A. Here, ξv/v is the cosine of the pitch angle and uav/(2vta) is the normalized velocity. For the spatial dimensions (r,α,θ), we write the distribution function ha(R,ξ,ua) in terms of the discrete indices (p, n) which are the Fourier wave numbers corresponding to the (r,α) coordinates, such that

ha(x,α,θ)=f0an=NNp=MMh̃a(n,p,θ)eipxeinα,
(14)

where x2π(rr0)/L[0,2π) and L is the radial domain size. In this notation, introducing unit vectors bB/B, exr/|r| and eyb×ex, we have the corresponding wavenumbers

kx(n,p)=ex·k=2πprLnr·ν|r|,
(15)
ky(n,p)=ey·k=nB|ψ|,
(16)

where k2=kx2+ky2. The related poloidal wavenumber kθ is defined so that it is a proper flux-surface function: kθnq/r, where q(ψ) is the safety factor. The field potential Ψa(R)15 can also be expanded as

eTeΨa(x,α,θ)=n=NNp=MMΨ̃a(n,p,θ)eipxeinα,
(17)

where

Ψ̃a=J0(γa)(δϕ̃effvcδÃ)+v2ΩcacJ1(γa)γaδB̃.
(18)

Here, γakv/Ωca. Note that we have defined a new effective electrostatic potential as

δϕ̃eff=δϕ̃U0c·δÃ,
(19)

which is just the potential in the rotating frame, using the Lorentz transformation for the electric field: δEeff=E+U0/c×B. The nonadiabatic distribution H̃a is defined in terms of Ψ̃a as

H̃a=h̃a+zaTeTaΨ̃a.
(20)

We write the gyrokinetic equation in dimensionless form by normalizing to the deuteron sound speed cs=Te/mD (where Te is the electron temperature and mD is the deuteron mass) and the midplane minor radius of the last closed flux surface, a. Thus, we define the normalized time as τ(cs/a)t. The dimensionless gyrokinetic equation for h̃a(n,p,θ,ξ,ua) becomes

h̃aτiΩsXh̃ai(Ωθ+Ωξ+Ωd)H̃aiΩ*Ψ̃a+2πaLqρsrh̃aΨ̃a=acsbCabL(H̃a,H̃b).
(21)

In Eq. (21), CabL is the linearized gyrophase-averaged collision operator.22 The spectral representation of the nonlinear term h̃aΨ̃a (a generalized convolution) is given by Eq. (56) of Ref. 15. The dimensionless operators are

streaming

iΩθ=vwsθ,
(22)

trapping

iΩξ=vtawsua2(1ξ2)lnBθξiΩξ,cf,
(23)

drift

iΩd=avtacsb×[ua2(1+ξ2)BB+cBvtaΦ(θ0)+ua2ξ28πB2(p)eff]·ikρaiΩd,coriΩd,cf,
(24)

gradient drive

iΩ*=[aLna+aLTa(ua232)azaeTadΦ0(θ0)dr]ikθρsiΩ*,trsiΩ*,cf.
(25)

In Eqs. (22)–(25), the species gyroradius is defined as ρavtamac/(zaeB), with the effective ion-sound gyroradius given by ρscsmDc/(eBunit), where Bunit(r)=(q/r)ψ is the effective magnetic field. We have also introduced an effective velocity wscs(JψB)/a and the density and temperature gradient inverse scale lengths, 1/Lna=dlnna(θ0)/dr and 1/LTa=dlnTa/dr.

Note that the rotation terms have been written explicitly in Eqs. (22)–(25). We first define the weak rotation terms, which are linear in Ma,

Coriolis drift

iΩd,cor=2MavcsaR0b×(RJψBRθφBtBR)·ikρa,
(26)

rotation shear drive

iΩ*,trs=γpavtavvtaRBtR0Bikθρs.
(27)

Here, γp=R0dω0/dr is the toroidal rotation shearing rate.23 Note that γp is sometimes referred to by the misnomer “parallel velocity shearing rate,” which reflects the historical usage in sheared slab geometry for which γp=R0d(U/R)/dr.24 The sheared E × B flow Ωs in Eq. (21) is also a weak rotation term, being linear in Ma. For inclusion of sheared E × B flow in Eq. (21), we first define the dimensionless E × B shearing rate

Ωsacsγs=kθL2πacsγE,
(28)

where γE=(r/q)dω0/dr is the Waltz shearing rate1 and γs is the domain shearing rate. In Eq. (21), the total shearing term is iΩsXh̃a, where X is a spectral shearing operator, defined by Eq. (21) of Ref. 25.

For the case of sonic rotation, the centrifugal terms must also be considered. These terms are quadratic in Ma and contribute to the trapping operator through

iΩξ,cf=12uaλaθ[vwsua+2vtaws(1ξ2)ξ],
(29)

the drift operator through

iΩd,cf=vtacsb×[Ma2aRRR02+acBvtaΦ*]·ikρa,
(30)

and the gradient drive operator through

iΩ*,cf={aLTa[zaeTaΦ*Ma22R02(R2R(θ0)2)]+Ma2aR(θ0)R02dR(θ0)dr+MaγpavtaR02(R2R(θ0)2)}ikθρs.
(31)

Note that centrifugal effects also modify the equilibrium density in the equilibrium distribution f0a in ha(x,α) in Eq. (14) and in the effective pressure gradient in the drift in Eq. (24). It is common in gyrokinetics to neglect the centrifugal effects due to their complexity and to instead only include the weak rotation terms, namely, the Coriolis drift, parallel velocity shear drift, and sheared E × B flow. The implications of this reduction will be explored in detail shortly.

The dimensionless, spectral gyrocenter-Maxwell equations are given by

(k2λD2+aza2TeTad3vf0ane)δϕ̃eff=azad3vf0aneJ0(γa)H̃a,
(32)
2βe,unitk2ρs2δÃ=azad3vf0anevcsJ0(γa)H̃a,
(33)
2βe,unitBBunitδB̃=ad3vf0anemav2TeJ1(γa)γaH̃a.
(34)

Here, we have introduced the Debye length λD=Te/(4πnee2) and the effective electron beta βe,unit=8πneTe/Bunit2. Note that with sonic rotation, we can solve the Maxwell equations for the effective electrostatic potential δϕ̃eff and thus the field potential given by Eq. (18). In the limit of weak rotation, δϕ̃effδϕ̃.

Finally, the turbulent radial particle flux Γa, energy flux Qa, and momentum flux Πa are defined in accordance with Ref. 8 as

Γa=n,pd3vH̃a*w̃a1,
(35)
Qa=Tan,pd3vH̃a*w̃a1(ua2+λa),
(36)
Πa=Rman,pd3vH̃a*[w̃a1(ω0R+BtBv)+w̃a2vta],
(37)

where

w̃a1=ikθρscsΨ̃aandw̃a2=kxkθk2BpBunitzaTeTavtaχ̃a.
(38)

The field quantity χ̃a is given by

χ̃a=γaJ1(γa)(δϕ̃effvcδÃ)v2γaΩcac[J1(γa)γaJ0(γa)]δB̃.
(39)

Note that the (n=0,p=0) component of the distribution, h̃a(0,0,θ), is zero in gyrokinetic theory. Thus, spatial averages of the distribution, such as the gyrokinetic bootstrap current, are zero.

To recap, CGYRO solves a 5 D nonlinear, electromagnetic gyrokinetic equation, Eq. (21), for each species a. These equations couple through the gyrocenter Maxwell equations, Eqs. (32)–(34), and also through the collision operator. As described in Ref. 15, the solver is both pseudospectral in velocity space dimensions (ξ,ua) and spectral in the radial and binormal spatial dimensions (x,α) and uses a 6th-order finite difference scheme in θ (the fieldline direction) with a 6th-order conservative dissipation algorithm.

The simulations use the General Atomics (GA) standard case parameters,26 considering a single ion species (deuterium) and electrons, with mi/me=60 with the following: R0/a=3,r/a=0.5, q = 2, s = 1, Ti = Te, a/Lni=a/Lne=1,a/LTi=a/LTe=3. In CGYRO, we specify the density and the density gradient at the outboard midplane, θ0=0. For the geometry, we use an unshifted Miller local circular equilibrium.27 We emphasize that the local equilibrium is not a model equilibrium, but rather an exact solution of the Grad-Shafranov equations on the specified flux-surface. The local equilibrium method is complicated to implement and the reader should consult Ref. 21.

To simulate a single linear eigenmode, we consider fixed kθρs=0.3. To vary the rotation frequency ω0 and rotation shear dω0/dr, we (equivalently) specify the main ion Mach number Mi=ω0R0/vti and the toroidal rotation shearing rate γp=R0dω0/dr. The species-dependent Mach numbers Ma are varied consistently with Mi. The poloidally varying equilibrium density, computed by CGYRO, is shown in Fig. 1 for varying Mi. Accumulation of ions at the outboard midplane (θ = 0) is observed as Mi increases due to the centrifugal force. To assess the impact of rotation, three cases are considered:

  1. Sonic rotation: All sonic rotation terms are included in the equilibrium and gyrokinetic equations.

  2. Weak rotation: Only rotation terms which are linear in Ma are included in the gyrokinetic equations; specifically, the Coriolis drift Ωd,cor, the toroidal rotation shear drive Ω*,trs, and the E × B term Ωs. No centrifugal terms are included, meaning Φ*=0 and na=na(0).

  3. Only centrifugal terms: Only rotation terms due to the centrifugal force are included in the gyrokinetic equations; specifically, Ωξ,cf,Ωd,cf, and Ω*,cf. The poloidal variation of the equilibrium density is included in f0a as well as in the pressure-gradient drift (p)eff. Weak rotation terms are neglected.

FIG. 1.

Poloidal variation of the ion equilibrium density relative to the value at the outboard midplane (θ = 0) for varying values of the ion Mach number Mi for the case of a pure plasma.

FIG. 1.

Poloidal variation of the ion equilibrium density relative to the value at the outboard midplane (θ = 0) for varying values of the ion Mach number Mi for the case of a pure plasma.

Close modal

As we have explained previously, the limit of weak rotation is commonly assumed in linear and nonlinear gyrokinetic simulations. It is equivalent to retaining only terms linear in the Mach number Ma. While this approximation is typically very accurate for deuterium ions at low rotation, it is not generally valid for heavy impurities.

To test the implementation of all sonic terms, CGYRO has been benchmarked with GKW11 using the GA standard case parameters specified earlier. The results are summarized in Figs. 2 and 3. The results show that the intercode agreement is excellent, without any significant deviation, in scans over Mi and γp in both low- and high-rotation regimes.

FIG. 2.

Linear growth rate γ and real frequency ωr versus ion Mach number Mi at fixed rotation shear γp = 0. The results from CGYRO are compared with those from GKW. Both codes include full sonic rotation.

FIG. 2.

Linear growth rate γ and real frequency ωr versus ion Mach number Mi at fixed rotation shear γp = 0. The results from CGYRO are compared with those from GKW. Both codes include full sonic rotation.

Close modal
FIG. 3.

Linear growth rate γ and real frequency ωr versus rotation shear γp at fixed ion Mach number Mi=0.4. The results from CGYRO are compared with those from GKW. Both codes include full sonic rotation.

FIG. 3.

Linear growth rate γ and real frequency ωr versus rotation shear γp at fixed ion Mach number Mi=0.4. The results from CGYRO are compared with those from GKW. Both codes include full sonic rotation.

Close modal

The effect of rotation on gyrokinetic linear stability of ion-temperature-gradient (ITG) modes and trapped electron modes (TEM) using the GA standard case parameters is further studied in Fig. 4. The results with sonic rotation are compared with those including only weak rotation terms and those including only centrifugal terms. At low Mi, the ITG mode is mostly affected by the stabilizing Coriolis drift Ωd,cor, while centrifugal terms are insignificant (i.e., red curve follows the blue curve). However, as Mi increases, the ITG mode transitions to a TEM, which is further destabilized as Mi increases. For the TEM, the centrifugal terms (green curve) dominate the Coriolis drift term (blue curve). Although the centrifugal force itself is small for the electrons due to their small mass, enhancement of the mirror force due to the induced poloidally varying electrostatic potential leads to an increased effective electron trapped fraction. For these scans, a resolution of (Nθ=Nξ=24) is required to resolve the enhanced particle trapping. This trapping effect has a strong destabilizing influence on the TEM, causing the mode transition to occur. This has been previously shown with the GKW code for similar parameters.12 However, the effect on TEMs can be quickly reduced by the detrapping of electrons due to collisions, as shown in Fig. 5. The results use the Sugama collision model. As the collision frequency increases, the growth rate of the TEM decreases and eventually mode transition to a dominant ITG mode occurs.

FIG. 4.

Linear growth rate γ and real frequency ωr versus main ion Mach number Mi at fixed rotation shear γp=0. The results with the complete sonic formulation are compared with those including only weak rotation terms and those including only centrifugal rotation terms.

FIG. 4.

Linear growth rate γ and real frequency ωr versus main ion Mach number Mi at fixed rotation shear γp=0. The results with the complete sonic formulation are compared with those including only weak rotation terms and those including only centrifugal rotation terms.

Close modal
FIG. 5.

Linear growth rate γ and real frequency ωr versus collision rate τee1 for a TEM at fixed rotation shear γp=0. Three different values of the main ion Mach number Mi are considered. The sonic rotation formulation was used in all cases.

FIG. 5.

Linear growth rate γ and real frequency ωr versus collision rate τee1 for a TEM at fixed rotation shear γp=0. Three different values of the main ion Mach number Mi are considered. The sonic rotation formulation was used in all cases.

Close modal

In contrast to the effect of the main ion Mach number, Fig. 6 shows that increasing the toroidal rotation shear γP destabilizes the ITG mode. Comparison with the Mi = 0 case indicates that, at large γp, the Miγp cross term explicitly in Ω*,cf and implicitly in Φ* in Ωd,cf enhances the destabilizing effect. In contrast to the ITG mode, the TEM is stabilized by the toroidal rotation shear, leading to TEM to ITG mode transition.

FIG. 6.

Linear growth rate γ and real frequency ωr versus rotation shear γp at various low (a-b) and high (c-d) main ion Mach numbers Mi. The results include full sonic rotation. The low Mi modes are ITG, while the high Mi modes transition from TEM to ITG as γp increases.

FIG. 6.

Linear growth rate γ and real frequency ωr versus rotation shear γp at various low (a-b) and high (c-d) main ion Mach numbers Mi. The results include full sonic rotation. The low Mi modes are ITG, while the high Mi modes transition from TEM to ITG as γp increases.

Close modal

To study the effects of rotation on turbulent transport, we consider the GA standard case parameters used above with the addition of a trace impurity species. For reactors, transport of heavy impurities from the material surface into the high temperature core must be minimized since accumulation can lead to fuel dilution and radiation losses. Both fully ionized carbon (zi = 6) and partially ionized tungsten (zi = 41) are considered. The trace limit is assumed, with nimp/ne=105. The impurity species is assumed to be in thermal equilibrium with the main ions. We consider a typical fixed rotation frequency of Mi=0.2 and varying rotation shear. The results compare inclusion of full sonic rotation effects with those including only weak rotation terms, neglecting centrifugal effects. The sonic and weak rotation cases converge when Ma = 0. The poloidal variation of the equilibrium density for the ions in the full sonic rotation case is shown in Fig. 7. Despite the small rotation frequency, the density asymmetry for the impurities is large, especially for tungsten. This is not surprising since MC=0.49 and MW=1.92.

FIG. 7.

Poloidal variation of the deuterium, trace carbon, and trace tungsten equilibrium densities relative to the value at the outboard midplane (θ = 0) for Mi=0.2. This corresponds to MC=0.49 and MW=1.92.

FIG. 7.

Poloidal variation of the deuterium, trace carbon, and trace tungsten equilibrium densities relative to the value at the outboard midplane (θ = 0) for Mi=0.2. This corresponds to MC=0.49 and MW=1.92.

Close modal

The results for the critical impurity density gradient scale length are shown in Fig. 8, comparing carbon and tungsten. The critical impurity density gradient length scale is defined as that for which the quasilinear turbulent radial particle flux is zero. For carbon, the case with only weak rotation terms (blue curve) qualitatively follows the trend of the full rotation case (red curve), with the latter predicting a larger critical gradient. The dominant rotation effect is the toroidal rotation shear drive Ω*,trs. The critical impurity density gradient decreases as γp increases, such that the impurity density profile is peaked for small rotation shear, γpa/cs<4, and hollow for large rotation shear. In contrast to carbon, for tungsten, including only weak rotation terms while neglecting the centrifugal terms produces incorrect results. This is not surprising considering the large poloidal density asymmetry due to the centrifugal force shown in Fig. 7. For γpa/cs>2, the case with inclusion of only weak rotation terms (blue curve) predicts a hollow impurity density profile, similar to the case with carbon for large γp. However, inclusion of sonic rotation (red curve) yields a critical impurity density profile which is peaked for the entire range of γp. The fact that both the weak and sonic rotation curves are close to the Ma = 0 limit (green curve) at γp=0 indicates that the main effect is due to the MaγP cross-terms rather than the Ma2 terms.

FIG. 8.

Critical impurity density gradient scale length a/Ln for which the quasi-linear turbulent radial particle flux is zero for carbon (a) and tungsten (b) versus rotation shear γp at Mi = 0.2. The results with full sonic rotation are compared with those including only weak rotation terms. The sonic and weak rotation cases converge when Ma = 0.

FIG. 8.

Critical impurity density gradient scale length a/Ln for which the quasi-linear turbulent radial particle flux is zero for carbon (a) and tungsten (b) versus rotation shear γp at Mi = 0.2. The results with full sonic rotation are compared with those including only weak rotation terms. The sonic and weak rotation cases converge when Ma = 0.

Close modal

Figures 9 and 10 show the corresponding nonlinear turbulent fluxes for a/LnW=0. The particle and energy fluxes are normalized by the gyroBohm factors

ΓGB=necsρ*2andQGB=neTecsρ*2
(40)

and by the relative density, such that {Γa,Qa}norm={Γa,Qa}/(na/ne). The deuterium ions and electrons are negligibly affected by the centrifugal force, and thus, the full sonic and weak rotation results are very similar for these species. However, similar to the quasilinear results, the tungsten particle flux is strongly affected by the centrifugal dynamics and thus including only weak rotation terms gives a large error. With full rotation, the tungsten particle flux is radially inward, with its magnitude increasing strongly with increasing rotation shear. This indicates a strong inward pinch, which leads to detrimental core accumulation in a reactor. Including only weak rotation terms incorrectly predicts an outward particle flux at large rotation shear, as shown by the blue curve in Fig. 10(b). Note that inclusion of E × B shear (not shown here) has a negligible effect on these results. Figure 10(a) shows that the tungsten energy flux is destabilized by the rotation shear, in contrast with the slight stabilizing effect observed for the deuterium ions and electrons in Fig. 9. The insignificance of centrifugal dynamics indicates that, unlike the tungsten particle flux, the tungsten energy flux is dominated by the toroidal rotation shear drive Ω*,trs. For comparison, the results for the carbon nonlinear turbulent fluxes for a/LnC=0 are shown in Fig. 11. As expected, the carbon fluxes are less affected by the centrifugal dynamics than tungsten but more strongly affected than the deuterium ions. Similar to tungsten, an inward pinch is predicted with full sonic rotation, and neglecting centrifugal dynamics incorrectly predicts an outward flux, as shown by the blue curve in Fig. 11(b). Ultimately, these results indicate that inclusion of full sonic rotation, including centrifugal effects, is critical for heavy impurities, such as tungsten, as well as for moderate impurities like carbon at large rotation shear.

FIG. 9.

Nonlinear turbulent radial energy flux and particle flux for deuterium ions (D) and electrons (e) versus rotation shear γp at Mi=0.2. The results with full sonic rotation are compared with those including only weak rotation terms.

FIG. 9.

Nonlinear turbulent radial energy flux and particle flux for deuterium ions (D) and electrons (e) versus rotation shear γp at Mi=0.2. The results with full sonic rotation are compared with those including only weak rotation terms.

Close modal
FIG. 10.

Nonlinear turbulent radial energy flux and particle flux for tungsten versus rotation shear γp at Mi=0.2. The results with full sonic rotation are compared with those including only weak rotation terms.

FIG. 10.

Nonlinear turbulent radial energy flux and particle flux for tungsten versus rotation shear γp at Mi=0.2. The results with full sonic rotation are compared with those including only weak rotation terms.

Close modal
FIG. 11.

Nonlinear turbulent radial energy flux and particle flux for carbon versus rotation shear γp at Mi=0.2. The results with full sonic rotation are compared with those including only weak rotation terms.

FIG. 11.

Nonlinear turbulent radial energy flux and particle flux for carbon versus rotation shear γp at Mi=0.2. The results with full sonic rotation are compared with those including only weak rotation terms.

Close modal

Finally, Fig. 12 compares the impurity turbulent and neoclassical fluxes. The neoclassical fluxes are computed with NEO,18 which also includes full sonic rotation. For low collision frequency typical in the tokamak core (ν¯e=0.1), the neoclassical fluxes are much smaller than the turbulent fluxes, by 1–2 orders of magnitude. However, since the magnitude of the neoclassical fluxes increases with increasing ν¯e while the magnitude of the turbulent fluxes decreases, at ν¯e = 1, which is typical in the H-mode pedestal, the neoclassical fluxes become competitive with the turbulent fluxes. In fact, for γpa/cs<2, the neoclassical fluxes dominate. Thus, both must ultimately be considered for accurate transport modeling of tungsten. The neoclassical tungsten particle fluxes are also directed inward, enhancing the pinch effect.

FIG. 12.

Radial particle flux for tungsten versus rotation shear γp at Mi = 0.2, comparing the nonlinear turbulent flux from CGYRO with the neoclassical flux from NEO at various collision frequencies, ν¯e=(a/cs)τee1. Both codes include full sonic rotation.

FIG. 12.

Radial particle flux for tungsten versus rotation shear γp at Mi = 0.2, comparing the nonlinear turbulent flux from CGYRO with the neoclassical flux from NEO at various collision frequencies, ν¯e=(a/cs)τee1. Both codes include full sonic rotation.

Close modal

The impact of sonic toroidal rotation on gyrokinetic stability is studied using the CGYRO code, which solves the gyrokinetic equation using the sonic rotation formalism of Sugama.8 This formalism is based on a rigorous asymptotic expansion of the Fokker-Planck equation and is valid for flow velocities on the order of the ion thermal speed. Centrifugal drifts and electrostatic trapping dynamics are included, in addition to the standard weak rotation terms included in most gyrokinetic codes: the E × B flow, Coriolis drift, and toroidal rotation shear. Because of their complexity, these sonic rotation effects—which are quadratic in the Mach number—are rarely considered in gyrokinetic simulations.

For a pure plasma, the ITG mode is primarily affected by the stabilizing Coriolis drift for low ion Mach number (Mi<1). At large Mi, the dominant effect of sonic rotation is on the electron dynamics. Though the centrifugal force is small for electrons due to their small mass, enhancement of the mirror force leads to an increased effective electron trapped fraction. This has a strong destabilizing influence on TEMs and leads to mode transition from ITG to TEM for increasing Mi. However, the effect on TEMs can be quickly reduced by the detrapping of electrons due to collisions. In contrast, the toroidal rotation shear γp is destabilizing for the ITG mode and stabilizing for the TEM, leading to a TEM-to-ITG mode transition. For heavy impurity turbulent transport at moderate Mi, including only the weak rotation terms (i.e., E × B flow, Coriolis drift, and toroidal rotation shear) while neglecting the centrifugal terms results in a large error. For analysis of the quasilinear critical tungsten density gradient, inclusion of only weak rotation terms predicts a hollow impurity density profile, while the complete sonic theory yields a critical impurity density profile which is peaked due to the MaγP cross-terms. In contrast, carbon is not significantly affected by rotation, and a hollow impurity density profile is predicted at large rotation shear.

Nonlinear turbulent fluxes for a/LnW=0 similarly show that the sonic formulation must be used, yielding a strong inward pinch, which can lead to detrimental impurity accumulation in the reactor core. Including only weak rotation terms incorrectly predicts an outward particle flux at large rotation shear. At a large collision rate, the neoclassical particle fluxes become competitive with the turbulent fluxes, enhancing the pinch effect. NEO simulations in conjunction with quasilinear GKW calculations have previously been used to analyze hybrid scenarios in JET. In this work, it was found that the neoclassical transport rather than the turbulent transport dominates the central part of the plasma due to enhancement of the neoclassical diffusive and convective components by rotation.28 More recent results, however, indicate that quasilinear simulations may be inaccurate for mixed turbulence regimes with comparable ion and electron transport.14 Unfortunately, impurity transport modeling that couples neoclassical with nonlinear gyrokinetic simulation and incorporates heavy impurity sources and radiative losses is still computationally impractical. Thus, it is important that reduced models of transport29,30 be extended to treat the sonic regime correctly in order to provide realistic assessments of heavy impurity accumulation.

This material is based on work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, Theory program, under Award No. DE-FG02-95ER54309, and by the Edge Simulation Laboratory project under Grant No. DE-FC02-06ER54873. This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of the authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.

1.
R.
Waltz
,
G.
Kerbel
, and
J.
Milovich
,
Phys. Plasmas
1
,
2229
(
1994
).
2.
K.
Burrell
,
Phys. Plasmas
4
,
1499
(
1997
).
3.
G.
Staebler
,
R.
Waltz
, and
J.
Kinsey
,
Phys. Plasmas
18
,
056106
(
2011
).
4.
R.
Waltz
,
G.
Staebler
,
J.
Candy
, and
F.
Hinton
,
Phys. Plasmas
14
,
122507
(
2007
).
5.
A.
Peeters
,
C.
Angioni
, and
D.
Strintzi
,
Phys. Rev. Lett.
98
,
265003
(
2007
).
6.
A.
Brizard
,
Phys. Plasmas
2
,
459
(
1995
).
7.
H.
Sugama
and
W.
Horton
,
Phys. Plasmas
4
,
2215
(
1997
).
8.
H.
Sugama
and
W.
Horton
,
Phys. Plasmas
5
,
2560
(
1998
).
9.
A.
Peeters
,
D.
Strintzi
,
Y.
Camenen
,
C.
Angioni
,
F.
Casson
,
W.
Hornsby
, and
A.
Snodin
,
Phys. Plasmas
16
,
042310
(
2009
).
10.
P.
Helander
and
D.
Sigmar
,
Collisional Transport in Magnetized Plasmas
(
Cambridge University Press
,
Cambridge
,
2002
).
11.
A.
Peeters
,
Y.
Camenen
,
F.
Casson
,
W.
Hornsby
,
A.
Snodin
,
D.
Strintzi
, and
G.
Szepesi
,
Comput. Phys. Commun.
180
,
2650
(
2009
).
12.
F.
Casson
,
A.
Peeters
,
C.
Angioni
,
Y.
Camenen
,
W.
Hornsby
, and
A.
Snodin
,
Phys. Plasmas
17
,
102305
(
2010
).
13.
C.
Angioni
,
F.
Casson
,
C.
Veth
, and
A.
Peeters
,
Phys. Plasmas
19
,
122311
(
2012
).
14.
C.
Angioni
,
R.
Bilato
,
F.
Casson
,
E.
Fabel
,
P.
Mantica
,
T.
Odstrcil
,
M. V. A. U.
Team
, and
J.
Contributors
,
Nucl. Fusion
57
,
022009
(
2017
).
15.
J.
Candy
,
E.
Belli
, and
R.
Bravenec
,
J. Comput. Phys.
324
,
73
(
2016
).
16.
E.
Belli
and
J.
Candy
,
Plasma Phys. Controlled Fusion
50
,
095010
(
2008
).
17.
E.
Belli
and
J.
Candy
,
Plasma Phys. Controlled Fusion
51
,
075018
(
2009
).
18.
E.
Belli
and
J.
Candy
,
Plasma Phys. Controlled Fusion
54
,
015015
(
2012
).
19.
F.
Hinton
and
S.
Wong
,
Phys. Fluids
28
,
3082
(
1985
).
20.
M.
Kruskal
and
R.
Kulsrud
,
Phys. Fluids
1
,
265
(
1958
).
21.
J.
Candy
,
Plasma Phys. Controlled Fusion
51
,
105009
(
2009
).
22.
E.
Belli
and
J.
Candy
,
Plasma Phys. Controlled Fusion
59
,
045005
(
2017
).
23.
R.
Waltz
,
G.
Kerbel
,
J.
Milovich
, and
G.
Hammett
,
Phys. Plasmas
2
,
2408
(
1995
).
24.
J.
Dong
and
W.
Horton
,
Phys. Fluids B
5
,
1581
(
1993
).
25.
J.
Candy
and
E.
Belli
,
J. Comput. Phys.
356
,
448
(
2018
).
26.
R.
Waltz
,
G.
Staebler
,
W.
Dorland
,
G.
Hammett
,
M.
Kotschenreuther
, and
J.
Konings
,
Phys. Plasmas
4
,
2482
(
1997
).
27.
R.
Miller
,
M.
Chu
,
J.
Greene
,
Y.
Lin-liu
, and
R.
Waltz
,
Phys. Plasmas
5
,
973
(
1998
).
28.
C.
Angioni
,
P.
Mantica
,
T.
Putterich
,
M.
Valisa
,
M.
Baruzzo
,
E.
Belli
,
P.
Belo
,
F.
Casson
,
C.
Challis
,
P.
Drewelow
,
C.
Giroud
,
N.
Hawkes
,
T.
Hender
,
J.
Hobirk
,
T.
Koskela
,
L. L.
Taroni
,
C.
Maggi
,
J.
Mlynar
,
T.
Odstrcil
,
M.
Reinke
,
M.
Romanelli
, and
JET Contributors
,
Nucl. Fusion
54
,
083028
(
2014
).
29.
G.
Staebler
,
J.
Kinsey
, and
R.
Waltz
,
Phys. Plasmas
12
,
102508
(
2005
).
30.
C.
Bourdelle
,
X.
Garbet
,
F.
Imbeaux
,
A.
Casati
,
N.
Dubuit
,
R.
Guirlet
, and
T.
Parisot
,
Phys. Plasmas
14
,
112501
(
2007
).