J. Plasma Fusion Res. SERIES, Vol. 9 (2010)

Alfvén wave propagation in initial value MHD simulation
in a tokamak geometry
Y. Nishimura and C. Z. Cheng
Plasma and Space Science Center, National Cheng Kung University, Tainan 70101, TaiwanFebruary 9, 2010

(Received: 30 October 2009 / Accepted: 10 February 2010)

Alfvén wave propagation in an initial value magnetohydrodynamic (MHD) simulation is reported.
The simulation result is compared with the theoretical shear Alfvén wave dispersion relation both in acylindrical geometry and in a toroidal geometry. Generation of the toroidicity induced Alfvén eigenmode(TAE) frequency gap is observed in the presence of 1/R variation of the toroidal magnetic ﬁeld. As apreliminary test, excitation of TAE mode is simulated in the presence of energetic particles. A numerical
scheme employed in the simulation (the second order moment of the kinetic ions and electrons fromparticle simulation replace the pressure evolution equation) is discussed.
Keywords: Toroidicity induced Alfvén eigenmode, kinetic-ﬂuid model, MHD simulation, gyrokinetic
theory, particle simulation

1. Introductiontion 3 presents simulation results. We summarize thiswork in Sec. 4.

Energetic particles can play important roles inburning plasmas.
For example, Toroidicity induced Alfvén Eigenmodes (TAE)[1, 2] can be excited when energetic particles (e.g. fusion born alpha particles),[3, 4] resonate with the phase velocityof the shear Alfvén wave within the frequency gap ofthe Alfvén continuum. Another example where theenergetic particles play important roles is the stabilization of sawtooth oscillations.[5, 6, 7] Long period
sawteeth accompanying the Alfvénic instabilities arereported. In this work, employing a kinetic-ﬂuid simulation model,[8] the Alfvénic instabilities are studied
in tokamak plasmas. The kinetic-ﬂuid model incorporates all the particle dynamics through the pressure
tensor by taking the second order moment of the particle simulation while the electromagnetic ﬁeld quantities are evolved in the ﬂuid equations.[8] The kineticﬂuid model retains the ion and electron wave-particle
interaction for both the bulk and the energetic particle
components. The model allows to study the interaction of the energetic particles with thermal kinetic ionsand electrons in contrast to the conventional kineticMHD hybrid models.[9, 10, 11, 12, 13] The advantage
of initial value approach is its application for nonlinear
simulation (and for studying resonant type instabilities such as RTAE).[14, 15] Our ﬁnal goal is to examine
the saturation mechanism of the instabilities incorporating geometrical eﬀects correctly. We approach stepby step toward our long term goal. We start fromexamining global Alfvén oscillation, continuum damping, and the generation of the TAE gap in the toroidal
geometry due to the magnetic ﬁeld inhomogeneity (inthe MHD limit). We then incorporate kinetic energetic particles to excite the instabilities. In Sec. 2 thebasic model for MHD simulation is discussed. Sec-

2. Basic model for computation
A nonorthogonal, straight-ﬁeld-line coordinate
system is employed in the present calculations, inwhich ρ is the ﬂux surface label, θ is the poloidallike angle, and ζ is the toroidal angle. The magnetic
ﬁeld in a tokamak is given by
B = Beq + B̃,(1)

Beq = ∇ψeq × (q∇θ − ∇ζ) ,(2)B̃ ≡ ∇ × Ã = ∇ × (−ψ∇ζ) = ∇ζ × ∇ψ,(3)where ψeq (ρ) stands for the equilibrium poloidal magnetic ﬂux and equilibrium poloidal current. Here,�ψ (ρ, θ, ζ) =m/n ψm/n (ρ) cos (mθ + nζ) is thepoloidal ﬂux function of the perturbed ﬁeld, where mand n are the poloidal and the toroidal mode numbers,respectively. We denote the total poloidal magnetic
ﬂux as ψt (ρ, θ, ζ) ≡ ψeq (ρ) + ψ (ρ, θ, ζ).
Numerical simulation has been conducted employing a reduced MHD formulation[16] which allows simulating both cylindrical and toroidal geometries. Toroidal geometry can enter the set of equations via metric elements obtained from an equilibrium solver.[17] The initial value simulation is conducted by employing the FAR code.[18, 19] In ourstudy the reduced MHD equations[16] are solved forthe magnetic ﬂux ψ, the toroidal component of thevorticity U ζ , and the pressure P . The relevant equations are: the toroidal component of ideal Ohm’s Law,
∂ψ
= −R2 (v × B) · ∇ζ,∂tthe vorticity equation,
ρmauthor’s e-mail: nishimura@pssc.ncku.edu.tw541(4)

��dU ζ
1= R2 B · ∇J ζ + ∇R2 × ∇P · ∇ζ,(5)dtc©2010 by The Japan Society of Plasma
Science and Nuclear Fusion Research

Y. Nishimura and C.Z. Cheng, Alfvén Wave Propagation in Initial Value MHD Simulation in a Tokamak Geometry

In these equations, R stands for the major radius,v ≡ R2 ∇ζ × ∇φ is the ﬂuid velocity [φ (ρ, θ, ζ) =�(ρ) sin (mθ + nζ) is the stream function],m/n φm/n��U ≡ ∇ × R2 v is the vorticity, and J = (c/4π)∇ × B
is the current density. The mass density is given byρm (we assume plasma to be incompressible and thus
ρm is constant in this paper). Superscripts denote thecontravariant components.
In the momentum balance (the vorticity equation), the pressure is given by
P = P⊥ (I − bb) + P� bb

In the δf simulation, the weight equation[23] carries the information of the perturbed ﬁeld (while thetrajectories are unperturbed ones for the linear simulation). As an initial test we provide the simulation
with only passing particles (we set μ = 0 to avoid complications such as particles lost at the plasma boundary).
We incorporate the electron inertia term, the Hallterm, and the gyroviscous term into the Ohm’s law, as
manifested in Ref.[8]. In this paper, however, we focus
on incorporating the pressure closure relation Eqs. (7)and (8).

(6)(the oﬀ diagonal tensor element is ignored for the moment)
�(7)P⊥j = μBδfj d3 v,
P�j =�mj v�2 δfj d3 v.

3. Simulation results in the MHD limitParameters used in the calculations are as follows.
The q-proﬁle (the safety factor) is taken as parabolic
q = 1.0 + ρ2 . The equilibrium pressure proﬁle has the2. Major radius is given by R = 5 m
form of P = ψeq
and the minor radius a = 1.25 m (thus the inverse
aspect ratio � = 1/4). A total of 200 equally spaced
mesh points was used in the radial direction. Themode spectrum was selected as m/n = 1/1 and m/n =2/1. Equations (4) and (5) are dissipationless. The
resistivity and the viscosity terms[19] are set to bezero. In Eqs.(4) and (5), we provide φ with an initial
perturbation of a sine parity.[20]In Fig.1 (a), shear Alfvén wave oscillation is tested
in a cylindrical limit. We set R = 1 in Eqs.(4) and (5).
In Fig.1, single Fourier mode of m = 1 and n = 1 istaken. Here, time is normalized by ωA = vA /q0 R0
which is the Alfvén frequency at the magnetic axis
(q0 = 1.0 is the safety factor and R0 is the major radius at the magnetic axis, respectively). The wave
equation (hyperbolic equation) and thus the shearAlfvén wave dispersion relation can be obtained by
Eqs.(4) and (5) in the absence of the curvature driven
term. The Alfvén wave frequency measured at diﬀerent radial locations satisﬁes shear Alfvén wave dispersion relation ω = vA k� (ρ), where k� = ρ2 /(1 + ρ2 )which is plotted onto Fig.1(b).
Figure 2 shows the shear Alfvén wave frequency
as a function of the radial coordinate ρ which isequivalent to Fig.1 of Ref.[4]. Due to the 1/R variation of the toroidal magnetic ﬁeld (R is the major radius), the cylindrical Alfvén continuum (dashedlines) breaks up and the frequency gap (or the frequency forbidden band) appears within the range of0.308 < ω/ωA < 0.366.
For the simulation in the toroidal geometry,the numerical equilibrium from a Grad-Shafranov
solver[17] is employed. In Fig.3, shear Alfvén wave
oscillation is tested in a toroidal geometry in the absence of kinetic particles [the dynamical pressure termin Eq.(5) is absent]. In Fig.3 (a), the m = 1 and them = 2 oscillations are synchronized each other and(8)which is obtained by taking second order moment ofthe particle simulation. Here, I is the identity matrix
and b = B/|B| is the unit vector along the magnetic
ﬁeld line. Toroidal curvature eﬀects are included in
the second term of Eq. (5). As suggested in Eqs. (7)and (8), the second order moment of the kinetic particles from particle simulation replace the pressure evolution equation. In Eqs. (7) and (8), the subscript j
stands for the particle species.
The MHD simulation in this study is based on
up-down symmetry of the equilibrium. Given a sinefunction as an initial condition for the stream function
φ, the sine terms in ψ and the cosine terms in φ remain zero.[20] In the same manner the kinetic pressure
retains a parity of cosine.
The kinetic ions and electrons are given by gyrokinetic δf Vlasov equation[21]∂f0∂δfj= −Ẋ · ∇f0 − v˙�.
∂t
∂v�

(9)

The particle positions are time advanced using a second order Runge-Kutta-Gill method.[22]Ẋ = v�

cB�+b × (μ∇B − qj E� )B��
qj B��

(10)v˙� = −B�
· (μ∇B − qj E� )mj B��

(11)

where X is the position, v� is the velocity parallelto the equilibrium magnetic ﬁeld, μ = mv⊥ 2 /2B is
the magnetic moment. Here, qj is the charge of theparticles and c is the speed� of light. � For the ﬁeld
quantities, B� = B + ∇ × cmj v� /qj b, and E� =−∇φ− c−1 ∂t Ã. Note that the B� component contains
the conventional curvature drift term. Here, B�� =b · B� .

542

Y. Nishimura and C.Z. Cheng, Alfvén Wave Propagation in Initial Value MHD Simulation in a Tokamak Geometry

behave as one global mode.
The frequency spectrum in Fig.3(b) is obtainedby a Fourier transform of the time series in Fig.3(a).
In Fig.3 (b), we ﬁnd the frequency peak at one of theTAE gap boundaries (the lower boundary). In Fig.3(c), the eigenmode of perturbed ﬂux function ψ atthree diﬀerent incidents are plotted, which reveals apurely oscillatory behavior. These are good evidencesthat the Alfvén wave in our simulation is responding
correctly to the toroidal geometry (and thus TAE).
In Fig.4, we show our preliminary simulation results incorporating Eqs.(7)-(11). We push energeticparticles. The TAE mode is excited in the presence of energetic particle ions. The energetic particles are given by Gaussian distribution with thethermal velocity set to the Alfvén velocity. Here,βα = 4πpα /B02 = 0.01 is taken (pα is the pressureof the energetic particles). With the βα value andassuming the equilibrium ion density to be 1020 m−3 ,the Larmor radius of a deuterium energetic particle
at the Alfvén velocity is given by 1.82% of the minor
radius. In this work, the background density gradientis parametrized by a (∂ρ log neq ) = −1.0. Here, neq
is the equilibrium density proﬁle. The mode amplitude of the ψ ﬂuctuation experiences an exponential
growth of the envelope together with the oscillation at
the TAE frequency. The linear growth rate is givenby γ/ωA = 0.0027 while the frequency is given by
ω/ωA = 0.316. In Fig.4(b), the frequency peak resides within the TAE gap. Figure 4(c) shows the TAE
eigenmode structure in the presence of the energetic
particles.

(a)

1e−07

ρ=0.25
ρ=0.50ψ1/1

0−1e−07
1000

1500
ωA t(b)0.2

AnalyticalSimulation

ω/ωA

0.15

0.1

0.05

00

0.1

0.2

0.3

0.4

0.5

ρ

Fig. 1 Alfvén wave propagation in a cylindrical limit : (a)shear Alfvén wave oscillation (of the perturbed ﬂux
function ψ) at ρ = 0.25 and ρ = 0.50, (b) frequency
versus the radial location. Frequencies from thesimulation results match with the analytical prediction.
1continuum (cylinder)
continuum (torus)
gap boundaries

0.8

4. Summary and Discussions
In this work, we have discussed Alfvén wave propagation in cylindrical and toroidal geometries. We
have observed the mode frequencies at the TAE gap
in a toroidal geometry. Compared to the majority of
the MHD simulation work[9, 10, 11, 12, 13, 18, 19]
which have focused on the instability with ﬁnite dissipation in the system, the present work is unique inthat it reveals the Alfvén oscillation in an ideal MHD
limit (and by an initial value approach). Furthermore,
linear excitation of TAE mode is demonstrated in thepresence of energetic particles.
The advantage of the initial value approach is its
application for nonlinear simulation. The kinetic-ﬂuid
model[8] can be a useful tool to study nonlinear saturation mechanism of Alfvénic instabilities. To conserve energy in the nonlinear simulation, we need todescribe the the equilibrium distribution function interms of invariant canonical ﬂux surfaces.[24, 25]
We recapitulate the limitations given in thepresent paper. For the kinetic simulation, we did notincorporate particle drifts (we set μ = 0) nor ﬁnite

ω/ωA0.6
0.4
0.2
000.2

0.4

0.6

0.81

ρFig. 2 Shear Alfvén frequency as a function of the radiallocation. The dashed curves are the continuum frequencies for the cylindrical limit, for m = 1 andm = 2 modes. The solid lines are the continuum
frequency with the toroidal geometry eﬀect. The
lower (upper) boundary of the upper (lower) curveis at ω/ωA = 0.366 (ω/ωA = 0.308). The ﬁgure isa recapitulation of Fig.1 of Ref.[4].

543Y. Nishimura and C.Z. Cheng, Alfvén Wave Propagation in Initial Value MHD Simulation in a Tokamak Geometry

(a)

2e−06

m=1
m=2(a)

1e−05

m=1
m=2

1e−06
5e−06
ψm

0

ψm

0

−1e−06−5e−06−2e−06

0

100

200
ωAt

300

400

−1e−05200

400

1Frequency Spectrum (a.u.)Frequency Spectrum (a.u.)

0.8
0.6
0.4
0.2

0.3
ω/ωA

0.4

0.8
0.6
0.4
0.2
00.2

0.3
ω/ωA(c)

2e−06

1e−06

0.4

(c)

0.0001

ψ1/1800(b)

100.2

600
ωAt(b)

8e−056e−05
ψ1/1

04e−05−1e−062e−05−2e−06

00.2

0.4

0.6

0.8

10ρ

00.2

0.4

0.6

0.8

1ρFig. 3 (a) The Fourier components of the magnetic perturbation ψ as a function of time. (b) The frequency
spectrum obtained from the time series of Fig.3(a).
The peak is found at one of the boundary of thegap, ω/ωA = 0.308. (c) the TAE eigenmode structure in the absence of drive, which goes throughpure oscillation with time.

Fig. 4 (a) The Fourier components of the magnetic perturbation ψ as a function of time. (b) The frequency
spectrum obtained from the time series of Fig.4(a),(c) the TAE eigenmode structure of ψ1/1 in thepresence of energetic particles.

544

Y. Nishimura and C.Z. Cheng, Alfvén Wave Propagation in Initial Value MHD Simulation in a Tokamak Geometry

Larmor radius eﬀects. The thermal ions and electrons
are absent as well (in the present paper).
This work is supported by National Cheng Kung
University Top University Project. The authors would
like to thank Dr.K. C. Shaing for useful discussions.
[1] C. Z. Cheng, L. Chen, and M. S. Chance, Ann. Phys.
(N.Y.) 161, 21 (1985).
[2] C. Z. Cheng and M. S. Chance, Phys. Fluids 29, 3695
(1986).
[3] C. Z. Cheng, Phys. Fluids B 2, 1427 (1990);
C. Z. Cheng, Phys. Fluids B 3, 2463 (1991).
[4] G. Y. Fu and J. W. Van Dam, Phys. Fluids B 1, 1949
(1989).
[5] S. von Goeler, W. Stodiek, and N. Sauthoﬀ,Phys. Rev. Lett. 33, 1201 (1974).
[6] D. J. Campbell et al., Phys. Rev. Lett. 60, 2148 (1988).
[7] S. Barnabei et al., Phys. Rev. Lett. 84, 1212 (2000).
[8] C. Z. Cheng and J. R. Johnson, J. Geophys. Res. 104,413 (1999).
[9] C. Z. Cheng, J. Geophys. Res. 96, 21159 (1991).
[10] W. Park et al., Phys. Fluids B 4, 2033 (1992).
[11] G. Y. Fu and W. Park, Phys. Rev. Lett. 74, 1594
(1995).
[12] S. Briguglio, G. Vlad, F. Zonca, and C. Kar, Phys.
Plasmas 2, 3711 (1995).
[13] Y. Todo and T. Sato, Phys. Plasmas 5, 1321 (1998).
[14] C. Z. Cheng, N. N. Gorelenkov, and C. T. Hsu, Nucl.
Fusion 35, 1639 (1995).
[15] N. N. Gorelenkov, C. Z. Cheng, and W. M. Tang,Phys. Plasmas 5, 3389 (1998).
[16] H. R. Strauss, Phys. Fluids 20, 1354 (1977).
[17] V. E. Lynch, B. A. Carreras, H. R. Hicks,
J. A. Holmes, and L. Garcia, Comput. Phys. Commun.
24, 465 (1981).
[18] L. A. Charlton, J. A. Holmes, H. R. Hicks,V. E. Lynch, and B. A. Carreras, J. Comput. Phys.
86, 270 (1990).
[19] Y. Nishimura, J. D. Callen, and C. C. Hegna, Phys.
Plasmas 6, 4685 (1999).
[20] J. A. Holmes et al., Phys. Fluids 25, 800 (1982).
[21] R. G. Littlejohn, J. Plasma Phys. 29, 111 (1983).
[22] Originally a fourth order Runge-Kutta-Gill method
is employed; See Y. Nishimura et al., Contrib. Plasma
Phys. 48, 224 (2008); Y. Nishimura et al., JAERI-M,
93-140 (1993).
[23] A. M. Dimits and W. W. Lee, J. Comput. Phys. 107,309 (1993).
[24] Y. Idomura, S. Tokuda, and Y. Kishimoto, Nucl. Fusion 43, 234 (2003).
[25] J. Hollenberg and J. D. Callen, Phys. Plasmas 2, 1917
(1995).

545