Critical evaluation of the neoclassical model for the equilibrium electrostatic field in a tokamak

0.34 MB·English

by  Robert W. Johnson

9
0
0 Critical evaluation of the neoclassical model for
2
v the equilibrium electrostatic field in a tokamak
o
N

4 Robert W. Johnson
2
Alphawave Research

] Atlanta, GA, USA
h
p
[email protected]
-
m
s November 24, 2009
a
l PACS: 28.52.-s, 52.30.Ex, 52.55.Fa
p
.
s
c
i Abstract
s
y
h The neoclassical prescription to use an equation of motion to de-
p terminetheelectrostatic fieldwithinatokamak plasmaisfraughtwith
[
difficulties. Herein, we examine two popular expressions for the equi-

3 librium electrostatic field so determined and show that one fails to
v
withstand a formal scrutiny thereof and the other makes a prediction
0
8 for the electron temperature profile which does not compare well to
7
that commonly seen in a tokamak discharge. Reconsideration of the
0
justification for the presence of the equilibrium electrostatic field in-
.
6
dicates that no field is needed for a neutral plasma when considering
0
8 the net bound current defined as the curl of the magnetization. With
0 any shift in the toroidal magnetic flux distribution, a dynamic electric
:
v field is generated with both radial and poloidal components.
i
X
r
a 1 Introduction
The neoclassical prescription to use an equation of motion to determine the
equilibrium electrostatic field within a tokamak plasma, rather than Gauss’s
law or Poisson’s equation [1, 2, 3], is fraught with difficulties alleviated by
insisting on the neutral fluid limit rather than the quasineutral approxima-
tion [4, 5]. Herein, we examine two popular expressions for the electrostatic
1
field so determined, the first established from electron momentum conserva-
tion and the second from the Pfirsch-Schlu¨ter current, and show that one
fails to withstand a formal scrutiny thereof and the other makes a predic-
tion for the electron temperature profile which does not compare well to that
commonly seen in a tokamak discharge. Reconsideration of the justification
for the presence of the equilibrium electrostatic field indicates that no field is
needed when all bound currents are considered and that the field must either
vanish or result from sources for a model based upon the physics of a neutral,
conducting fluid. Shifting of the toroidal magnetic flux density results in a
circulating electric field with both radial and poloidal components.
Atissueisthevalidityofthequasineutralapproximation, whichallowsfor
a divergenceful electric field in the absence of a non-vanishing space charge
density ρ n e , formally expressed as ∇ E = 0 for ρ = 0, and
e ≡ s s s · 6 e
requires the determination of the electrostatic potential from an equation
P
of motion. However, the quasineutral approximation does not respect the
mathematics of electrodynamic field theory,
∇ E = ρ /ǫ . (1)
e 0
·
Fromaparticlephysicist’s field-theoreticperspective [6,7,8,9,10], thegauge
invariant Maxwell field tensor Fµν ∂µAν ∂νAµ is known to carry only 3
≡ −
scalar degrees of freedom in media, not 3 for each of the electric and mag-
netic fields, embodied by the four-potential Aµ (Φ/c,A) subject to the
≡
gauge condition and coupled to sources given by the conserved four-current
Jµ (cρ ,J) through the inhomogeneous Maxwell equations ∂ Fµν = µ Jµ,
e ν 0
≡
and the homogeneous equations are recognized as the Bianchi identity for
electromagnetism given by the field equation for the dual tensor ∂ Fµν = 0
ν
and are satisfied identically when written in terms of the electromagnetic po-
tential hence do not determine any degrees of freedom, thus the electerostatic
field is determined by the space charge density ρ and not by an equation
e
of motion. The equations by Maxwell may be expressed succinctly using
intrinsic, geometric notation as d∗dA = J in terms of the exterior derivative
d, the Hodge dual ∗, the connection 1-form A, and the current 3-form J,
as given in many standard quantum field theory texts, such as Ryder [6],
or more esoteric monographs, such as Davis [7]. What this shows is that
Gauss’s law may not be isolated from the remainder of the source bearing
Maxwell field equations, ∇ B µ ǫ ∂E/∂t = µ J, as they are but one
0 0 0
× −
unit of truth.
2
Note that we are not criticising the kinetic approach to plasma calcula-
tions, which in its original inception as the Vlasov-Maxwell system [11] fully
respects the microscopic electrodynamic field theory, but rather the neoclas-
sical (non-classical) fluid model based on the quasineutral approximation,
which does not respect the macroscopic electrodynamic field theory. While
various approximations are made in the drift kinetic equations [12, 13], most
practicalnumericalevaluationsaddressPoisson’sequationdirectly[14]. With
this manuscript, we examine in detail the mathematical difficulties one en-
counters when following the neoclassical prescription. Similar discussions on
the existence of the whistler oscilliton in geophysical plasmas [15, 16, 17, 18]
and of the helicon wave in propulsion devices [19, 20] are noted. Criticism of
the quasineutral approach in a cosmological context [21, 22] has also recently
appeared, and it is time for the fusion community to address these difficulties
head on.
2 Electrostatic field from the species equa-
tions of motion
Some geometric nomenclature and vector identities are given in the Ap-
pendix. Thefirstequilibrium∂/∂t 0electrostaticfieldunderconsideration
→
is one commonly used in the analysis of tokamak experiments [23, 24, 25, 26],
determined by integration of the electron poloidal and ion radial equations
of motion in the large aspect ratio, concentric circular flux surface approx-
imation. The use of concentric circular flux surfaces with a toroidal inte-
grating measure is pursued herein to remain consistent with the model as
presented in the literature, as are the expansions of the electron density
n = n0(r)[1+nc(r)cosθ +ns(r)sinθ] and electrostatic potential. Refer-
e e e e
ence [23] states that “the electron momentum balance can be solved for
Φc,s Φc,s/ε = nc,s/ε(eΦ0/T ), which represents the poloidal asymmetry in
≡ e e
the electrostatic potential.” Let us examine that statement in detail.
e This model [23, 24, 25] writes the equilibrium poloidal equation of motion
for arbitrary species s as
[n m (V ∇)V +∇ Π ] θˆ+∂p /r∂θ F +n e (V B E ) = 0 ,
s s s s s s sθ s s sr φ θ
· · · − −
(2)
where p = n T for T k T and F is the friction term, and takes the
s s s s B s s
←
poloidalcomponent oftheelectrostatic field ona fluxsurface atr inCoulomb
3
gauge as
∂Φ(r,θ) ∂
E = Φ0(r)[1+Φc(r)cosθ+Φs(r)sinθ] (3)
θ
≡ − r∂θ −r∂θ
Φ0(r)
= [Φs(r)cosθ Φc(r)sinθ] , (4)
− r −
whereΦistheelectrostaticpotential,indicatinganexpansionaroundΦ0(r)
rdrE0 = 0 for a last closed flux surface at r = a, where the radi≡al
− a r 6
electrostatic field is calculated from an ion equation of motion [26]. The
R
resulting evaluation of the flux surface unity, cosine, and sine moments of the
electron poloidal equation of motion with ∂T /∂θ = 0 (where other terms are
e
assumed negligible at equilibrium),
T ∂n /∂θ = en rE , (5)
e e e θ
−
defined by the expressions A dθ 1,cosθ,sinθ (1+εcosθ)A/2π,
{U,C,S}
h i ≡ { }
yields three equations which have only trivial solution. Specifically, we have
H
the system of equations
U : εnsT = eΦ0(εΦs +ncΦs nsΦc) , (6)
e e e − e
C : nsT = eΦ0(4Φs +3εncΦs εnsΦc)/4 , (7)
e e e − e
S : ncT = eΦ0(4Φc +εncΦc εnsΦs)/4 , (8)
e e e − e
valid ε,nc,ns,andT . Thetermswithfactorsofε r/R abovearestrictly
∀ e e e ≡ 0
due to the toroidal geometry and would disappear for a cylindrical plasma
column R ; one cannot address the extension to toroidal geometry
0
→ ∞
of other aspects of the model [27] without addressing the extension here.
Solution in pairs given finite (fixed) Φ0 yields inconsistent values of Φc,s and
an overdetermined system, which therefor has no solution, thus the poloidal
electrostatic field in this neoclassical model, which fails to consider the O(ε)
terms within the cosine and sine moment equations, is unphysical. Failing to
include the O(ε) terms indicates expressions applicable only on the magnetic
axis r = 0, where ε = 0 for R = , yet this neoclassical model is commonly
0
6 ∞
used to address the physics near the edge of the confinement region [23].
This system may be put into linear, homogeneous form Ax = 0 by dividing
through by Φ0,
εnsT /e ns ε+nc 1/Φ0 0
− e e − e e
4nsT /e εns 4+3εnc Φc = 0 , (9)
 − e e − e e    
4ncT /e 4+εnc εns Φs 0
− e e e − e
    
4
thus its only exact solution for ε = 0 is trivial (1/Φ0,Φc,Φs) (0,0,0),
6 ≡
which we interpret to mean exactly what it says, that it is solved when
Φ0(r) = , displaying its unphysical definition when determined from
±∞
the electron poloidal equation of motion. The matrix has rank 3, thus the
solution has 0 free parameters [28]. For a cylindrical column ε 0, one
→
recovers a matrix of rank 2 and the solution Φc,s = nc,s(T /eΦ0) with Φ0 a
e e
free parameter. If one assumes that the smallness of the coefficients may
be represented by ε, eg nc,s εnc,s, then the offending terms near the edge
e ≡ e
where ε approaches 1/2 represent up to a 19% correction to the leading order
equations. As non-vanishing Φc,s are an integral part of the development of
e
this neoclassical model and appear in its remaining equations without a pref-
actor of Φ0 through substitution, the validity of its conclusions is in jeopardy.
Note that a putative non-vanishing radial electrostatic field without poloidal
variation demands the existence of no poloidal electrostatic field, else the
poloidal variation to the potential ruins the poloidal symmetry of the radial
field; if that radial field is determined from a radial equation of motion then
the associated poloidal field is determined by the poloidal dependence of that
equation.
One might think to alleviate the difficulty by invoking the logarithmic
derivative, writing the electron poloidal equation of motion as ∂(lnn
e
−
eΦ/T )/∂θ = 0, with solution C(r) = n (r,θ)exp[ eΦ(r,θ)/T (r)]. Expand-
e e e
−
ing the exponential gives, for < eΦ/T < ,
e
−∞ ∞
∞ ( eΦ/T )k eΦ 1 eΦ 2 1 eΦ 3
e
C(r) = n − n 1 + +... ,
e e
k! ≈ − T 2 T − 6 T
k=0 " e (cid:18) e(cid:19) (cid:18) e(cid:19) #
X
(10)
and to each order in eΦ/T , taking the flux surface moments yields three
e
equations to solve. Using the expansions for n and Φ above, the zeroth or-
e
der equations give C(r) = n0 and nc,s = 0, and the first order equations yield
e e
ε-dependent solutions which, upon restriction to the magnetic axis ε = 0, are
C(r) n0[1 (nc)2/2 (ns)2/2](1 eΦ0/T ) and Φc,s nc,s(T /eΦ0 1).
→ e − e − e − e → e e −
Agreement is found with the previous Φc,s = nc,s(T /eΦ0) provided eΦ0/T
e e e ≪
1. Considering now a ratio eΦ0/T 10 1, such as may be found in the
e
∼ ≫
plasma core [24, 29, 30], the expansion in Equation (10) requires 10 orders
before the absolute magnitude of the terms starts decreasing, indicating that
a low-order approximation is not an exact solution. Considering the second
order equations leads to solutions for C(r) and Φc,s only when ε = 0 which
diverge when nc,s = 0, and no solution has yet been found to the nonlin-
e
5
ear equations for general ε. The difficulty with Φ0 encountered above has
simply been shifted to the “undetermined” function C(r), which is perfectly
determinable in principle from the system of equations. We remark that,
assuming the model’s system to be well-determined (ie equal numbers of in-
dependent equations and degrees of freedom) before invoking the logarithmic
derivative, the introduction of C(r) without an additional equation for its de-
termination leaves the system under-determined—if C(r) is identified as an
independent degree of freedom, then it needs its own independent equation,
else if not, then one of those degrees of freedom must be determined by the
unity moment equation here. Pursuing the argument, as the electron density
n0 divides out of Equation (5), its poloidal variations nc,s are determined by
e e
continuity ∂n /∂t + ∇ n V n˙ = 0 for particle source rate n˙, and
e e e C,S
h · − i
the thermal energy T is determined by the heat equation and given as in-
e
put from experimental measurement for the analysis, the remaining degree
of freedom Φ0 must be determined by C(r) as given by the poloidal force
balance system, Equation (9).
Notwithstanding the above, let us now consider an alternate expansion of
the exponential which extracts the θ dependence using ξ Φccosθ+Φssinθ,
≡
so that
C(r) = n exp( eΦ0/T )(expξ)−eΦ0/Te , (11)
e e
−
0 2 −eΦ0/Te
n exp( eΦ /T ) 1+ξ +ξ /2+... , (12)
e e
≈ −
n0exp( eΦ0/T )(1+nccosθ+nssinθ) (13)
≈ e − e (cid:0) e e (cid:1) ×
eΦ0 1 eΦ0 2
2
1 ξ + ξ +... , (14)
− T 2 T
" e (cid:18) e (cid:19) #
noting the use of low-order truncations to the two infinite series (which can
never be as exact as the algebraic method above) and that technically the
expansion found in Equation (10) is more democratic in its treatment of
the powers of eΦ0/T . Taking the flux surface moments as before, we find
e
nearly the same first order solutions, whereupon reduction to the magnetic
axis (or cylindrical geometry) yields C(r) n0exp( eΦ0/T )[1 (nc)2/2
→ e − e − e −
(ns)2/2] and Φc,s nc,s(T /eΦ0), and for the second order equations the
e → e e
solution on the magnetic axis diverges and for general ε cannot be found.
Inserting factors of ε at the cost of the loss of the cylindrical interpretation
of ε 0, which still identifies the magnetic axis in toroidal geometry, yields
→
6
Equation (11) in terms of the tilded coefficients (which need not be small),
C(r) n0exp( eΦ0/T )[1+ε(nccosθ+nssinθ)] (15)
≈ e − e e e ×
eΦ0 eΦ0 2
1 ε(Φccosθ +Φssinθ)+ 1 ε2(Φccosθ +Φssinθ)2 ,(16)
− T 2 T e e
" e (cid:18) e (cid:19) #
e e e e
and affords additional control over the terms retained in the expansion. The
resulting system isoftheform0 = U +U ε2+U ε4 = C ε+C ε3 = S ε+S ε3
0 2 4 1 3 1 3
and has not solution until reduced to order O(ε2) which, as a factor of ε
cancels out of the ... moments, is the first appearance of an explicit ε
C,S
h i
dependence in the equations. Its expression for Φc contains terms with 1/ε2,
thusthephysicalΦc whenε = 0. Onlywheneachequationisreducedto
→ ∞
itsleadingtermdoesthesystemreturnthesolutioensC(r) = n0exp( eΦ0/T )
e − e
and Φc,s = nc,s(T /eΦ0). Comparing the leading coefficients in Equation (16)
e e
gives the ratios 1 : 10ε : 55ε2 for eΦ0/T 10, and allowing for eΦ0/T 1
e e
∼ ∼
yields 1 : ε : ε2 which still represents up to a 25% correction. The lesson here
is that knowledge of a series’ convergence need not imply that a low-order
approximation is an adequate representation of the function. Concern over
the expansion in no way detracts from the observation that the physics of
the situation is embodied by the algebraic Equations (6) which contain no
exponential factor thus represent a more exact method of solution.
Finally, we consider the consistency of using two different species’ equa-
tions of motion to determine the species independent electrostatic potential,
taking the radial electrostatic field from the equation of motion for arbitrary
ion species j as
1 ∂p
j
E = +V B V B , (17)
r φj θ θj φ
n e ∂r −
j j
and dropping the convective term as is standard practice in the field [24, 25,
26]. Theusualevaluationofthatexpressionfromexperimentalmeasurements
neglects any poloidal dependence; however, when the intrinsic variation of
quantities induced by the geometry is considered [31], given in the large
aspect ratio ε 1, concentric circular R R , flux surface approximation
r 0
≪ ≡
B = (0,B ,B ) by
θ φ
B = B0/(1+εcosθ) , V = V0 /(1+εcosθ) , V = V0 (1+εcosθ) , n = n0 ,
θj θj φj φj j j
(18)
where A0 is the average of the values of A on the vertical midplane, a geo-
metric dependence is introduced. Using these values (r dependence implied),
7
we find
p′0
j 0 0 0 0 2
E (r,θ) = +V B V B /(1+εcosθ) , (19)
r n0e φj θ − θj φ
j j
for species pressure gradient p′ ∂n T /∂r, thus E (θ) E0. Expanding
j ≡ j j r 6≡ r
the denominator reveals a power series in εcosθ,
p′0
j 0 0 0 0 2 2 3 3
E (r,θ) +V B V B 1 2εcosθ+3ε cos θ 4ε cos θ(20,)
r ≈ n0e φj θ − θj φ − −
j j
p′0 (cid:0) (cid:1)
j 0 0 0 0 0 0
+V B V B +2εV B cosθ (21)
≈ n0e φj θ − θj φ θj φ
j j
0 1
= E (r)+E (r)cosθ , (22)
r r
and as integration with respect to r does not affect the θ dependency, the
potential associated with the radial electrostatic field relative to its central
r
value, Φ (r,θ) drE (r,θ), may be written as a cosine series,
Er ≡ − 0 r
Φ (r,θR) Φ0 +Φ1 cosθ = Φ0 1+Φc cosθ , (23)
Er ≈ Er Er Er Er
where Φc = rdrE1/ rdrE0 = 0 in general,(cid:0)noting that t(cid:1)he potential on
Er 0 r 0 r 6
the last closed flux surface at r = a is not single valued. Returning now
R R
to the electrostatic potential appearing in the poloidal equation of motion,
expressed to leading order as
Φ (r,θ) = Φ0 (r)[1+Φc (r)cosθ +Φs (r)sinθ] , (24)
Eθ Eθ Eθ Eθ
where Φ0 (r) rdrE0 = 0 is the potential relative to that of the last
Eθ ≡ − a r 6
closed flux surface, with solution Φc,s = nc,s(T /eΦ). No loss of generality
R Eθ e e
ensues if oneredefines the potentialrelative to itscentral value. Withvanish-
ing extrinsic poloidal dependence to the electron density, nc,s 0 as above,
e →
the potential retains no explicit poloidal dependence, Φ (r,θ) Φ0 (r).
Eθ → Eθ
Thus, we conclude that the electrostatic potentials associated with the ra-
dialandpoloidalelectrostaticfieldsevaluatedfromtheionandelectronequa-
tions of motion by this neoclassical model are inconsistent, Φ = Φ , as
Er 6 Eθ
Equation (23) does not equal Equation (24) in the case of vanishing density
asymmetries nc,s = 0.
e,j
8
3 Equilibrium field from the Ohm’s law equa-
tion
3.1 Derivation
Examining the expression of another leading contender for the equilibrium
electrostatic field [3] evaluated from the Ohm’s law equation and the Pfirsch-
Schlu¨ter current, one may put its poloidal component
E B /B B2 E B 1/B B2 1
E = h φ φ θi φ φ +RB p′η h θi , (25)
θ B2/B B − B φ k B2/B B − B
h θi θ θ (cid:18)h θi θ θ(cid:19)
into the form
E = Ec 2εcosθ (ε2/2)cos2θ , (26)
θ θ −
whentheShafranovshiftisneg(cid:2)lected, asintheconcen(cid:3)triccircularfluxsurface
approximation of above, upon application of Stokes’ theorem to Faraday’s
law[32], ie by requiring dθE = 0. Notethat thisneoclassical modelforthe
θ
poloidal electrostatic field differs distinctly from that of the previous section
H
in detailed functional form. Inserting Equation (26) into Equation (5) and
taking the flux surface Fourier moments yields three equations which have a
nontrivial solution only when expanded to order O(ε3), given by
nc ε3/(6ε2 8)
e −
ns = ε√3ε4 168ε2+192/(18ε2 24) , (27)
 e   ± − − 
Ec 4(T /eR )/ε√3ε4 168ε2+192
θ e 0
± −
   
thus the presence of a poloidal electrostatic field of that form should be
accompanied by a potentially measurable shift in the electron density profile.
Notethatthederivationimmediately preceding isslightlyinconsistent, asthe
associated electrostatic potential Φ = Φ (2εsinθ (ε2/4)sin2θ) is of the
axi 0
−
correct harmonic form [33, 34, 35, 36] for axial geometry, as is easily verified
in (Z = rsinθ,R = R +rcosθ,z) coordinates via application of the axial
0
Laplacia−n ∇2 ∂2/∂Z2 + ∂2/∂R2 to Φ = Φ Z[(R R )/2R2 2/R ],
axi ≡ axi 0 − 0 0 − 0
yet the flux surface average is done in toroidal geometry. Note that this Φ
0
is not the Φ0 of the preceding section but is a unit bearing constant which
sets the scale.
In order to achieve the correct harmonic form for tokamak geometry,
the potential must satisfy the toroidal Laplacian, which in cylindrical co-
2 2
ordinates [32, 37, 38] is given by ∇ ∇ + ∂/R∂R (note that the
tor ≡ axi
9
expression for ∆⋆ given by Hopcraft [2] is not the toroidal Laplacian and
differs by the sign of the additional term), and the form of the geometric
term hints at the solution. Direct integration yields the toroidal potential
Φ Φ (R lnR), from which E = Φ [(lnR R )/2R2 2/R ] and
tor axi Z 0 0 0 0
≡ → − − −
E = Φ Z/2RR2, noting that the introduction of the logarithm breaks the
R 0 0
−
usually obvious relation between the symbol for the magnitude of a quantity
and the units associated with that quantity—carefully pulling the units be-
side the leading coefficients of expressions ensures that they are respected.
From these, we determine the poloidal field to be
1∂Φ Φ ∂Z ∂ ∂R ∂ Φ E E
E = − 0 + = Ec (R R ) Z Z R ,
θ ≡ −r ∂θ r ∂θ ∂Z ∂θ ∂R Φ θ − 0 Φ − Φ
0 0 0
(cid:18) (cid:19)(cid:18) (cid:19) (cid:20) (cid:21)
(28)
where we identify Ec Φ /r, and the corresponding radial field is
r,θ ≡ − 0
∂Φ Φ ∂Z ∂ ∂R ∂ Φ E E
E = − 0 r + = Ec Z Z (R R ) R .
r ≡ − ∂r r ∂r ∂Z ∂r ∂R Φ r − Φ − − 0 Φ
0 0 0
(cid:18) (cid:19) (cid:18) (cid:19) (cid:20) (cid:21)
(29)
In(r,θ,φ)coordinates, wehave Φ = Φ rsinθ[ ln(R +rcosθ)+5R ]/2R2,
tor 0 0 0 0
−
and
Φ r2sin2θ
0
E (r,θ) = rcosθ[ln(R +rcosθ) 5R ] , (30)
θ 2rR2 0 − 0 − R +rcosθ
0 (cid:26) 0 (cid:27)
Φ r2cosθsinθ
0
E (r,θ) = rsinθ[ln(R +rcosθ) 5R ]+ . (31)
r 2rR2 0 − 0 R +rcosθ
0 (cid:26) 0 (cid:27)
As this is an electrostatic field within a neutral medium (ie one for which
the net charge on a differential volume element vanishes) at equilibrium,
Maxwell’s equations ∇ E = 0 and ∇ E = 0 are satisfied.
× ·
As electrostatic fields necessarily require a supporting charge density [6,
32, 39, 8, 9], one may very well ask where these charges are. As the solution
to Laplace’s equation is uniquely determined by the boundary condition, we
will find them on the boundary of the region under consideration, which in
this case is the R/R weighted circle representing our outermost flux sur-
0
face at normalized minor radius r/a = 1 upon collapse of the toroidal di-
mension, giving us an inverse Dirichlet problem, which is usually defined as
solving for the potential given the boundary condition. The source of these
charges is intended to be the accumulation resulting from the divergence of
the pressure gradient driven diamagnetic current in toroidal geometry, and
10