Table Of ContentGEM – An Energy Conserving Electromagnetic Gyrofluid Model
Bruce D. Scott∗
Max-Planck-Institut fu¨r Plasmaphysik, Euratom Association, D-85748 Garching, Germany
(Dated: February 2, 2008)
The details of fluctuation free energy conservation in the gyrofluid model are examined. The
polarisation equation relates ExB flow and eddy energy to combinations of the potential and the
5 density and perpendicular temperature. These determine the combinations which must appear
0 underderivativesinthemomentequationssothatnotonlythermalfreeenergybutitscombination
0 withtheExBenergyisproperlyconservedbyparallelandperpendicularcompressionaleffects. The
2
resultingsystemexhibitsthesamequalitativeenergytransferpropertiesascorrespondingBraginskii
n or Landau fluid models. One clear result is that the numerical model built on these equations
a is well behaved for arbitrarily large perpendicular wavenumber, allowing exploration of two scale
J phenomenalinkingdynamics at theion and electron gyroradii. When thenumerical formulation is
3 done in the globally consistent flux tube model, the results with adiabatic electrons are consistent
2 with the“Cyclone Base Case” results of gyrokinetic models.
] PACSnumbers: 52.65.Tt,52.35.Ra,52.30.-q,52.25.Fi
h
p
-
m
s I. INTRODUCTION – GYROFLUID ENERGY CONSERVATION IN GENERAL
a
l
p
Gyrofluid models, whose most prominent application has been to tokamak core turbulence as exemplified by the
.
s Cyclone Base Case [1], were originally constructed to incorporate finite ion gyroradius effects at arbitrary order into
c
simplecomputationsofturbulenceoccurringinlargelytwodimensionalfluidexperiments[2]. Totreationtemperature
i
s gradient (ITG) turbulence the temperatures were incorporated and the model acquired severalnew advection terms,
y
producing nonlinearities as well as drift frequency corrections, resulting from the effect of temperature fluctuations
h
on the gyroaveraging operator [3, 4]. However, although the nonlinearities were incorporated, much of the analysis
p
[ involvedlinearfrequenciesandgrowthratesandthenonlinearitieswereaddedlargelyasanafterthought. Specifically,
there was no complete energetic analysis of the type familiar from drift wave turbulence work [5, 6, 7, 8]. The GEM
1
modelwasintroducedpreviouslyinthecontextofenergeticconsiderations,includingthecorrespondencebetweenfluid
v
drift and gyrofluid models under drift ordering [9]. We develop the energetics for the six-moment model including
4
temperatureandparallelheatfluxdynamicsherein,placingthismodelonasecureenergeticfootingforthefirsttime.
2
1 Underdriftordering[10],thisenergeticsinvolves“fluctuationfreeenergy,”inwhichthethermalfreeenergyentersas
1 the averagesquaredamplitude ofthe density fluctuations [5], with additionalcontributions from the averagesquared
0 amplitude of the temperature fluctuations in the appropriately generalised models [7]. The rest of the free energy is
5 made up of contributions due to the ExB energy involving the electrostatic potential, the magnetic energy involving
0
the parallel magnetic potential, and the parallel free energy in not only the parallel velocities but also the parallel
/
s heat fluxes [8]. A properly constructed model should conserve this free energy in all processes except those involving
c
clearly identifiable sources (gradients) and sinks (dissipative processes such as resistivity, thermal conduction, or
i
s Landau damping). Many models neglect this consideration because small errors in the energetics lead simply to
y
negligible contributions to the growth rate in a linear model. However, if the model is to be useful in a turbulent
h
setting, the energetics must be consistent in order to achieve a reliable saturated state in which the salient energetic
p
: processes and the turbulent transport can be statistically measured.
v
The processes linking the thermal free energy to the ExB turbulence are of physical importance because the free
i
X energysourcesandsinksarenotintheperpendicularequationofmotion. Theenergysourcegivenbythegeneralprofile
gradient, is in the equations for the thermal state variables (density, temperatures). The dissipation processes are in
r
a the equations for the parallel flux variables (current, heat fluxes). In saturation, the ExB energy itself is maintained
as a statistical balance between various conservative transfer effects which connect to these sources and sinks in the
otherpartsofthe dynamics(this neglects certainrotationdamping processeswhichareoftennotconsidered). Inthis
balance it is important that the conservative nature of such processes such as shear Alfv´en dynamics or interchange
effects is maintained in the model. The simplest example is an isothermal two dimensional magnetohydrodynamical
∗email: [email protected];
URL:http://www.rzg.mpg.de/~bds/
2
(MHD) interchange model, given by
n M c2 ∂
i i 2
φ= T (n ) (1)
B2 ∂t∇⊥ − eK e
e e
∂n
e +v (n +n )=n (φ) (2)
E e e e
∂t ·∇ K
e
where the perpendicular Laplacian and curvature operator areedefined bye
2 = [b (b )] (3)
∇⊥ −∇· × ×∇
c
= (B ) (4)
K ∇· B2 ×∇
h i
respectively, with B=Bb the equilibrium magnetic field. The ExB velocity is given by
c
v = b φ (5)
E
B ×∇
The curvature terms in Eqs. (1,2) respectively represent quasistaetic compression of the diamagnetic current and the
ExB velocity.
Under drift ordering [10], the background parameters are constants except where operated upon by v , the
E
·∇
ExB advection. The magnetic field is treated as constant except for the existence of (). In the advection term, the
K
velocity v is treated as divergence free; the finite ExB divergence is accounted for by (φ). If we multiply these
E
K
equationsby φand(T /n )n , respectively,andintegrateoverthe entirespatialdomain,we find, neglectingsurface
e e e
−
terms, e
e e ∂ n M c2 2
i i
dΛ φ = dΛT n φ (6)
∂t 2 B2 ∇⊥ − e eK
Z (cid:12) (cid:12) Z (cid:16) (cid:17)
(cid:12)(cid:12) e(cid:12)(cid:12) e e
2
∂ n T n
dΛ e e e = dΛT n v logn + dΛT n φ (7)
e e E e e e
∂t 2 n − ·∇ K
Z (cid:18) e(cid:19) Z Z (cid:16) (cid:17)
e
where dΛbyitselfgivesthetotalvolume. ThesetwolinesegivetheevolutionoftheEexBdrieftandthermalfreeenergy,
respectively. The interchange effect represented by () transfers free energy between these two pieces conservatively,
R K
as we would expect from a compressional process, because () is at once a total divergence and a first order linear
K
differential operator. The source is given by the advection of the background gradient, proportional to the average
flux. This model will saturate only if there is a loss process at the boundary, or some nonlinear effect enters to cause
the source to go to zero, or some additional effect is considered in the model by which an explicit sink term appears
to balance the source (a detailed analysis of how this functions in a three dimensional drift Alfv´en model is given in
Ref. [8]).
Inagyrofluidmodel,thereisnoequationforthevorticityexplicitlyinvolvingtimederivatives. Instead,atthesame
levelofsophisticationasabove,thereisanevolutionequationforeachofthespecies’gyrocenterdensities,asopposed
tospacedensities. Ineachdensityequation,thepotentialisgyroaveragedusingasuitableconvolutionoperatorwhich
is described by a kernel in Fourier space:
φG =G φ = Gk⊥φk⊥eik⊥·x (8)
k
(cid:16) (cid:17) X⊥
e e e
The usual form for Gk⊥ is Γ10/2(bi), which is an average of the single particle form J0(k⊥v⊥/Ωi) averaged over the
perturbed distribution function [3]. The argument is b = k2ρ2, where ρ is the thermal gyroradius and Ω is the
i ⊥ i i i
gyrofrequency. Unlessphenomenabelowtheionscaleρ areconsidered,electrongyroradiuseffectsareusuallyignored,
i
sothatGk =1forthem. Thepotentialisgivenbyapolarisationequationwhichequatesthetwospacedensities,each
⊥
given by a combination of the gyrocenter and polarisation densities, the latter involving the gyroscreened potential.
The equations appear as [2]:
∂n T
e +v (n +n )=n φ e (n ) (9)
E e e e e
∂t ·∇ K − e K
(cid:16) (cid:17)
e
e e e
3
∂n T
i +u (n +n )=n φ + i (n ) (10)
E i i i G i
∂t ·∇ K e K
(cid:16) (cid:17)
e
e e e
G(n ) e n
i 2 2 e
+ ρ φ= (11)
n T i∇⊥ n
i i e
e e
Gyroscreening is distinct from gyroaveraging; in general this eis given by another operator whose usual kernel in
k⊥-space is Γ0(bi)−1. The low-k⊥ form of this is just −bi, which in real space yields the ρ2i∇2⊥ form used in Eq. 11.
It is important to note however that the same operator is involved in the gyroaveragingof the potential (Eq. 8) as
in the conversion of the gyrocenter density, n , to the ion space density given by the left side of Eq. (11). For the
i
ions,theoperatorisG; fortheelectrons,itisunity. Consistentwiththis,the ExBadvectionofthe iondensityoccurs
with the gyroaveragedpotential, with velocity
e
c
u = b φ (12)
E G
B ×∇
while the electrons are advected by the “bare” version, the samee v as in Eq. (5). The curvature operator acts on
E
the total force potential of eachspecies, respectively,i.e., the diamagnetic velocity divergences are necessarilykept in
the model, which cannot take the MHD form if both pressures are to contribute to charge separation.
The point to be made here is the wayin which the polarisationequation, Eq. (11), determines the acceptable form
for many of the differential operators in the moment equations. The ExB energy and the electron and ion thermal
free energies are given by
n M c2 2 n T n 2 n T n 2
i i e e e i i i
U = dΛ φ U = dΛ U = dΛ (13)
E 2 B2 ∇⊥ e 2 n i 2 n
Z (cid:12) (cid:12) Z (cid:18) e(cid:19) Z (cid:18) i(cid:19)
respectively. Using Eq. (11) and the(cid:12)(cid:12)Herme(cid:12)(cid:12)itian property of G, we meay recast U as e
E
φ n φn
G i e
U =e e (14)
E
2 − 2
e e ee
which shows the ion and electron contributions separately. By similar means, the time derivatives are given by
∂U ∂n ∂n
E i e
=eφ eφ = dΛT n φ dΛT n φ (15)
G e e i i G
∂t ∂t − ∂t − K − K
Z (cid:16) (cid:17) Z (cid:16) (cid:17)
e e
e e e e e e
∂U
e = dΛT n v logn + dΛT n φ (16)
e e E e e e
∂t − ·∇ K
Z Z (cid:16) (cid:17)
e e e
∂U
i = dΛT n u logn + dΛT n φ (17)
i i E i i i G
∂t − ·∇ K
Z Z (cid:16) (cid:17)
From this we can see that, since both polarisatioen and thermal energy foreeachespecies now follow from its density
equation, the density and correspondinggyroaveragedpotential must appear together under spatialderivatives in all
conservative processes in order for the energetics to remain consistent. For the ions this is φ +(T /n e)n , while for
G i i i
the electrons it is φ (T /n e)n . Thus, for the ions, the combination φ +(T /n e)n multiplied by the curvature
e e e G i i i
−
term in Eq. (10) yields a term which vanishes under the integration, conserving the free eneergy, and theesame occurs
for the electrons upeon multiplyieng the curvature term in Eq. (9) by φ e(Te/nee)ne. e
−
This is relativelytrivialfor suchaone-momentgyrofluidmodel, butthe centralpointis clear: the gyrofluidclosure
arising from gyroaveraging must be done the same way in the polareisation and ein the moment equations. That is,
in this case, the operator G in the gyroaveragedpotential for a given species must be the same as the one operating
upon the corresponding density in the polarisation equation. In a model which incorporates the temperatures, this
extends to anothergyroaveragingoperatoracting upon the temperature, and the same one acting upon the potential
producing extra finite gyroradius advection terms, plus corrections to the temperature wherever the latter appears
under parallelgradients or curvature terms. This is much more involvedand will be treated in the next two sections,
one discussing the problems with the currently standard gyrofluid model, and the next one formulating the GEM
model. Then, the last section shows the applications of GEM to the damping of kinetic shear Alfv´en waves, to the
hyperfine electron gyroradius scale turbulence problem, and to the Cyclone Base Case which provides the standard
benchmark.
4
II. DRIFT ORDERING AND NORMALISATION CONVENTIONS
With many constant factors containing the physical units to be carried about while manipulating the equations,
it is convenient to go to a system of normalised units. While this is somewhat arbitrary, two examples are useful in
elucidating the salient units for gyrofluid turbulence. One of these is the prospect of force balance in the electrons
along the magnetic field. There are several effects which affect the evolution of the electric current in the electron
Ohm’s law, but these all act to mediate the response of the electrons to the two static parallel forces: the pressure
gradient and the static part of the parallel electric field, given by
e φ
k
n e φ p =p ∇ logp (18)
e ∇k −∇k e e T −∇k e
(cid:18) e (cid:19)
assuming small disturbances from the equilibrium, and a field line geometry in which the parallel gradient for finite
sized disturbances is not allowed to vanish (see below), a quasistatic balance between these two forces yields the
following relationship between the disturbances:
eφ p
e
= (19)
T p
e e
e e
If the parallel responses also equalises the temperature along the field lines, we then have the combination
eφ n T
e e
= =0 (20)
T n T
e e e
e e e
This situation is called adiabatic electrons, and the response to the force in Eq. (18) is called the adiabatic response.
The importance of the pressure force is what departs gradient driven turbulence in general from the world of MHD.
The adiabatic response indicates the useful units for the electrostatic potential and all the thermodynamic state
variables, which will be scaled in terms of T /e, n , or T following the forms in Eq. (20).
e e e
The otherusefulexampleis closertothe idea ofthe gyrofluidformulation: the relationshipbetweenioninertiaand
polarisation. In the low-k limit, the ExB energy is the same as the fluid one,
⊥
v2 n M c2 2
U =n M E = i i φ (21)
E i i 2 2 B2 ∇⊥
(cid:12) (cid:12)
given in Eqs. (6,13). The functional derivative of UE with respect to(cid:12)(cid:12) φ leea(cid:12)(cid:12)ds to the polarisation density
n M c2
Ω= i i 2φ e (22)
B2 ∇⊥
given in Eq. (11). Expressing φ in terms of T /e and the densitiesein terms of n e, we find
e e
e Ω ni 2 2 eφ
= ρ (23)
n e n s∇⊥T
e e e
e
where ρ is the drift scale given by
s
T M c2
2 e i
ρ = (24)
s e2B2
This is equally well described in terms of the ion gyroradius,
Ω n T eφ
i e 2 2
= ρ (25)
n e n T i∇⊥T
e e i e
e
The role of the drift scale remains, however,even if T 0, because it is a purely inertial phenomenon. With several
i
→
species present, T and ρ make good choices for the basis of the normalisationscheme, and we will use them herein.
e s
We workin terms ofanarbitrarynumber ofchargedfluids, eachwitha particularbackgrounddensity andtemper-
ature, andmass and chargestate per particle. In terms of the electrondensity and temperature,n andT , the mass
e e
of a main ion species M , and the unit of charge, e, we have a normalised backgroundcharge density
i
n Z
z
a = (26)
z
n
e
5
temperature to charge ratio
T
z
τ = (27)
z
ZT
e
and mass to charge ratio
M
z
µ = (28)
z
ZM
i
for each species labelled by z. For electrons, these are a = τ = 1 and µ = m /M . For each species, the
e e e e i
background mass density is given by a µ , pressure by a τ , and s−quared gyrora−dius by ρ2 = µ τ , in units of
z z z z z z z
n M , p = n T , and ρ2, respectively. In fusion applications the main ion is often considered to be deuterium,
e i e e e s
and so the latter mass M is used for M . However, with discharge experiments run at fusion-relevant normalised
D i
parameters using several disparate ion types becoming increasingly important [11], it is important not to make this
choice automatic.
We work under drift ordering, also called gyrokinetic ordering. There are two basic perpendicular length scales:
the drift scale ρ and the background profile scale L . Drift ordering assumes their ratio to be small,
s ⊥
δ =ρ /L 1 (29)
s ⊥
≪
where δ is called the drift parameter. It also assumes that the relative amplitude of all the disturbances is small by
this same order, for example,
eφ n T
z z
=O(δ) (30)
T ∼ n ∼ T
e z z
e e e
and similarly for the parallel flux variables (velocities, heat fluxes) in terms of the sound speed c given by
s
2
c =T /M (31)
s e i
It also takes a “maximal” ordering with respect to perpendicular wavenumbers and the drift scale,
k ρ 1 (32)
⊥ s
∼
while assuming a flute/drift ordering for the parallel wavenumber,
k /k =O(δ) (33)
k ⊥
The significance of these statements taken together is that the dynamics can be very nonlinear, in the sense that
v ∂/∂t, even though the disturbance amplitude remains small. It therefore follows that all nonlinearities are
E
drop·p∇ed∼except for the quadratic ones represented by ExB advection (v in a fluid model, additionally all its
E
·∇
finitegyroradiusrelativesinagyrofluidmodel)and,inanelectromagneticmodel,the“magneticflutter”nonlinearities
representedbythecontributionsduetothemagneticdisturbancesintheparallelgradient . Withthefluteordering
k
∇
we assume that the disturbed and undisturbed parallel gradient pieces are of similar size.
Animportantimplicationofdriftorderingonthetreatmentofthegeometryconcernsthedivergencesofthevarious
perpendicular driftfluxes. The ExBvelocityin aninhomogeneousmagneticfieldhas afinite divergence,sothat both
v n and n v would enter a fluid model. However, under drift ordering, both the profile and disturbance are
E E
adv·e∇ctedatthe∇sa·meorderbyv ,whilethefactorofnmultiplyingthedivergenceistreatedasaconstantparameter.
E
This forces an explicit split between the advection and divergence terms. Due to their various cancellations [12, 13],
the diamagnetic fluxes enter only in the divergences, so this splitting concerns solely the ExB velocity in a fluid
model, plus the associated finite gyroradius forms in a gyrofluid model. Under drift ordering, then, the advection
termappearsasapure Poissonbracketform(betweenthetwoperpendicularcoordinatesinafieldalignedtreatment,
or between all three pairs in general), multiplied by a constant coefficient, and the divergence term appears as a
curvature term (through , as above), whose properties are that (1) it is a pure divergence, (2) it is a first order
K
differential operator with the linearity property, and (3) all curvature terms are linear in the dependent variables.
The further impact of drift ordering on the treatment of the magnetic geometry is summarised in the statement
that if the three coordinates are aligned to the magnetic field such that one of them is parallel (s), one is radial (x,
acrossflux surfaces,downthe gradient),andthe third(y)hasvanishingprojectionstoboththe equilibriummagnetic
field and the background gradient, then the variation of the geometry is retained only in s. Field aligning means
that only one contravariantcomponent of B is nonvanishing, in this case Bs. These are the basic statements of field
6
aligned coordinates and magnetic flux tube geometry, as explained elsewhere [14, 15]. The principal consequence is
that the perpendicular Laplacian 2 involves only x and y, and the metric coefficients depend only upon s. Note
∇⊥
especially that this includes the variationof the field strengthB, whose variationalong the magnetic field is retained
but commutes with 2 (hence, B2 appears as a coefficient in the normalised gyroradius in all gyroaveraging and
∇⊥
gyroscreening operations). Although the equations of the standard gyrofluid model and of GEM are expressed in
covariant terms, the flux tube geometry enters when they are discretised in a numerical representation, and their
essential coordinate dependence is reflected in the abovementioned split between advection by and divergence of the
drift velocity. We will save further comment on this for the sections (below) describing the numerical scheme.
Thermodynamic state variables are normalised in terms of their background quantities, the electrostatic potential
in terms of Te/e, and the parallel magnetic potential in terms of B0ρs. Additionally, a factor of δ is folded into all
the normalisations, leaving the ExB advective nonlinearities coefficientless. This is expressed as
eφ A
φ δ−1 A δ−1 k (34)
← Te k ← B0ρsβe
e e
for the potentials, and e e
n T
n δ−1 z T δ−1 z (35)
z z
← n ← T
z z
e e
e e
u q
u δ−1 zk q δ−1 zk (36)
zk ← c zk ← n T c
s z z s
e e
for the state and parallel flux variablees respectively, where e
4πp
e
β = (37)
e B2
istheelectrondynamicalbeta(thisentersthroughAmpere’slaw). Thedensityandtemperaturearethereforetreated
on an equal footing with respect to their flux variables, the velocity and heat flux, respectively. Additionally, in
bothGEMandthestandardgyrofluidmodel,parallelandperpendiculartemperaturesandparallel-parallelandperp-
parallelheat fluxes aretreated separately,leavinga six momentmodel, for both ions andelectrons as previously[16],
and for all species herein.
We leave normalised in terms of L . Hence, the size of the contravariant magnetic unit vector component bs
k ⊥
∇
is comparable to L /qR, where q is the magnetic pitch parameter (toroidal/poloidal magnetic field, contravariant
⊥
componentratio),Risthetoroidalmajorradius,and2πqR givestheconnectionlengthalongthemagneticfieldlines.
The parameters governing core or edge turbulence result from competition between ExB turbulence and parallel
dynamics, and so the scale ratios enter the ion inertia and curvature parameters according to
2
qR 2L
⊥
ǫˆ= ω = (38)
B
L R
(cid:18) ⊥(cid:19)
respectively. These lead to the parameters governing the parallel electron dynamics,
m 0.51ν
βˆ=β ǫˆ µˆ= e ǫˆ C = e µˆ (39)
e
M c /L
i s ⊥
These are the drift Alfv´en parameter, the electron inertia parameter, and the drift wave collisionality parameter,
respectively. For core turbulence C and µˆ are small, and βˆ can be in the range from small values to several times
0.1 in high performance tokamaks. For edge turbulence C and µˆ are larger than unity, and the turbulence becomes
electromagnetic for βˆ near or larger than unity. Ideal and resistive ballooning regimes occur when βˆω > 1 or
B
Cω > 1, respectively. The magnetic shear parameter is sˆ, nominally given by dlogq/dlogr where r is the minor
B
radius but generalisable to arbitrary geometry (cf. Ref. [15]). It is usually of order unity. Under the shifted metric
fluxtube geometry, sˆenters only through the shifts incurredin the y-coordinatewhile taking derivatives with respect
to s [20]. Further details on the significance of these parameters may be found in Ref. [9].
III. ENERGETIC PROBLEMS WITH THE STANDARD GYROFLUID MODEL
Most of the problems with the standard gyrofluid model [3, 4] arise from the finite Larmor radius (FLR) terms.
Theclosureapproximationsaredonetermbyterminareasonableway,butwithsubtle differencesinthe polarisation
7
equation and in the gyrofluid moment variable equations. In the polarisation equation the gyroaveraging of the
distribution function is computed from density andtemperature moments of a perturbed Maxwellianapproximation.
In the moment equations the gyroaveragingis done onthe potential (concentrating here upon φ), which is subject to
derivatives and then to the moment integrals. These two procedures are in general different unless one approaches
them simultaneously with energy conservationin mind. For this reason,the result that the modeel has inconsistencies
isnotsounreasonable,particularlyconsideringthatitsmotivationstartswithlineartheoryandonlysecondarilyadds
the nonlinearities one needs for turbulence. Additional to this is an inconsistency in treating the higher moments
which arise from the curvature terms: a perturbed Maxwellian model is used for those terms while the model itself
retains the parallel heat flux moments as dynamical variables. These difficulties are treated in turn, with the FLR
terms first. The standard gyrofluid model is sufficiently well constructed that these repairs are in the end a minor
matter.
There are three gyroaveragingoperators used in the model, given by Eqs. (20,26,27)of Ref. [4],
φ =Γ1/2φ 1 2φ =b∂Γ10/2φ 2φ =b ∂2 bΓ1/2 φ (40)
G 0 2∇⊥ G ∂b ∇⊥ G ∂b2 0
(cid:16) (cid:17)
assuming a field aligned coeordinate seystem abndea Fourier reperesenbbtatieon in the two perpeendicular coordinates, such
that b = k2ρ2 (the model concentrates on a single ion species and leaves the electrons to be adiabatic). We relabel
⊥ i
these in terms of Γ1, Γ2, and Γ3, given by
Γ1 =Γ10/2 Γ2 =b∂∂Γb1 Γ3 = 2b∂∂b22 (bΓ1) (41)
and recast the gyroaveragedand FLR corrected potentials as
φG =Γ1φ ΩG =Γ2φ ΩG =Γ3φ (42)
for clarity. Note the factor of two iensertedeinto the deefinition feor Γ3, so tebhat ΩG aend ΩG both have the same low-k⊥
limit.
The polarisation equation [17] is given by Eqs. (7,48), rewritten as Eq. (93)e, all fromeb Ref. [4],
ni (b/2)Ti⊥ Γ0 1
n = + − φ (43)
e 1+b/2 − (1+b/2)2 τ
i
e e
preserving the normalisation of φ in terems of T , where T is the perpendiceular temperature disturbance defined as
e i⊥
the normalised(v2 1) moment of the perturbed distribution function. Hence three gyroaveragingoperatorsappear
⊥−
in the moment equations, but onely two appear in the poelarisation. Moreover, the Pad´e approximate forms are used
in the latter but not in the former. Even if this were repaired, using
Γ1/2 1 (44)
0
→ 1+b/2
in the gyroaveragingof the potential (cf. Sec.III.C.4 of Ref. [3]), it wouldbe impossible to reconcile the fact that the
2 operator does not appear in the polarisation equation, which in terms of the gyroaveragingoperators reads
∇⊥
bb ne =Γ1ni+Γ2Ti⊥+ Γ0−1φ (45)
τ
i
ForthepurposesofenergeticconsistencyitdeoesnotematterehoworwhetheertheΓ0 1screeningtermisapproximated,
−
only which operators appear in the n and T terms.
i i
The ExB energy is given by
e e
1 Γ0φ2
dΛ − (46)
τ 2
Z i
e
Using Eq. (45) we may replace this according to
1 Γ0φ2 φne Γ1ni+Γ2Ti⊥
dΛ − + = dΛφ (47)
τ 2 2 2
Z i ! Z !
e ee e e
e
8
Using the Hermitian property of the gyroaveraging operators and the definitions of the potentials, and placing the
electron contribution on the right side, we find
1 Γ0φ2 φGni+ΩGTi⊥ φne
dΛ − = dΛ (48)
Z τi 2 Z 2 − 2 !
e e e e e ee
The drift energy is therefore expressed as a combination of potentials multiplied by thermal state variables.
If the electrons are adiabatic due to fast parallel dynamics on closed flux surfaces, the electron density is itself
replaced by the potential according to
n =φ φ (49)
e
−
D E
where the angle brackets denote the flux surface (“zeonal”e) avereage. In this case the electron contribution is a sort of
field energy which is combined with the proper drift energy to form the potential energy given by
φ φ
− 1 Γ0φ φGni+ΩGTi⊥
dΛφ + − = dΛ (50)
Z e 2DeE τi 2 Z 2 !
e e e e e
e
The flux surface average is subtracted because the adiabatic state arises through the large value of the parallel
wavenumberk combinedwith the electronthermalvelocity,in comparisonwiththe dynamicalfrequencies,while for
k
the zonalcomponent there is no action by . The adiabatic electronapproximationdoes not hold in general,but it
k
∇
is assumed in Refs. [3, 4], and keeping it within this Section makes this discussion more transparent.
The thermal free energy is given by the fluid moment variables in quadratic combination. The state variable free
energy is given by
n2+T 2 +T 2
i i⊥ ik
dΛ τ (51)
i
2
!
Z e e e
while the flux variable free energy (the “generalised parallel kinetic energy”) is given by
2 2 2 2
u +q + q
k i⊥ 3 ik
dΛ (52)
2
!
Z
e e e
Comparisonofallthese pieces leadsto the conclusionthatthe combinationsτ n +φ andτ T +Ω shouldappear
i i G i i⊥ G
together under first derivatives in linear terms in the moment equations in order that in combination among all the
energy pieces the various terms reduce to terms involving single first derivativee opeerators acteing upeon combinations
of the variables,which have the form of total divergences,either B or . It is clear that since no operationby the
k
∇ K
third gyroaverage Γ3 appears in the polarisation equation with only densities and temperatures kept in the closure
for the total space density, there is no place for the third gyroreduced potential Ω in the moment equations. This
G
is the first and most obvious inconsistency of the standard gyrofluid model, and it impacts the parallel dynamics,
the curvature terms (quasistatic compressible part of the drift dynamics), and thebe magnetic pumping process, each
of which we presently examine in turn. Finally, the additional inconsistency in the treatment of higher moments in
the curvature terms in the equations for q and q is addressed.
ik i⊥
Wefirstlookattheparalleldynamics. Thisinvolvesconservativeenergytransferinthesoundwavesandconductive
heatfluxes. Thesimplestcaseofasoundwaveintheabsenceofanyeffectsduetothepotentialisofaparallelgradient
e e
in the pressure causing a parallel flow, and the corresponding parallel compression of that flow acting to restore the
pressure disturbance. The total divergence term expressing energy conservation in such a local model as this one is
B (pu /B), up to numerical constant factors dependent on how the temperatures are described. This divergence
k k
∇
represents a transport process, in this case advection of thermal energy by u . When the potential is also involved,
k
charge currents become part of this, but the structure is the same.
ee
In the standard gyrofluid model, the part of the dynamics involving sound waves is
e
∂n u
i k
= B (53)
∂t − ∇kB
∂uk e e
= τ n +T +φ (54)
∂t −∇k i i ik G
h (cid:16) (cid:17) i
e 1∂eTik =e B euk (55)
2 ∂t − ∇kB
e e
9
in additionto the polarisationequation. Evolutionof the potential energy,the thermalstate variable energy,andthe
thermal flux variable energy pieces involved in this is given by
∂ni φGuk
φ = B +u φ (56)
G ∂t − ∇k B k∇k G
uk∂∂utek =e−uk∇k τi nei+eTike+φGe (57)
eτinei∂∂nti =e−B∇hkτ(cid:16)ineBiuk +e u(cid:17)k∇keτinii (58)
1τ T ∂eTike= B τiTikeuek +ue τ Te (59)
2 i ik ∂t − ∇k B k∇k i ik
e e e
Whenthe piecesaresummed,the transfeerterms denotedbyuk k caencel,leeavingthe totaltransportdivergenceterm
∇
denoted by
e
τ n +T +φ u
i i ik G k
B (60)
∇kh (cid:16) B(cid:17) i
e e e e
that is, just the divergence of a transport flux given by the force potential in the equation for u multiplied by
k
u B/B. Here we note that T is not involved, so there is no action by Ω . We find by this analysis that the sound
k i⊥ G
wave dynamics is in order and does not require modification. The same conclusion results from exeamination of the
neon-closure part of the heateconduction dynamics, for which the equatione parts are
1∂T q
ik ik
= B (61)
2 ∂t − ∇k B
∂qik e 3 e
= τ T (62)
∂t −2∇k i ik
h i
e
and e
∂T q
i⊥ i⊥
= B (63)
∂t − ∇k B
∂qi⊥ e e
= τ T +Ω (64)
∂t −∇k i i⊥ G
h i
e
and the energetics parts are e e
1 ∂Tik τiTikqik
τ T = B +q τ T (65)
2 i ik ∂t − ∇k B ik∇k i ik
e e 32qik∂e∂qitek =−eqik∇khτiTeiki (66)
h i
e
and e e e
τ T +Ω q
∂T i i⊥ G i⊥
i⊥
τ T = B +q τ T +Ω (67)
i i⊥ ∂t − ∇kh B i i⊥∇k i i⊥ G
e e e∂q e h i
e qi⊥ ∂it⊥ =−eqi⊥∇k τieTi⊥+eΩG (68)
h i
e
Here we note the factor of two difference in the definiteion of qik heree (conformeal wieth the Braginskii definition [18])
and in the standard model.
There are also magnetic pumping terms in the gyrofluid parallel dynamics, due to the combination of parallel flow
e
and conduction, magnetic moment conservation at the gyrokinetic level, and the parallel gradient in the strength of
the magnetic field. Here, the standard model has the terms in the right places except for a single occurrence of the
“forbidden” potential Ω :
G
eb ∂uk = τ T T +Ω logB (69)
∂t − i i⊥− ik G ∇k
h (cid:16) (cid:17) i
e
e e e
10
1∂Tik
= q +u logB (70)
2 ∂t − i⊥ k ∇k
∂eT (cid:0) (cid:1)
i⊥ = qe +ue logB (71)
∂t i⊥ k ∇k
∂q e (cid:0) (cid:1)
i⊥
∂t =− τi Ti⊥−Tik +2ΩeG−ΩeG ∇klogB (72)
(cid:20) (cid:16) (cid:17) (cid:21)
e e e eb e
If we merely replace Ω by Ω we restore consistency,
G G
eb e ∂qi⊥ = τ T T +Ω logB (73)
∂t − i i⊥− ik G ∇k
h (cid:16) (cid:17) i
e
Whenthe freeenergypiecesareconstructed,theseeffeectsfoermtranesferchannelswhichthenproperlyconserveenergy
becauseτ T andΩ occurincombination. ThereasonΩ is“forbidden”isthatitdoesn’tappearinthepolarisation
i i⊥ G G
equation. If it is present, then the conservation in the exchange between potential and thermal energy is broken.
Asimilaerprobleme appearsinthe curvatureterms. Forebthe thermalstate variablesinthe standardmodel,these are
∂ni ΩG pik+pi⊥
= φ + +τ (74)
G i
∂t K 2 2 !
e e e e
1∂Tik e φG+τipik
= +τ T (75)
2 ∂t K 2 i ik!
e e e
e
∂T φ +Ω +τ p 3τ T +Ω +2Ω
i⊥ G G i i⊥ i i⊥ G G
= + (76)
∂t K 2 2
e e e e e e eb
Again, we need merely replace Ω by Ω to restore consistency,
G G
eb ∂Ti⊥e φG+ΩG+τipi⊥ τiTi⊥+ΩG
= +3 (77)
∂t K 2 2 !
e e e e e e
so that τ T and Ω again occur in combination. Once more, thermal free energy was already conserved in the
i i⊥ G
standard model, but due to Ω a mismatch in the FLR part of the transfer between potential energy and thermal
G
free energyeremaineed.
The curvature terms in theebfluid moment flux variables present a different problem, the only inconsistency in the
standard model which is not a FLR effect. There are no curvature terms involving the potential in the equations for
u , q , and q , but there is a closure treatment at the level of the fifth moments which appears in the equations for
k ik i⊥
the third moments (q and q ), arising from the factors of v2 and v2 in the grad-B and curvature drift terms in
ik i⊥ ⊥ k
teheegyrokineetic equation. In the standard model the closure for the 4th and 5th moments is taken from a perturbed
Maxwellian (cf. its Eeqs. 81 aned 82). However, the model retains qik and qi⊥ as dynamical variables, and so the 5th
moment should include contributions from pressures times conductive heat fluxes. For example, parts of the v2v3
⊥ k
moment is provided under drift ordering by p q and p q , aned part oef the v4v moment is provided by p q .
i⊥ ik ik i⊥ ⊥ k i⊥ i⊥
In the normalisation, the factors of p and p are replaced by unity. A way to do this systematically is to express
ik i⊥
the perturbed distribution function as a general six-term polynomial in which each of the coefficients is represented
e e e
by one of the fluid moment variables retained in the six-moment model. Then, the fifth moments are computed by
evaluating the integrals over v2v or v3 times the perturbed distribution function. The result of this calculation is a
⊥ k k
combination which automatically conserves free energy within the fluid moment system:
∂u τ
k i
= 4u +2q +q (78)
∂t 2K k ik i⊥
e ∂qik (cid:0) τi (cid:1)
∂t = 2eK 3uek+8eqik (79)
∂eqi⊥ τi (cid:0) (cid:1)
∂t = 2K uek+6qei⊥ (80)
e (cid:0) (cid:1)
e e
