ITR/P1-34
1
Benchmark of gyrokinetic, kinetic MHD and gyrofluid codes for the linear calculation of fast particle driven TAE dynamics A. K¨onies1 , S. Briguglio2 , N. Gorelenkov3 , T. Feh´er4 , M. Isaev5 , P. Lauber4 , A. Mishchenko1 , D. A. Spong6 , Y. Todo7 , W. A. Cooper8 , R. Hatzky4 , R. Kleiber1 , M. Borchardt1 , G. Vlad2 and ITPA EP TG9 1
Max Planck Institute for Plasma Physics, EURATOM Association, Wendelsteinstr. 1, 17491 Greifswald, Germany 2 Associazione EURATOM-ENEA sulla Fusione, C.P. 65-00044 Frascati, Rome, Italy 3 Princeton Plasma Physics Laboratory, Princeton University, Princeton, New Jersey 08543, USA 4 Max Planck Institute for Plasma Physics, EURATOM Association, Boltzmannstr. 2, 85748 Garching, Germany 5 Tokamak Physics Institute, NRC Kurchatov Institute, 123182 Moscow, Russia 6 Oak Ridge National Laboratory, Oak Ridge, Tennessee, USA 7 National Institute for Fusion Science,322-6 Oroshi-cho, Toki 509-5292, Japan 8 Ecole Polytechnique F´ed´erale de Lausanne, Centre de Recherches en Physique des Plasmas, Association Euratom-Suisse, CH1015 Lausanne, Switzerland Corresponding Author:
[email protected] Abstract: Fast particles in fusion plasmas may drive Alfv´en modes unstable leading to fluctuations of the internal electromagnetic fields and potential loss of particles. Such instabilities can have an impact on the performance and the wall-load of machines with burning plasmas such as ITER. A linear benchmark for a toroidal Alfv´en eigenmode (TAE) is done with about ten participating codes with a broad variation in the physical as well as the numerical models. A reasonable agreement of around ±20% has been found for the growth rates.
1
Introduction
Fast particles in ITER may originate from the fusion process itself or from external heating, such as Neutral Beam Injection (NBI). It is well known that those non-thermal populations of fast particles may interact with otherwise stable Alfv´en waves in the bulk plasma driving them unstable. This process takes place as a resonance phenomenon that requires a kinetic treatment of the fast particles but not necessarily a kinetic treatment of the bulk plasma.
ITR/P1-34
2
The oscillating electro-magnetic field in the plasma may lead to a loss of supra-thermal particles. As a consequence damage to in-vessel components of the machine is possible. In the last decades, much effort has been invested in the development of theory and codes that can be used to describe and explain the related phenomena. However, up to now, there is no well-understood standard case that these models have been tested against quantitatively. After studying mode damping [5], by providing the first comprehensive quantitative code comparison for the fast particle drive, the ITPA Energetic particle Topical Group is contributing to the design activity of the ITER operation scenario. The benchmark of the codes for the energetic ion driven modes is necessary for the accurate prediction of plasma behavior in ITER. The international benchmarking effort between a variety of codes shall ensure scientific quality and reliability when predictions for ITER are made.
2
Theory
As it has been outlined in the previous section, a kinetic description of fast particles or an appropriate closure of the fluid equations is necessary. In this section, different physical models in different implementations are described.
2.1
Kinetic MHD/ hybrid approach
Historically, the waves have been obtained from ideal MHD theory while the interaction with the wave has been treated subsequently with a drift or gyro-kinetic model. We call these type “kinetic MHD” or “hybrid” models. If the calculation is done linearly, the linear drift- or gyro-kinetic equation for the fast particles is solved in the given field of a precalculated MHD wave. The growth rate is calculated from the power transfer in relation to the energy stored in the mode (see e.g. [8]). In some codes, the kinetic equation is solved with a particle-in-cell (PIC) method. Here, numerical marker particles are distributed in phase space and represent a large number of real particles. As the Vlasov equation is of first order, it can be solved using the method of characteristics, i.e. by following particle orbits. PIC codes have advantages in difficult geometries such as stellarators but particle losses and noise issues have to be properly addressed. The following table summarizes the most important properties of the linear kinetic MHD codes participating in the benchmark: code NOVA-K CAS3D-K AE3D-K CKA-EUTERPE VENUS
MHD code NOVA CAS3D AE3D CKA CAS3D
MHD model ideal ideal reduced reduced ideal
2D 3D 3D 3D 3D
PIC no no yes yes yes
FOW yes no yes yes yes
FLR yes no no yes no
reference [4] [8] [12] [3] [2]
While in the first theoretical models the radial extent of the fast particle orbits has not been considered, meanwhile almost all codes consider the full orbit width (FOW). However, the effects of the finite Larmor radius (FLR) of the fast particles still are not commonly accounted for.
ITR/P1-34
3
Non-linear kinetic MHD hybrid codes, such as MEGA [13] and HMGC [1], solve the time-dependent non-linear MHD equations. The fast particles contribute to the plasma pressure and thus influence the evolution of the modes. This contribution is calculated from the non-linear drift- (HMGC) or gyro-kinetic (MEGA) equation for the fast particle species, where the MHD field enters as an external force.
2.2
Completely gyro-kinetic
There are also completely kinetic codes which solve the gyrokinetic equations for all ion species and the electrons together with the Poisson equation and Ampere’s law. The LIGKA code [9] is an eigenvalue solver which uses precalculated orbits from the HAGIS code to integrate the kinetic equations. Furthermore, there are fully kinetic PIC codes solving the time dependent kinetic equations with marker particles but with a consistently evolving field from the solution of the Maxwell equations. The GYGLES code [10] is two-dimensional in space and uses an ad-hoc equilibrium while the EUTERPE code ([7], still under construction) is able to use a realistic 3D equilibrium generated with the VMEC code [6]. Note, that the EUTERPE code is participating in this benchmark as a part of CKA-EUTERPE hybrid model, as well as a completely gyrokinetic model.
2.3
Gyro-fluid approach
The TAEFL code [11] solves the time dependent reduced MHD equations for the bulk plasma in a tokamak where the fast particles contribute with their pressure gradient. For the energetic particles a gyrofluid model with Landau closure is used. Due to the approximations used FOW and FLR effects are not included in the description.
3
Benchmark setup
A circular large aspect ratio tokamak (A = 10, R = 10 m) has been chosen as a test case. This was dictated by restrictions of the participating codes with respect to geometry or numerical properties. The minor radius is a = 1 m. The profile of the rotational transform is given by q(r) = 1.71 + 0.16(r/a)2, and the magnetic field is 3 T in the center. A TAE mode with m = 10, 11 and n = −6 is calculated in a hydrogen plasma with a flat density profile for both electrons and ions, i.e. ne = ni = 2.0 · 1019 m−3 , while Te = Ti = 1 keV holds for the electron and ion temperatures. Note that in fully kinetic codes like GYGLES and EUTERPE, the electron density is calculated from the quasi-neutrality condition if fast particles are present. That is not done within the hybrid codes. The pressure has to be zero at the boundary for an MHD equilibrium code. Therefore the bulk pressure decays towards the edge for the hybrid models: p(s) = (7.17 · 103 − 6.811 · 103 s − 3.585 · 102 s2 ) Pa while for the kinetic model the bulk pressure which can be calculated out of the background distribution functions is taken to be constant to avoid gradient driven modes to become important p(s) = pi + pe = 6408.0 Pa. Near the axis the values are quite close and the plasma beta is at around 0.2%.
ITR/P1-34
4 (b)
(a)
50
30
-1
20
3
γ/10 s
30
3
γ/10 s
-1
40
CAS3D-K (ZOW) GYGLES (ZLR) TAEFL (gyro fluid) AE3D-K (ZLR) NOVA-K (ZLR) HMGC (ZLR) MEGA (ZLR) CKA-EUTERPE (ZLR) VENUS (ZLR)
20 10 0 0
200
400 T/ keV
600
CAS3D-K (ZOW) GYGLES (FLR) CKA-EUTERPE (FLR) MEGA (FLR) NOVA-K (FLR) LIGKA (FLR) EUTERPE (FLR)
10
0
800
0
200
400 T/ keV
600
800
FIG. 1: Growth rates from calculations without FLR effects (a) and with FLR effects (b). The dashed line from CAS3D-K is valid in the limit of zero orbit width (small energies) and is shown for comparison. The shaded grey area marks the ±15% margin around the mean value for (a) and ±20% for (b). The influence of the fast particle pressure on the equilibrium has not taken into account in the MHD calculations. The fast particle (deuterons) distribution is taken to be a Maxwellian and is varied in the temperature range from 0 keV to 800 keV. The fast particle density profile (in m3 ) is given by √ s − c0 c2 (1) n(s) = n0 c3 exp − tanh c1 c2 with n0 = 1.44131 · 1017 m−3 , s = Ψtor /Ψtor (a) and the coefficients c0 = 0.49123, c1 = 0.298228, c2 = 0.198739, c3 = 0.521298.
4
Results
In Fig. 1a the growth rates have been calculated without the effects of finite Larmor radius. Although this limit is appropriate only for small energies, it provides a good comparison as not all participating codes are able to consider FLR effects. The TAEFL code deviates from the kinetic models but shows an energy scaling which is comparable to CAS3D-K as both codes do not consider FOW effects. It has been found from analysis of the kinetic MHD results of HMGC and MEGA that the resonances at vf /vA = 1/3 and at vf /vA = 1/5 contribute mostly to the energy transfer between fast particles and waves (Fig. 2).
4.1
Eigenfunctions
The qualitative agreement between the eigenfunctions is good: Some examples are shown in Fig. 3 and also in Fig. 4 for comparison.
ITR/P1-34
5
HMGC: power transfer (0.448< r/a 400 keV) the grid would need to be refined in order to resolve properly the lost-particle boundary that influences significantly the EP drive. This convergence study was not carried out here and therefore no LIGKA results for T > 400 keV are available. 4.3.3
Numerical damping in time dependent MHD codes
Those codes that solve time dependent fluid equations like MEGA and HMGC inherently have numerical damping. Varying the fast particle density in growth rate calculations, the damping value can be extrapolated. The growth rates in this work have been corrected by adding the damping value. Note, that damping rates vary with temperature. For example, for T = 200 keV without FLR, the damping rate from the MEGA code is −1.2 · 103 s−1 , and for T = 400 keV without FLR, −4.4 · 103 s−1 . They are comparable to the collisionless damping rate from LIGKA and GYGLES.
5
Conclusions
A linear benchmark for the calculation of the growth rate of a TAE mode has been performed by a number of codes using fully gyro-kinetic, kinetic MHD and gyro-fluid models. The importance of FLR effects has been illustrated, which is in agreement with earlier research [4]. The MEGA code has been upgraded for this benchmark, to include FLR effects and an extended version of HMGC is under construction [14]. The overall agreement of the codes is satisfactory for fast particle energies below 400keV and lies within ±20% for the codes including FLR effects and ±15% for those without such effects.
ITR/P1-34
8
For higher energies, the orbits become large and numerical problems (lost particles, orbits outside the plasma boundary) become more severe. With respect to calculations for burning machines with large orbits the question of lost particles should be addressed in further research, especially for PIC codes. The results mark the starting point for a non-linear benchmark and linear calculations of stability boundaries due to fast particle effects in ITER.
References [1] S. Briguglio, G. Vlad, F. Zonca, and C. Kar, Phys. Plasmas 2 (1995), 3711. [2] W. A. Cooper, J. P. Graves, S. Brunner, and M. Yu. Isaev, Plasma Phys. Contr. Fusion 53 (2011), no. 2, 024001. [3] T. Feher, A. K¨ onies, and R. Kleiber. IAEA TM on Energetic Particles Austin, TX (2011). [4] N. N. Gorelenkov, C. Z. Cheng, and G. Y. Fu, Phys. Plasmas 6 (1999), 2802–2807. [5] ITPA group on energetic particles, D. Borba, A. Fasoli, S. G¨ unter N. N. Gorelenkov, Ph. Lauber, N. Mellet, R. Nazikian, T. Panis, S. D. Pinches, D. Spong, D. Testa, and JET-EFDA contributors, 23st IAEA Fusion Energy Conference, 2010, pp. TWP/P7–08. [6] S. P. Hirshman, W. I. van Rij, and P. Merkel, Comp. Phys. Comm. 43 (1986), 143–155. [7] R. Kleiber, Theory of Fusion Plasmas, 2006, pp. 136–146. [8] A. K¨onies, A. Mishchenko, and Roman Hatzky, Theory of Fusion Plasmas, 2008, pp. 133–143. [9] Ph. Lauber, S. G¨ unter, A. K¨ onies, and S. D. Pinches, J. Comp. Phys. 226 (2007), 447–465. [10] A. Mishchenko, A. K¨ onies, and R. Hatzky, Phys. Plasmas 16 (2009), no. 8, 082105. [11] D. Spong, B. Carreras, and C. Hedrick, Phys. Fluids B 4 (1992), no. 10, 3316–3328. [12] D. Spong, E. D’Azevedo, and Y. Todo, Contr. Plasma Phys. 50 (2010), 708–712. [13] Y. Todo and T. Sato, Phys. Plasmas 5 (1998), 1321–1327. [14] X. Wang, S. Briguglio, L. Chen, C. Di Troia, G. Fogaccia, G. Vlad, and F. Zonca, Phys. Plasmas 18 (2011), no. 5, 052504.
