Computer Physics Communications 40 (1986) 153—157 North-Holland, Amsterdam
153
FOKKER-PLANCK CODE FOR THE QUASI-LINEAR ABSORPTION OF ELECTRON CYCLOTRON WAVES IN A TOKAMAK PLASMA RI. MEYER, G. GIRUZZI and V. KRIVENSKI Laboratoire de Physique des Mi/leux lonises (C.N.R.S. 1/A. 835), Universith Nancy I — B.P. 239 — 54506 Vandoeuvre les Nancy Cedex, France and Association EURA TOM C,E.A. sur Ia Fusion, Centre d’Etudes Nuc/eaires — B.P. No. 6 — 92260 Fontenay-aux-Roses, France
We present the solution of the kinetic equation describing the quasi-linear evolution of the electron momentum distribution function under the influence of the electron cyclotron wave absorption, Coulomb collisions and the dc electric field in a tokamak plasma. The solution of the quasi-linear equation is obtained numerically using a two-dimensional initial value code following an ADJ scheme. Most emphasis is given to the full non-linear and self-consistent problem, namely, the wave amplitude is evaluated at any instant and any point in space according to the actual damping. This is necessary since wave damping is a very sensitive function of the slope of the local momentum distribution function because the resonance condition relates the electron momentum to the location of wave energy deposition.
I. Introduction
2. Basic equations
The theory of electron cyclotron resonance heating in tokamak plasmas is now well developed. The linear computations of wave damping are generally performed using a Maxwellian isotropic distribution [1—3]and the quasilinear treat-
The combined effects of the electron cyclotron wave absorption and the dc electric field are described by the equation
ment is mostly formulated and investigated ignoring the evolution of the wave energy in space [4—6].In fact, the wave amplitude varies in space and this results in an inhornogeneous deformation of the electron distribution. In this paper we examine the case of the combined effects of the electron cyclotron wave absorption and the dc electric field parallel to the tokamak magnetic field B0 on the momentum distribution for propagation near normal to B0. The paper is organized as follows. In section 2 we formulate the physical problem. In section 3 we present the numerical method for the solution of the Fokker—Planck equation. In section 4 we cornment some numerical results and the conclusions are given in section 5.
31 0T
—~-
/ =
a
I\ / 31 \ / 31 —i+ \B’T/cy —i-+ \3T/~ —i\3T1E
J
)
)
,
(1)
0~~
where f is the electron momentum distribution function averaged over the magnetic (assumed cylindrical) surfaces and the Larmor rotation phase, i.e., f—f (p11, p~, r, T), p1~and p~ are the cornponents of the momentum p parallel and perpendicular to the tokamak 4n magnetic field B0, respectively, T Pet ~‘e = 2‘r~e 0Am 1/2 ~j-’— 3/2, n0 is the plasma density, e is the electron charge, m is the electron rest mass, 7~is the electron temperature and A is the Coulomb logarithm. We study the evolution of the electron tail which occurs in a relatively fast time scale r < 100. The body of the distribution function is maintamed at constant temperature 7~while the tail experiences the parallel acceleration (or retardation) due to the inductive electric field E and the perpendicular heating due to the ordinary mode
0010-4655/86/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)
-~
154
R. L. Meyer et a!.
/ Absorption
wave absorption. In this model the collision operator used is the asymptotic form of the Boltzmann operator, i.e., Z+1 3 Sifl ~~con
=
00f 2
0/10f
)
(2) U 1~2,u and where is thecoordinates ion charge,with u =p/(mT~) O are Z polar the reference axis along B 0, and equal ion and electron temperature are assumed. The influence of the dc electric field E is described by the operator U
(.~ 3f\) E
~_)(cos3fl 7~ Ofu
of x along which wave damping takes place results from the process of average over the surface 2rR
2’rr
0.
For small values of 7~one has for the ordinary mode
0 ~
U
of electron cyclotron waves
~u ~j’ 3f \ (3) is the Dreicer field. In
k,—~
—
~~ 2
8mc
(6) where kr = (co/c)(1 w~/w2)”~2and k~are the real and imaginary parts of the propagation vector —
in the x direction, respectively. Eqs. (5) and (6) exhibit two important properties of the physical problem.
where E~= 2 ‘rrn 0e this paper we consider values of EII/EC for which the parallel energy of the electron tail is below the critical value for appreciable runaway production. We consider a narrow wave-packet propagating in the tokamak equatorial plane from the low B0 side in the x direction normal to B0 [8]. For a gaussian wave-spectrum centered around N~= 0 and with an angular spread of half-width ~N11the quasilinear diffusion operator is given by
an First important of all factor the diffusion governing coefficient the evolution D which of the is momentum local wave distribution, energy represented is a linear by function the factor of the P 0 exp(—2f~’k1(r’)dr’)in eq. (5). Secondly as shown by eq. (6) the damping is a sensitive function of the slope of the local momentum distribution and must be evaluated simultaneously with it. And therefore we must solve a full non-linear and self consistent problem, namely the wave energy has to be evaluated at any instant and any point in space according to the actual damping which depends upon the momentum distribution function.
(-~) =----~--u~D-~-, 0 0 0
3. Numerical method
I
T
C)’
=
(El’
U1~ U
—
(4)
U~
where 2 \ 1/2 w w ) D=~1—--U /
i03° u
The equation we solve has the general shape 2f 02f
0 rRo(L~N,,)wAno~i
Of
3
r
/
02f
xPOexp(—2Jki(r,)dr,) 0
(y
Of
Of
+a 3—~+a4~+a5—,
w~(r)\2 2(~N)2] I /u1 wJ Xexp{—~ 2/T~, y = (1 + u 2/p.) 1/2
Ou where f=f(u, 0, x,
—
.
(5)
p. = mc is the electron plasma frequency, w is the launched wave frequency, R 0 is the tokamak major radius and P0 is the total power delivered by the wave beam. The presence of the radial coordinate r instead
FJu
~) and where the coefficients
a, are very complicated functions of the distribution function via the wave damping. They have to be evaluated numerically at each point in space and time. We have developed a two-dimensional initial value code following an alternating direction implicit scheme [9] to investigate the time evolution of the distribution function from an initial
R. L. Meyer et al.
/ Absorption
Maxwellian distribution at temperature T~(x)given by the temperature profile in the tokamak. Specifically f was computed as a function of x and T and simultaneously k1 (x, T) was evaluated according to eq. (6) with the momentum distribution obtained at the preceding time step. For sufficiently small time step, the momentum distribution is so little modified that this is a fairly good approximation. The space (u, 0) extends from 0 to ‘rr in the 0 direction and from 0 to Um~ in the u direction. Full velocity space boundary conditions are applied namely f(u = 0, 0) is independent of 0, Of/O0(u, 0=0)=Of/B0 (u, 0=~)=0, Of/Ou(u=0, O~r/2)0,
ofelectron cyclotron waves
155
tion. Owing to the presence of a ~ function in the integrand of eq. (6), the computation of k1 requires a one-dimensional integration over 0, the value of u depending upon space location. This value of u never cd~ncideswith a grid point. The needed interpolation is made using a cubic spline routine.
4. Numerical results We consider the following profiles 2/a2), 2 n0(x) = n0(0)(1 x —
B 0(x)=B0(0)/(1+x/R0),
and assuming the momentum distribution to remain Maxwellian for a sufficiently high velocity modulus we add the boundary condition: f(Umax,
0)
fMaxwei1ian(~’max)~
The velocity space is divided by a grid. Each point of the grid is labeled by two indexes i and j. The first one extending from 1 to I corresponds to the u direction and the second one j ranging from 1 to J represents the 0 direction. The various derivatives of f are approximated by finite differences at each point of the grid excepted the four boundaries. For each i (or j) this method leads to a tridiagonal system of linear equations with the index j (or i) varying from 2 to J— 1 (or 2 to I—i). The system is solved by a Gauss elimination method which has a simple form in the case of tridiagonal matrix. The values of f at the boundaries are extrapolated by using the boundary conditions. In the ADI method the time step z~Tis splitted in two parts. During the first half time step the derivatives of f are approximated implicitely in one direction (for example u) and explicitly in the other direction (0). The situation is reversed during the second half time step. (By this device both u and 0 direction are treated on the same ground). We carefully checked that the momentum distnbution remainsduring normalized unity with a veryfunction good accuracy the timetointegra-
and the following parameters a = 40 cm, R0 = 150 cm, 13
~ie(0)3X10 ~(0) = 1 keV,
~N
=
10-s
w= =
.,
II
cm
-3
w~(0),
EII/EC
10—2,
F(~)
a1/ ‘\
.11
0
/1
S\
20
\ \~100\ \ \ \~ \ “\ \ \‘~
I/f
“.
.j~
..~
‘~
j
~
.7/ / /
,~‘
io
b
~
x
MW
~
\ a
2
0
~
\
=
..~ ..~
ji
// //
b
0
;uat
~ r=O, 100 for2,rno(O)= =1.25 ~ cm and
a=40 a) P cm, R0=150 cm, E~1/E~=2x10 0 =0, b) P0 1 MW.
=
156
R. L. Meyer et a!.
/ Absorption of electron cyclotron waves
The ordinary wave is launched from the low magnetic field side in the equatorial plane with the above mentioned parameters. The cyclotron energy is deposited in a region extending from x = 0 cm to x —5 cm. This narrow range allows us to neglect the spatial dependance of EII/EC. The distribution function was computed as a function of x and T with a space step of 0.5 cm from x = 0 cm to x = 5 cm and a time step of 0.05 from T = 0 to T = 100 for which a quasi-steady state is reached. Eq. (1) was integrated by using a 50 x 50 grid in the (u, 0) space. From the resonance condition and the spatial deposition of the cyclotron energy we fixed the upper U boundary at Umax = 10. The collision operator was set up at u~011 0.5. It is worth noting that the results obtained are weakly dependent on the choice of the value UCOIJ <1. Each run requires about three hours on a CII DPS8 computer whose power is 1.6 MIPS. The results obtained are illustrated by the following figures. In fig. 1 we show F(u,,) = 2.71f u1du1f
for ‘r = 0 (the initial Maxwellian) and T = 100 without and with electron cyclotron power for x = 1.25 cm. The increase of the net current in the direction of 8~due to the cyclotron power appears clearly from fig. 1. Since the wave-packet is symetric around ~ = 0 there is no net current due to the electron cyclotron wave in the absence of dc electric field. In fig. 2 we show i
T1 (U11)
/
/
/ / //
3.75
1
\//2~75
A
~ \
~ I/ V ~
I I \~\ / I
:,
I,,,
versus u1,
k1 (cm~)
4
~
T~2.rrJ0 du1 ~~U1/2f/F(U11)
versus
\/\ .~
\
at T = 100 for several values of x. T1 (U11) is found to be a slightly asymmetric function of U,, which indicates that the perpendicular enegy is mainly determined by electron cyclotron wave absorption. Moreover T1 (u,,) depends strongly upon space since wave absorption varies in space. From fig. 3 it appears that the wave damping experiences two contrasting effects namely a quasilinear reduction for x <2 cm and a small rise due to the increased superthermal population for x> 2 cm.
T~(keV)
(‘I
=
3
.75
.
3
fr=100
:
175 2
/ \t~o \
/
I
/
‘~
/i ~~~
r\~s
Fig. 2. T~(u11) versus u11 at r = 100 for several values of r and the conditions of fig. ib).
-
2:
Fig. 3. k, versus r at r
2
=
3
~
r(c~
100 for the conditions of fig. lb).
R. L. Meyer et al. 5.
/ Absorption
Conclusion
The time evolution of the electron momentum distribution function in a tokamak under the combined effects of the tokamak dc electric field and electron cyclotron wave absorption is described by a nonlinear equation. The problem is self consistent in the sense that the damping of the wave is determined by the actual momentum distribution. In order to solve this problem we developed a two-dimensional initial value code based on an alternating direction implicit scheme. This code takes into account the plasma inhomogeneities and calculates the damping simultaneously with the momentum distribution function.
of electron cyclotron waves
157
References [1] V.V. Alikaev, Yu.D. Dnestrovskii, V.V. Parail and G.V. Pereverzev, Fia. Plasmy 3 (1977) 230, Soy. J. Plasma Phys. 3 (1977)Litvak, 127. G.V. Permitin, E.V. Suvorov and A.A. Fraj[21A.G. man, NucI. Fusion 17 (1977) 659. [3] I. Fidone, G. Granata, G. Rampom and R.L. Meyer, Phys. Fluids 21(1978) 645. [4] J. Rowlands, Sizonenko K.L. Stepanov, Zh.661. Eksp. Teor. Fiz. 50 V.L. (1966) 994, Soy. and Phys. JETP 23 (1966) [5] C.F. Kennel and F. Engelmann, Phys. Fluids 9 (1966) 2377. [6] V.V. Alikaev and V.L. Vdovin, Fiz. Plazmy 9 (1983) 928. [7] J.C. Wiley, Duk-In Choi and W. Horton, Phys. Fluids 23 (1980) 2193. [8] I. Fidone, R.L. Meyer and G. Granata, Phys. Fluids 26 (1983) 3292.and K.D. Marx. in: Methods of Computational [9] J. Killeen Physics, vol. 9, eds. B. Adler, S. Fernbach and M. Rotenberg (Academic Press, New York, 1970) p. 421.