24D02
Kinetic Integrated Modeling of Toroidal Plasmas
トロイダルプラズマの運動論的統合モデリング
Atsushi Fukuyama, Sadayoshi Murakami, Hideo Nuga
福山 淳，村上定義，奴賀秀男
Department of Nuclear Engineering, Kyoto University
Yoshida-Honmachi, Sakyo-ku, Kyoto 606-8501, Japan
京都大学大学院工学研究科原子核工学専攻 606-8501 京都市左京区吉田本町
In order to self-consistently describe heating and current drive and various influences of energetic particles intoroidal plasmas, we have been developing a kinetic integrated modeling code TASK-3G. This modeling isbased on the behavior of the momentum distribution function of each particle species. The time evolution ofthe momentum distribution function is described by a Fokker-Planck component TASK/FP and the influenceof energetic particles on global stability is studied by a full wave component TASK/WM. Self-consistentanalysis of multi-scheme heating in a ITER plasma is demonstrated including radial transport and fusion
reaction rate calculated from the modified momentum distribution function. The linear stability of global
eigen modes in the presence of energetic particles is also discussed.

1. Introduction
In order to accurately predict the behavior of fusion plasmas and to develop reliable schemes controlling them, development of integrated simulation
codes for burning plasmas is urgently needed. Fusion reactions and external plasma heating and control generate energetic particles and modify the momentum distribution functions of plasma species.
The existence of these energetic particles may drive
global instabilities, such as resistive wall modes, internal kink modes, Alfvén eigenmodes and so on.
The occurrence of global instabilities usually redistributes the energetic particles. Previous analyses of momentum distribution functions usually assume background species with Maxwellian distribution and the temperature of the bulk component cannot be changed during the analyses. These assumptions are not satisfied in a strong heating case where
the bulk temperature changes rapidly. Therefore integrated core-transport simulation based on the momentum distribution is required for.
2. Kinetic Integrated Modeling of Tokamak Plasmas
For these purposes, we have been developing akinetic integrated simulation code, TASK-3G, basedon the time evolution of the momentum distributionfunctions. It is an extension of the integrated tokamak modeling code TASK [1] which includes thecomponents for MHD equilibrium, diﬀusive transport, ray tracing, full wave and the Fokker-Planck
analyses. In TASK-3G, the Fokker-Planck component TASK/FP describes the behavior of the bulk
component including the radial transport, and selfconsistently simulates strong heating cases where themomentum distribution function is strongly modified. The modification is taken into account in thefull wave component TASK/WM which describes
the ICRF heating and low-frequency global instabilities.
3. Multi-species Fokker-Planck analysis
The relativistic bounce-averaged Fokker-Planck
component TASK/FP [2] has been extended to describe the time evolution of the multi-species momentum distribution function f s (p, θ, ρ) where s, p,θ and ρ are particle species, magnitude of momentum, pitch angle and normalized minor radius, respectively. In this modeling, axisymmetry, time scale
longer than the particle bounce time, and zero bounce
orbit width are assumed. The Fokker-Planck equation includes nonlinear Coulomb collision, quasilinear wave-particle interaction, parallel electric fieldacceleration, radial diﬀusion, and particle source.
Figs. 1, 2 and 3 show typical simulation results
for multi-scheme heating in ITER plasma. We assumed a radial diﬀusion coeﬃcient Drr (ρ) = 0.1(1 +9ρ2 ) [ m2 /s] and an inward pinch term to keep the initial density profile. The ICRF waves heat tritons bysecond-harmonic cyclotron damping and electrons
by Landau damping. Deuterons are heated by NBI
and alpha particles are generated by DT fusion reaction. Fig. 1 indicates the power transfer between
diﬀerent species. Fig. 2 illustrates the contours off s at various radial position, 1 s after heating starts.
Initial electron density and temperature on the magnetic axis are 1020 m−3 and 20 keV. Fig. 3 indicate the radial profiles at t = 1 s of the absorbed
power density of each particle species, average kinetic energy of electrons, D and T, collisionally trans-0.25

e0.85

0.5

[1] A. Fukuyama et al., Proc. of 20th IAEA FEC (Villamoura, 2004) IAEA-CSP-25/CD/TH/P2-3.
[2] H. Nuga, A. Fukuyama, J. Plasma and Fusion Res.
Series, 8, 1125 (2009)[3] A. Fukuyama et al.: Proc. of 21th IAEA FEC(Chengdu, 2006) IAEA-CN-149/TH/P6-10.

DTHe121212

4

88

82

4

40-10

0-20

0-20

6

8

10-16

-12-8-4

0

4

8

1216

20

4-16-12-8-40

4

81216

20

0-20

1020

20

20

816

16

16

6

12

12

12

4

88

8

2

4

4-8-6-4-2

0

2

4

68

100-20-16-12-8-4

0

4

81216

20

0-20

-12-8-40

481216

20

0-20

20

20

20

816

16

166121212

4882

4

4-6-4-202468

100-20-16-12-8-4

0

4812

16

20

0-20-8-40

481216

20

0-20

20

20

81616

16

61212

12

4

88

8

2

4

40-10

0-20

0-20

-4-2

0

2468

10-16-12-8-4

0

4812

16

20-12-8-4

0

4812

16

20

0-20

20

20

20

816

16

16

6

12

12

12

4

8

8

8

2

4

4-8-6-4-2

0

24

68

100-20-16-12-8-4

0

4812

16

20

0-20-4048121620-16-12-8-4048121620-16-12-8-4048121620-16-12-8-4048121620-16-12-8-40481216204-16

100-10-84-1220-6-128-16

10-8-164-1610-8

1.0References64

ρ

Acknowledgments
This work is supported in part by Grant-in-Aid for
Scientific Research (S) (20226017) from JSPS.

16

20.5

4. Full wave analysis including the influence of energetic particles The ICRF wave electric field is calculated by the full wave component TASK/WM [3]in the case shown in Fig. 1. Using the kinetic dielectric tensor calculated by numerical integration in momentum space, TASK/WM can deal with arbitrary f s
in analyzing ICRF heating and current drive. Theenergetic particles also aﬀect the stability of Alfvéneigen modes and low-frequency global eigen modes.
The analysis of global eigen mode usually studied
with MHD model, e.g. resistive wall modes and internal kink modes, is planned including the influenceof energetic particles.

2000Fig. 3. Radial profiles of (a) absorbed power densityand (b) average kinetic energy

16-201.0ρ20-4

e1016-630

D0020-8T20

α80-100.65

1.0

100-100.45

40
T [keV]

Pabs [MW/m3]

DIn the calculation shown in Fig. 3, the radial diﬀusion broadens the kinetic energy density, especially T
heated by ICRF waves. The region of triton temperature increase becomes broader and is shifted inward
due to the inward pinch. We found that the heating
and current drive profiles are sensitive to the radialtransport model. We have introduced several kindsof turbulent diﬀusion coeﬃcient and√inward pinch,
without energy dependence, with 1/ 1 + p2 dependence, fixed radial profile. Validation of the radialtransport model with experimentally observed radial
profiles of energetic particles is under way.
Electron50

TFig. 1. Power transfer among plasma species

ρ
0.05

(b)

(a)

ferred power density, and time evolution of collisionally transferred power. The distribution f s is calculated on 50 magnetic surfaces by full-implicit parallel solver on 50 CPU cores. At t = 200 ms, thealpha particle density is only 0.5% of the electrondensity. The tritons produced by DD reaction aﬀect
fT near the normalized momentum p ∼ 10, about
1 MeV. Collisional power transfer between species
is also calculated self-consistently by using the nonlinear collision operator which conserves momentum
and energy with reasonable numerical accuracy.

4-16-12-8-40481216

200-20Fig. 2. Radial dependence of momentum distribution functions for electron, D, T, He.

2