第 13 回
核燃焼プラズマ
統合コード研究会
13th Burning Plasma Simulation
Initiative (BPSI) Meeting
日時：2015 年 12 月 10 日(木)－11 日(金)場所：九州大学筑紫キャンパス
応用力学研究所 2 階大会議室

第 13 回核燃焼プラズマ統合コード研究会
13th Burning Plasma Simulation Initiative (BPSI) Meeting(Ver.3)
日時：2015 年 12 月 10 日(木)－11 日(金)場所：九州大学筑紫キャンパス 応用力学研究所 2 階大会議室(18 min talk+7 min discuss or (15 or 10) min talk+5min discuss)

12 月 10 日(木)9:00 – 9:10 はじめに
福山 (京大)(座長：糟谷)9:10 – 9:35 講演 1-1
矢木 (原子力機構)
Simulation study on transport barrier in tokamak plasmas
9:35 – 10:00 講演 1-2
林 (原子力機構)
Integrated modeling of impurity transport in core and SOL/divertor
plasmas10:00 – 10:15 講演 1-3
前田 (京大)Heat transport simulation of H/He experiment plasma in LHD
10:15 – 10:40 講演 1-4
村上 (京大)Integrated simulation of deuterium plasma experiment of LHD
10:40 – 10:55 休憩(座長：小菅)10:55 – 11:20 講演 1-5
福山 (京大)Analysis of EC wave propagation in high-density tokamaks
11:20 – 11:35 講演 1-6
池田 (京大)Particle simulation of ECR plasma production in tokamaks
11:35 – 11:50 講演 1-7
二口 (京大)Kinetic linear stability analysis of Alfven eigenmode in tokamaks
11:50 – 13:00 昼休み(座長：本多)13:00 – 13:25 講演 2-1
登田 (核融合研)
Progress in applying gyrokinetic heat diffusivity model to transport
simulations for helical plasmas1

13:25 – 13:50 講演 2-2
沼波 (核融合研)
Sensitivity and prediction of plasma profile in turbulent transport ofhelical plasmas
13:50 – 14:15 講演 2-3
仲田 (核融合研)
Multi-species ITG-TEM turbulence simulations in ITER D-T-He plasma
14:15 – 14:40 講演 2-4
糟谷 (九大)Simulation study of turbulence and diagnostic in toroidal plasmas
14:40 – 14:55 休憩(座長：林)
14:55 – 15:20 講演 2-5
本多 (原子力機構)
Development of the fluid-type transport code on the flux coordinates intokamaks
15:20 – 15:45 講演 2-6
Nicolas (核融合研)
Bifurcation of the interchange mode growth rate and rotation frequency
due to the perpendicular heat conductivity
15:45 – 16:10 講演 2-7
奴賀 (京大)Fokker-Planck simulation for runaway electron generation including the
hot-tail effect
16:10 – 16:25 講演 2-8
永井 (京大)Burning plasma start-up simulation by the TASK code
16:25 – 16:40 休憩(座長：福山)16:40 – 17:05 講演 2-9
大澤 (九大)Hydrogen isotope trapped in tungsten vacancy
17:05 – 17:30 講演 2-10
佐々木 (九大)
Eigenmode analysis of geodesic acoustic modes driven by fast ions
17:30 – 17:55 講演 2-11
小菅 (九大)Analysis of triplet term in parallel momentum flux
17:55

散会

19:30 – 22:00 懇親会 (十徳や博多駅筑紫口店にて)

2

12 月 11 日(金)9:00 – 9:05 事務連絡
(座長：小関)9:05 – 9:30 講演 3-1
河村 (核融合研)Transport analysis of impurity in the LHD peripheral region during neon
gas-puff
9:30 – 9:55 講演 3-2
滝塚 (阪大)Bohm criterion and virtual diverter model for SOL-divertor simulation
9:55 – 10:15 講演 3-3
東郷 (東大)Numerical analysis of the effects of radial plasma diffusion on the
detached divertor plasmas
10:15 – 10:30 休憩(座長：矢木)
10:30 – 10:55 講演 3-4
横山 (核融合研)
Statistical analysis for deducing "transport model" in LHD plasmas(continued)
10:55 – 11:20 講演 3-5
小関 (原子力機構)
Predictive simulation of fusion plasma by data assimilation
11:20 – 11:45 講演 3-6
藤田 (名大)Development of PHZ pinch model with plasma rotation effects on particle
orbit
11:45 – 12:00 まとめ糟谷 (九大)12:00 – 13:00 昼休み
13:00 – 17:00 核融合エネルギーフォーラムサブクラスター会合
17:00

散会

3

Simulation study on transport barrier in tokamak plasmas
M. Yagi, H. Seto, Y. W. Cho1, and T. S. Hahm1
Japan Atomic Energy Agency, Rokkasho, Aomori, Japan
1

Seoul National University, Seoul 151-744, Republic of Korea

1. Introduction
To understand physical mechanism of transport barrier formation is a crucial issue fordesigning advanced operation scenario in JT-60SA and ITER. There are 4 approaches to solve
this problem. Namely, (1) 1 dimensional transport model combined with local gyro-kinetic
model such as GS2, TGLF. TRINITY code is an example for this approach[1]. (2) Two-fluid
model taking account of neoclassical transport. CENTRI code is an example for thisapproach[2]. However, it only covers MHD scale, practically (n < 20, where n is toroidal
mode number). (3) Global gyro-fluid model where the neoclassical transport is usuallysimplified. (4) Global gyro-kinetic model where full electron dynamics is still difficult tohandle. In JAEA, the project for multi-scale simulation with transport/MHD is going on. (1)
TRESS(transport code)+GKV[3], (2) LHTOK Project (IFERC-CSC): The neoclassical flows
are implemented into BOUT++ code to investigate L/H transition and ETB formation[4]. (3)MS_Turb Project (IFERC-CSC): The neoclassical flows are implemented into GF-ITG code
to revisit ITB formation investigated by Tokunaga et al[5]. (4) Near future. In this paper, thethird part is briefly reported.
2. Model equations
GF-ITG code consists of three field equations(ion continuity equation, parallel momentum
equation and ion temperature evolution equation)[5]
dW
1 ∂Φ
q 1 ∂
+ κn
+ A∇ //V = ε aω̂ d F + ρ* µ1(rU p ) − ρ*2 µ∇ ⊥4 Fdtr ∂ϑε r ∂r

(1)

dV
12
21
= −A∇ // F + 4µ∇ ⊥2 V − µ1U p − A
π k// V + A
∇ //T
dt
τ 5
5τ

(2)

3 # dT1 ∂Φ & # dn
1 ∂Φ & 5
92+ κTA ∇ // T + A∇ //V + χ ⊥ ∇ ⊥2 T + ST δ0,0 (3)%( − % + κn
( = ε aω̂ d T −2 $ dtr ∂ϑ ' $ dtr ∂ϑ ' 2τ55 π
where W = n − ∇ ⊥2 F,

n =Φ− Φ ,F = Φ+ p τ,

p = n + T . The magnetic curvature

operator is given byω̂ d = 2 cosϑ1 ∂∂+ 2sin ϑ
r ∂ϑ
∂r

(4)

The poloidal flow is given by
U p = V + ρ*q ∂( F − κ ncT )
ε ∂r

(5)

The neoclassical viscosities are given by[6,7]µ1 =0.66ε 1/2ν i(1+1.03ν *i1/2 + 0.31ν *i )(1+ 0.66ε 3/2ν *i )(6)%
1 " 1.17 − 0.35ν *i1/2− 2.1ν *i2 ε 3 '2 3$
1/2
1+ ν *iε # 1+ 0.7ν *i&(7)

κ nc =For the normalization, Bohm time t B = a 2 χ B , χ B = cTe eB and the minor radius a areused[5].
3. Analytical solutions of neoclassical flow
For simplicity, we assume the cylindrical geometry without turbulence to obtain the steadystate solution. The equations are given by'∂1 ∂ ∂
1 q 1 ∂ $q ∂r (Φ + T ) = −µ1
r & V + ρ*
Φ + (1− κ nc ) T )((∂t r ∂r ∂r
ρ* ε r ∂r %
ε ∂r

(8)

$
'
∂V
q ∂
= −µ1 &V + ρ*
Φ + (1− κ nc ) T )
(
%
(
∂t
ε ∂r

(9)

∂T
1 ∂ ∂T=χr+ ST∂tr ∂r ∂r(10)

where ST = χ̂ T0 J 0 ( j0,1r),2χ̂ = j0,1
χ,

j0,1 is the solution of J 0 ( j0,1 ) = 0 . The analyticalsolutions are written as##µ ncκ nc & − χ̂ t &Φ(r, t) = J 0 ( j0,1r) %% c1 + c2 e−( µ1+µnc )t + %1+( T0 e (($ χ̂ − µ1 − µ nc '$'V (r, t) = J1 ( j0,1r)ρ*

%
q "
ε2
µ1κ ncT0
j0,1 $ c1 + T0 (1− κ nc ) − 2 c2 e−( µ1+µnc )t −e− χ̂ t 'ε #qχ̂ − µ1 − µ nc&

T (r, t) = T0 J 0 ( j0,1r) (1− e− χ̂ t )

(11)
(12)
(13)

where µ nc = q 2 ε 2 µ1. Imposing initial conditions: Φ(r, 0) = V (r, 0) = 0, the steady statesolution of the electric field is given by
Er (r,∞) = −%∂κ nc (
Φ(r,∞) = j0,1J1 ( j0,1r)T0 ' −1+*∂r1+ ε 2 q 2 )&(14)

where κ nc = 1.17(banana), −0.5( plateau), −2.1(collisional) . The electric field will change thesign from positive to negative according to banana-plateau transition.
4. Simulation results

In nonlinear simulation by GF-ITG, we choose the source ST = 8 χ (1− 2r 2 ), and initial value
T0,0 (r, t = 0) = 0.5((1− r 2 )2 + T1 ), T1 = 10 −3 . Figure 1 shows q profile. The minimum point islocated at r = 0.5 . Figure 2 shows the time evolution of fluctuating internal energy. The
dominant Fourier modes are shown. It is seen that the energy of (0,0) and (1,0) modes
decrease. Figure 3 shows time evolution of mean ion temperature profile.

Fig.1 q profile

Fig.2 Time evolution of internal energy

The mean ion temperature decreases and the steady state is not obtained in this simulation.

Fig.3 Time evolution of mean ion

Fig.4 Contour plot of

temperature profile

T=6.

at

Figure 4 shows the contour plot of T (r, ϑ , z = 0) at T = 6 . Outside of q minimum point, thelarge structure develops. The heating source is used to maintain the large structure. Moresimulation study is necessary to understand flow effect on ITB formation.
5. Discussion
In the previous study, we have compared flux driven simulation and fixed gradient simulation
without neoclassical flows. In these cases, we have observed ITG growth and saturation forboth cases. In the present case, flux driven simulation gives no ITG growth instead low nmode growth and saturation. We will perform fixed gradient simulation with neoclassicalflows and investigate evolution of each mode. If we have saturation of ITG turbulence for thiscase, then we will apply the perturbative source and observe the plasma response. It will help
us to understand the effect of neoclassical flows on ITB formation. It is left for a future work.
Acknowledgments
The author would like to thank Drs. N. Miyato, A. Matsuyama and T. Takizuka for usefulcomments. The computation was carried out using the HELIOS supercomputer system at
IFERC-CSC.
References[1] M. Barnes et al., Physics of Plasmas 17 056109(2010).
[2] P. J. Knight et al., Computer Physics Communications 183, 2346(2012).
[3] M. Honda and M. Nakata, 32th JSPF annual meeting, Nagoya, 24aE06P(2015).
[4] H. Seto, M. Yagi and X.Q. Xu, 15th International Workshop on Plasma Edge Theory in
Fusion Devices, Nara, Japan, P1-10(2015).
[5] S. Tokunaga et al., Nucl. Fusion 49 075023(2009).
[6] J. D. Callen et al., Plasma Physics and Controlled Nuclear Fusion Research (International
Atomic Energy Agency, Vienna, 1987), Vol. I, p. 157; µ1 is given here.
[7] F. L. Hinton and R. D. Hazeltine, Rev. Mod. Phys. 48 239(1976); κ nc is given byeq.(6.136).

	
 Integrated modeling of impurity transport in core and SOL/divertor plasmas	
 N. Hayashi, K. Shimizu, K. Hoshino, M. Honda
Japan Atomic Energy Agency
Impurity seeding plays a key role to reduce the heat load on divertor plates in tokamak plasmas.
However, the accumulation of seeded impurity in the core may reduce the core confinement byenhancing the radiation and dilution. In order to study the transport of impurity in SOL/divertor andcore regions and its self-consistent effect on the divertor heat load and the core energy confinement,
the integrated modeling code for core and SOL/divertor plasmas including impurity is one of theeffective methods. For this integrated modeling, it is necessary to couple a core impurity transport
code with a SOL/divertor impurity transport code. In this paper, we have developed a model tocouple 1-dimensional (1D) fluid code for core impurities with 2D Monte-Carlo (MC) code for
SOL/divertor impurities.
We have been developing an integrated modeling code TOPICS (formerly called TOPICS-IB[1]), in which a 1.5D core transport solver is coupled with a 2D divertor code SONIC [2]. The 1.5D
core transport solver in TOPICS, however, treats only fully striped ions and can not treat many
charge states of impurity and resultant radiations. To treat core impurity physics in TOPICS, wedeveloped a 1D core impurity transport code IMPACT, which calculates ionization/recombination
processes, radiations, radial diffusive/convective transports of multi impurity species, and couple
IMPACT with TOPICS. One of the features of SONIC is the MC approach in an impurity-part code
IMPMC, which has a lot of flexibility in impurity modeling such as interactions with walls andkinetic effects. However, the IMPMC run including the core region takes much calculation time andthus requires coupling with IMPACT. At the interface boundary of IMPMC, only inward diffusive
flow is obtained because MC codes need a region beyond the interface to evaluate the net flow.
Outflow is a direct boundary condition for MC codes. On the other hand, the fluid code IMPACT
uses the net flow and thus can not use directory the inward diffusive flow of IMPMC as a boundary
condition. Therefore, we develop a interface model to couple the codes as shown in Fig.1. When the
diffusive transport is considered, low charge state Z impurities in a low temperature T region go tohigh-T region and high-Z impurities in a high-T region go to the low-T region, vice versa. Thus, thedependence of radial impurity particle flux on Z is expected in the upper figure of Fig.1. A
procedure to couple is as follows; (1) Evaluate impurity densities of each charge state near theboundary in IMPMC, Set them to boundary densities in IMPACT and Run IMPACT and Evaluate a
net flow F, (2) Evaluate inflow F− in IMPMC, (3) Evaluate outflow F+ from a relation of F = F+ −

F−, Give it to IMPMC and Run IMPMC. The
above procedure is iterated to obtain astationary state solution, or used in the time
evolution. To verify the applicability of thismodel, we calculate a simple case where Ar is
used as an impurity and only the constantanomalous diffusivity is assumed. Figure 2
shows the dependence of radial particle flux on
the charge state at the boundary of ρ=0.97
obtained from IMPACT. The Z dependence
almost agrees with the expectation in Fig.1, butdistorts largely. The peak at Z=16 due to its
largest accumulation in the core edge region.

Fig.1 Schematic view of dependence of radial
impurity particle flux on charge state Z (upper) andmodel procedure to couple impurity fluid code
IMPACT and impurity MC code IMPMC(bottom)

The obtained flux profile as well as profiles of
impurity densities almost agree with thoseobtained from IMPMC including the core region
beyond ρ=0.97 without the interface model. As aresult, the applicability of interface model isconfirmed in the IMPCAT side. Next, we will run
IMPMC with test particles generated at theboundary of IMPMC side on the basis of outflow
F+ in Fig.2 and confirm the model applicability in
the IMPMC side by comparing obtained resultswith those obtained from IMPMC including thecore region without the interface model. Then, we
will analyze more complicated cases including theFig.2 Dependence of radial Ar particle flux on
charge state Z at interface of ρ=0.97 in casewhere only anomalous diffusivity Dan=0.3m2/s
is assumed in whole region.

neoclassical transport.
References[1] N. Hayashi et al Phys. Plasmas 17 (2010) 056112.
[2] M. Yagi et al Contrib. Plasma Phys. 52 (2012) 372.

Integrated heat transport simulation of
the helium-hydrogen plasma of LHD by TASK3D

S. Maeta, S. Murakami, H. Yamaguchi, A. Fukuyama, K. Nagaokaa , H. Takahashia ,H. Nakanoa , M. Osakabea , K. Tanakaa , M. Yokoyamaa , and LHD Exp Groupa
Department of Nuclear Engineering, Kyoto University, Kyoto 615-8540, Japana National Institute for Fusion Science,322-6 Oroshi-cho, Toki 509-5292, JapanThe heat transport of the helium-hydrogen experiment plasma is studied to makeclear the ion species dependence on the heat transport. We evaluate the NBI heat deposition of multi-ion plasma by GNET code[1, 2] and study the heat transport by TASK3D.
It is found that the turbulent transport improvement is necessary to explain the increasein ion temperature in the high fraction helium plasma. We apply the turbulent transport
improvement models and obtain relatively good agreements.

Introduction
Various physical phenomena with wide ranges of time and space scales are connectedin the fusion plasma. Thus, in order to predict the performance of the fusion plasma, anintegrated simulation combining various physics models is required. Many integrated
simulation code has been developed for tokamak and helical plasmas. TASK3D is an
integrated simulation code for helical plasma and has been applied to the analyses ofthe plasma transport of LHD plasmas[3, 4].
On the other hand, the deuterium experiment project from 2017 is planned in LHD,
where the deuterium NBI heating beams with the power more than 30MW are injectedto the deuterium plasma. Main objects of this project are to make clear the isotope eﬀect
on the heat and particle transport in the helical plasma, and to study energetic particle
confinement in a helical magnetic configuration measuring triton burn-up neutrons.
In the 2014 campaign of LHD we have performed the helium-hydrogen plasma experiment to make clear the ion species dependence on the heat transport, which would
be related to the isotope eﬀect on the plasma confinement. We have found that the central ion temperature goes up as the helium fraction increases in the helium-hydrogenplasma.
In this paper we study the heat transport of the helium-hydrogen experiment plasma
to make clear the ion species dependence on the heat transport.

1

(a)

(b) H/(H+He)=0.73(c) H/(H+He)=0.34

Figure 1: Comparisons of the ion and electron temperatures of the simulation and experiment results. (a) the helium density rate dependency, (b) the low helium density ratecase, and (c) the high helium density rate case. Solid lines and dots show the experimentand simulation results, respectively.

Heat transport simulation of hydrogen-helium plasma
We evaluate the NBI heat deposition of multi-ion plasma by GNET code[1, 2] and study
the heat transport by TASK3D.
In the previous studies we found that the turbulent transport models for electron, thegyro-Bohm χTe B , and for ion, the gyro-Bohm(gradT) χTi B , show relatively good agreements with the experimental results of the hydrogen experiment plasma[4]. We, first,
apply these turbulent transport model to the helium-hydrogen plasma. Figure 1 shows
the heat transport simulation results based on the turbulent transport model for hydrogenplasma. We can see that the ion temperature does not agree with the experimental results as the helium density rate is increased. Thus, we considered three other turbulent
transport models; the conventional model, the Ze f f depending model and the strongly
Ae f f depending model.
To reproduce the experiment results of the ion temperature, a conventional transport
improvement model is assumed as χTi B = χTi B × A−α
e f f , where the eﬀective mass number∑
is given by Ae f f = i ni Ai /ne . Figure 2 shows that the of the ion temperatures of thesimulation and experiment results assuming the conventional improvement model with
α = 0.5 and 10. It is found that the conventional model can not fit the strong helium
density rate dependency observed in the experiment plasma.
Next, we assume the turbulent transport improvement model depending on the Ze f f ,which is introduce to explain the transport improvement in the high Ti plasma withthe carbon pellet injection[5]. Figure 3 shows the comparisons of the simulation andexperimental results assuming the Ze f f dependent turbulent improvement model. We
can see a slight diﬀerent tendency with the experiment results.
Finally, we assume a new model strongly depending on the eﬀective mass number2(a) α = 0.5(b) α = 1.0

Figure 2: Comparisons of the ion temperatures of the simulation and experiment results
assuming the conventional improvement model with α = 0.5 and 10.
χTi B = χTi B × exp(−kF(Ae f f − c)). We determined constant numbers c = 1.4 and k = 2.0
to minimize the diﬀerence between the temperatures of experiments and using leastsquares method. Figure 4 shows the comparisons of the simulation and experimental
results assuming the strongly Ae f f depending model. We obtain good agreements in the
H/He dependences and in radial profiles of T c .

Conclusions
We have developed the integrated simulation code TASK3D to predict the plasma performance in the 3D magnetic configuration plasma. We have studied the heat transportof H/He plasma by TASK3D comparing with the experimental results of LHD to clarify
the ion temperature increase and the isotope eﬀect. Three turbulent transport models
are assumed; the conventional model, the Ze f f depending model and the strongly Ae f fdepending model. A relatively good agreement is obtained with the strongly Ae f f depending model (c=1.4 and k=2.0) and Ze f f model.

References[1] H. Yamaguchi and S. Murakami, Plasma Fusion Res. 9 (2014) 3403127.
[2] S. Murakami, et al., Nuclear Fusion 46 (2006) S425.
[3] S. Murakami, et al., Plasma Phys. Control. Fusion 57 (2015) 054009.
[4] A. Sakai, et al., Plasma Fusion Res. 9 (2014) 3403124.
[5] S. Murakami, et al., Proc. IAEA Fusion Energy Conf. 2014, Russia (2014) TH/P638.
3

(a)(b) H/(H+He)=0.73(c) H/(H+He)=0.34

Figure 3: Comparisons of the ion and electron temperatures of the simulation and experiment results assuming the Ze f f dependent improvement model, (a) the helium density
rate dependency, (b) the low helium density rate case, and (c) the high helium density
rate case.

(a)(b) H/(H+He)=0.73(c) H/(H+He)=0.34

Figure 4: Comparisons of the ion and electron temperatures of the simulation and experiment results assuming the strongly Ae f f depending model, (a) the helium density
rate dependency, (b) the low helium density rate case, and (c) the high helium density
rate case.

4

Integrated Simulation of Deuterium Experiment Plasma in LHD
S. Murakami, H. Yamaguchi, M. Homma, S. Maeta, A. Sakai, A. Fukuyama, K. Nagaoka∗ ,H. Takahashi∗ , H. Nakano∗ , M. Osakabe∗ , M. Yokoyama∗ , K. Tanaka∗ , K. Ida∗ , M. Yoshinuma∗ ,M. Isobe∗ , H. Tomita∗∗ , K. Ogawa∗ and LHD Exp Group∗
Department of Nuclear Engineering, Kyoto Univ., Nishikyo, Kyoto 615-8540, Japan∗ National Institute for Fusion Science, 322-6 Oroshi, Toki 509-5292, Japan∗∗ Department of Quantum Engineering, Nagoya Univ., Nagoya, Japane-mail: murakami@nucleng.kyoto-u.ac.jp
The deuterium experiment plasma of LHD is investigated by applying the integrated simulation code
TASK3D and the 5-D drift kinetic equation solver GNET. (i) More than 20% of ion temperature
increment is obtained in the deuterium plasma (nD /nH + nD = 0.8) due to the isotope eﬀect
assuming the turbulent transport model based on the H/He plasma experiment of LHD. (ii) The
triton burn-up simulation shows the triton slowing down distribution and more than 7.0 × 1011
m−3 s−1 of 14MeV neutron production by the D-T fusion reaction at the plasma center.
Introduction
The deuterium experiment project from 2017 is planned in LHD, where the deuterium NBI
heating beams with the power more than 30MW are injected into the deuterium plasma. Principal
objects of this project are to clarify the isotope eﬀect on the heat and particle transport in the helicalplasma and to study energetic particle confinement in a helical magnetic configuration measuringtriton burn-up neutrons.
We have developed an integrated transport simulation code, TASK3D[1], and applied to the highion temperature plasma of LHD[2]. We obtain relatively good agreements between the simulated andexperimental radial profiles of the density and temperature in the steady state hydrogen plasmas.
Also, to study the kinetic behavior of the energetic particle in helical plasmas, we have developed the
GNET[3] code, in which the drift kinetic equation of energetic particles is solved in five-dimensional
phase space. In this paper, we predict the plasma performance of the deuterium experiment plasma
in LHD by applying the TASK3D and GNET codes.
Prediction of high-Ti deuterium experiment plasma
We perform the integrated heat transport simulation of deuterium plasma by TASK3D code. A
typical density profiles are assumed as ns (r/a) = ns0 {1 − (r/a)8 } where ns is the s-species density∑
profile satisfying the charge neutrality ne = i ni . We assume two central density cases; ne0 = 2.0
and 3.0 × 1019 m−3 , and the deuterium-hydrogen density ratio as nD /nD + nH = 0.8.
We evaluate the heat deposition profiles for the multi-ion species plasma (e, H, He, C) by using
GNET code. Three tangential NBI heating (#1-3) with the beam energy Eb = 180keV and theheating power PtN BI = 14MW, and two perpendicular NBI heating (#4 and 5) with Eb = 60−80keV
and PpN BI = 18MW are injected to the deuterium plasma of LHD. Figure 1-(left) shows the heat
deposition profiles of the electron, deuterium and hydrogen ions. We can see peaked heating profilesof the electron and deuterium ion heating.
1

One-dimensional (1-D) diﬀusive heat transport equation with multi-ion species (H, He, C) issolved using the heat deposition profiles by GNET. We assume that the heat transport consists of theneoclassical and turbulent transport. The neoclassical transport database, DGN/LHD, evaluates
the neoclassical heat transport for all species, and the radial electric field is determined by the
ambipolar condition of neoclassical electron and ion fluxes.
In our previous studies, we have assumed the gyro-Bohm model for the electron heat conduction
and the gyro-Bohm grad-Ti model for the ion one[4] as the turbulent transport model for thehydrogen plasma. However, it is necessary to consider the isotope eﬀect in the deuterium plasmacase. To solve this problem we have performed the H/He experiment in the 2014 campaign of
LHD. We have varied the H-He ratio, ξH = nH /nH + nHe , keeping the electron density profile. It
is found that the ion temperature increases as ξH decreases while the electron temperature shows
almost no change. To take into account this ion temperature increment we assume a turbulent
transport model, in which we multiply an improvement factor depending on the eﬀective mass∑
numberAef f = i (Mi /Zi ) to the gyro-Bohm grad-Ti model.
Figure 1-(right) shows the deuterium ion temperatures in the deuterium experiment plasma of
LHD with and without the isotope eﬀect on the turbulent transport by the Aef f depending model.
It is found that the deuterium ion temperature reaches more than 6 keV with the isotope eﬀect inthe deuterium experiment plasma. On the other hand, the ion temperature reaches about 5 keV if
we assume a pure hydrogen plasma. This result indicates that we will obtain about 20% higher ion
temperature than that of the hydrogen plasma in the deuterium experiment of LHD if we assume
an isotope eﬀect on the turbulent transport based on the He/H experiment results.
Triton burn-up simulation
We perform the triton burn-up simulation of the deuterium experiment of LHD and evaluate
the D-T fusion reaction rates to compare with the experimental results of the 14 MeV neutron
diagnostic system. The radial profile of the 1 MeV triton production rate due to the D-D fusion
reaction between D-beam and D-thermal ions is evaluated by GNET using the radial profiles of TD
and nD obtained by TASK3D. As a result we obtain a center peaked triton production rate about
1.4 × 1014 m−3 s−1 in the ne0 = 2.0 × 1019 m−3 case.

Figure 1: Radial profiles of the NBI heat depositions (left) and ion temperatures (right) in the deuterium
experiment plasma with ne0 = 2.0 × 1019 m−3 and nD /nH + nD = 0.8.

2

-1

-0.5

0

0.5

6

19 -3

ne=2.0×1019m-3
ne=3.0×10 m

5

4

-3 -1

T(d,α)n fusion rates [10 cm s ]v⊥/v1MeV

28
26
24
22
20
18
16
14

1
0.8
0.6
0.4
0.2
0
1

4
3
2
1
0
0

v||/v1MeV

0.2

0.4

0.6

0.8

1

r/a

Figure 2: Slowing down distribution of 1MeV tritons in the velocity space (left) and radial profiles of 14MeV
neutron production rate by the D-T fusion reaction[PtN BI = 1MW] (right).

Next we solve the 5-D drift kinetic equation for the 1MeV tritons by GNET and evaluate theslow down distribution in the real and velocity space. Figure 2-(left) shows the slowing down velocity
space distribution of the 1MeV tritons in LHD. We can see a large loss region near the helicallytrapped particles. About 25% of tritons escape with almost the initial energy of 1 MeV by theprompt orbit loss due to drift motion immediately after their birth (t < 10−5 s). After the prompt
orbit loss, the collisionless diﬀusive loss (10−5 < t < 10−2 ) and, then, the collisional diﬀusive loss(t > 10−2 s) become dominant. We evaluate the radial profile of 14MeV neutron production rateusing the triton distribution by GNET and TD , nD profiles by TASK3D [Fig. 2-(right)]. It is foundthat more than 7.0 × 1011 m−3 s−1 of 14MeV neutrons are generated by the D-T fusion reaction at
the plasma center.

[1] M. Sato et al., Plasma Fusion Res. 3, S1063 (2008).
[2] S. Murakami et al., Plasma Phys. Control. Fusion 57, 054009 (2015).
[3] S. Murakami et al., Nuclear Fusion 46, S425 (2006).
[4] Wakasa, A. et al., Proc. 39th EPS Conf. and 16th Int. Cong. Plasma Phys., 2012, P2.028.

3

Particle simulation of plasma production by electron cyclotron waves
in Tokamak
Tomohide Ikeda, Atsushi Fukuyama
Graduate School of Engineering, Kyoto university

1

Introuction

Plasma production and start up by electron cyclotron waves in tokamak has been successfully applied
to reduce the magnetic flux comsumption and the peak value of toroidal electric field, though the physics
mechanism has not been understood well. Various phenomena are coupled with each other - electron accelation at the electron cyclotron resonance, ionization by collisions with neutral particles, current drive
by toroidal electric field and electron cyclotron waves, suppresion of particle losses by vertical magnetic
field, generation of closed magnetic configuration, and so on. In order to analyze these pheomena selfconsistently, we are developing a particle simulational module, TASK/PIC. In this report, we show thestatus of the two-dimentions electromagnetic particle simulation code and the results on the excitationand propagation of the electron cyclotron waves.
2 Numerical model
The governing equations in TASK/PIC are the equation of motion for electron and ion (1), Gauss’s
law for the scalar potential ϕ (2), time evolution of the vector potential A (3) and equations for the electricfield E and the magnetic field B (4):mdu
= q(E + u × B),
dt

(1)

ρ
,
ϵ0

(2)

−∇2 ϕ =

∇2 A + µ0 j − µ0 ϵ0
E = −∇ϕ −

∂A
,
∂t

∂∇ϕ
∂2 A
− µ0 ϵ0 2 = 0,∂t∂t
B = ∇ × A.

(3)

(4)

We normalize the physical quantities by the electron plasma frequency ω pe for the average electron
density n0 , ω2pe = n0 e2 /me ϵ0 , and the grid spacing ∆ of the order of the Debye length.
In our alogrithm, ū and Ā are defined at the integer mesh in time, while x̄, ϕ̄, Ē and B̄ are at thehalf-integer mesh in time. We assume uniformity in the z direction and solve eq.(2) by FFT or finite
diﬀerence method. For the particle shape, we use the spline function of the second order. The spline
function of the first order is used for the half-integer grid, and the spline function of the second order in
the integer grid. We consider a rectangle space in the x-y plane, and employ the boundary condition inperiodic or of the conducting wall.

3 Simulation result
The initial particle distribution is uniform in space, and the velocity distribution are Maxwellian.
Simulation parameters are as follows: total number of particles for each species (electron, ion): 160000,
number of grids: 100 × 100, β ≡ ω pe ∆/c = 0.1, the time interval ∆t = 0.02. Further, the externallyapplied fields are as follows: the static magnetic field with linear gradient B̄y0 = B̄ymin + (Bymax − Bymin ) ·x/100 and the wave guide elctric field at x = 0 and 25 < y < 75, Āywg (Āzwg ) = Āwg sin ω̄t¯ whichcorrespond to Ēywg (Ēzwg ) = −∂Āywg (Āzwg )/∂t¯ = −ω̄Āwg cos ω̄t¯.
First, we show the time evolution of total energy with periodic boundary condition up to t¯ =Fig.1(a): B̄ = 01000 in Fig.1. When By0 = 0, the potential energy PEtot increase and the kinetic energy KEtot
decrease up to around t¯ = 200. After that they almost saturated, but the total energy Etot gradually
increase from the begining. When uniform magnetic field B = 6, PEtot , KEtot and Etot gradually
increase all along. We have to resolve this problem.
Next, we excite the electromagnetic waves
with Āywg = Āwg sin ω̄t¯ (O-mode) or Āzwg =Fig.1(b): B̄ = 6
Āwg sin ω̄t¯ (X-mode) with Āwg = 1 and B̄ = 6.
In Fig.2, in the case of X-mode, absorption of theelctromagnetic wave occurs near the electron cyclotron resonance ω̄ = B̄ = 6. As is theoretically
predicted for the ordinary-mode, the absorption
does not occur. It is because finite Larmor radiuseﬀect is small with the present parameters. If wemake ∆ smaller than cyclotron radius vte /ωce , theabsorption will be confirmed.
Fig.1: Time evolution of total energy
with periodic boundary conditiont¯ = 2100

O-modeEy

MAX = 6.4996E+00

t¯ = 4
STEP = 1.4444E+00

100

Ey

MIN = -8.0000E+00

MAX = 8.0000E+00

t¯ = 6
STEP = 1.7778E+00

100

Ey

MIN = -8.0000E+00

MAX = 8.0000E+00

t¯ = 8STEP = 1.7778E+00

100

Ey

MIN = -6.0164E+00

MAX = 6.0164E+00

t¯ = 10
STEP = 1.3370E+00

100

80

80

80

80

80

60

60

60

60

60

40

40

40

40

40

20

20

20

20

20

0

0

0

0

0

100

X-modeMIN = -6.4996E+00

Ez

20

40

MIN = -6.0000E+00

60
MAX = 6.0000E+00

80

100

STEP = 1.3333E+00

0

100

Ez

20

40

MIN = -6.5747E+00

60
MAX = 6.5747E+00

80

100

STEP = 1.4610E+00

0

100

Ez

20

40

MIN = -6.6722E+00

60
MAX = 6.6722E+00

80

100

STEP = 1.4827E+00

Ez

20

40

MIN = -6.5216E+00

60
MAX = 6.5216E+00

80

100

STEP = 1.4493E+00

0
100

80

80

80

80

80

60

60

60

60

60

40

40

40

40

40

20

20

20

20

20

0

0

0

0

0

20

40

60

80

100

0

20

40

60

80

100

0

20

40

60

80

100

MIN = -6.0087E+00

MAX = 6.0087E+00

STEP = 1.3353E+00

0

0

100

Ey

0

20

40

60

80

100

0

Ez

0

20

40

MIN = -6.8453E+00

20

40

60
MAX = 6.8453E+00

60

80

100

STEP = 1.5212E+00

80

100

Fig.2: Wave propagation with ω̄ = 6 and B̄ = 4 ∼ 8Finally, we show the time evolution of electron temperature profile in Fig.3. T ∥ is the parallel components of the electron temperature, T ⊥ is the perpendicular component of the electron temperature and T e
is the electron temperature T e = (T ∥ + 2T ⊥ )/3. In the case of X-mode excitation, T ⊥ around x = 40 ∼ 50
increases up to about 7 times the initial temperature.

O-modeT ∥ (x)T ⊥ (x)T e (x)Tperpe(x)Tparae(x)

Te(x)

0.8
1.0

0.3
0.6

0.8

0.2

0.6

0.4

0.4
0.1

0.2
0.2

0.0

0.0
0

20

40

60

80

100

0.0
0

20

40

60

80

100

0

20

40

60

80

100

60

80

100

X-modeTparae(x)T ∥ (x)Tperpe(x)T ⊥ (x)6

0.25

T e (x)

Te(x)
6

0.20
0.15
0.10

4

4

2

2

0.05
0.00

0

20

40

60

80

100

0

0

20

40

60

80

100

0

0

20

40

Fig.3: Time evolution of electron temperature profile
4 Summary and Outlook
We have been developing a two-dimentional electromagnetic PIC simulation code TASK/PIC, inorder to analayze the plasma production near the electron cyclotron resonance. We use the parallel processing by MPI only for particle motion. After discussion on the conservation of total energy, we haveshown the propagation of the electromagnetic waves excited by the waveguide model. We have shown
the propagation and absorption of ordinaly and extraordinaly waves propagation and confirmed electron
heating near the electron cyclotron resonance. Owing to small finite Larmor radius eﬀect, the absorption
cannot be confirmed for ordinary mode.
The remaining issues include to improve phisycal models, such as relativistic eﬀect, absorbing
boundary condition, plasma production by ionization, to enhance computional performance, full-scaleparallelization, introdution of region splitting, and the use of OpenMP, and to extend to toroidal configuration.

Linear kinetic analysis of Alfvén eigemodes in tokamak plasmas
トカマクプラズマにおけるアルヴェン固有モードの運動論的線形安定性解析
Taisei Futakuchi, Atsushi Fukuyama
二口泰成，福山 淳
Graduate School of Engineering, Kyoto University
Kyoto-Daigaku-Katsura, Nishikyo-ku, Kyoto 615-8540, Japan
京都大学大学院工学研究科 615-8540 京都市西京区京都大学桂The excitation of Alfvén eigenmodes driven by fast ions is studied in tokamak plasmas using the full wave
code TASK/WM. The dielectric tensor for species with anisotropic velocity distribution function is numerically calculated and complex eigenmode frequency is obtained.

1

taken into account to describe the wave-particle resonant interaction. Eigenmode with complex wave
frequency is obtained by maximizing the wave amplitude for given external current proportional to theplasma pressure.
Since the response of plasmas is determined by
the dielectric tensor, the model of the dielectric tensor is important. However, the more general themodel is, the less tractable.
For simplicity, the model based on the drift kinetic equation system[1] is adopted in the present
analyses. This model which consists of the linear
current density and the drift kinetic equation can describe the poloidal mode coupling due to the interaction between drift motion of the fast particles andAEs. The induced current density is calculated by integrating perturbed distribution function of fast ionsin velocity space∫jp =dv eα vd fα1 = ↔
σ E = −i0 ω↔
χE(2)Introduction

In tokamak plasmas, Alfvén eigenmodes are excitedby fast ions generated as a result of fusion reaction,
ICRF heating or NBI heating. Resonant interaction
between fast ions and the destabilized Alfvén eigenmodes (AEs) may cause rapid transport of fast ions,which leads to a loss of fast ions, decreasing heatingeﬃciency and damaging plasma facing components.
Therefore stability analysis of Alfvén eigenmodes isimportant and it has been extensively explored theoretically, experimentally and numerically.
There are various AEs which have diﬀerent frequencies, such as toroidicity-induced Alfvén eigenmode (TAE), beta induced AE (BAE) and so on.
Some of them are aﬀected by mode conversion (e.g.
kinetic Alfveén eigenmode (KAE)) and lower frequency modes are coupled with other waves (e.g. ion
acoustic wave). Therefore careful study of AE using
appropriate model which can properly describe thedynamics in the focused scale is required.
We have analyzed the linear stability of AEs by
the full wave code TASK/WM. The model used in
the analysis and some examples of analyses with thecode is presented in this report.

and the induced current is proportional to the wave
electric field in a linear analysis. Then, the electricsusceptibility tensor ↔
χ is calculated by using perturbed distribution function which is described by thedrift kinetic equation

2 Full wave analysis code TASK/WM[

In the full wave code TASK/WM, Maxwell’s equation for the wave electric field E with frequency ωω2∇×∇×E= 2 ↔
 · E + i ωµ0 jext
c

]
∂
eα
∂
+ vk ∇k + (vd + vE ) · ∇ +(vk Ek + vd · E)
fα = 0,∂tmα∂εwhere vE is the E × B drift velocity and vd is themagnetic drift velocity. It is necessary to linearize
this equation in terms of the perturbed distributionfunction. Since it is a system of diﬀerential equations
including the term vd · ∇ f , it is complicated to solve
the equation directly. Even if the radial derivative(1)

is solved with the kinetic dielectric tensor ↔
 forvarious plasma models. With poloidal and toroidal
mode expansion, parallel wave number is correctly1term of vd · ∇ f is neglected, the drift kinetic equation
is infinite dimensional system of equations in Fourier
expansion. Considering only several sideband harmonics by truncating terms satisfying k · vd /ω  1,or assuming low drift frequency, the problem becomes tractable. Based on these assumptions, approximate electric susceptibility tensor ↔
χ is calculated in Fourier expansion form as
↔
χ mm0
+
+
+
+

The mode whose frequency is around 80kHz is
toroidaicity-induced Alfvén eigenmode (TAE) andits growth rate has good agreement (Fig.3) with theone calculated by approximate formula by G. Y. Fu
[4]
[ (
) ( )
]
9
ω∗αm 1
vA
vA
'
βα
−
F
− βe
. (6)
ωT AE 4
ωT AE 2
vT α
ve
γ



i
0 
 1
=  i −1 0  Pm−1,m−2 δm0 ,m−2
0
0
0


0
0
i Qm−1,m−1 00− Qm−1,m−1  δm0 ,m−1i Qm,m−1 − Qm,m−10i (Pm−1,m − Pm+1,m )0 −(Pm−1,m + Pm+1,m ) −i (Pm−1,m − Pm+1,m ) −(Pm−1,m + Pm+1,m ) δm0 ,m
0
0
0
−Rm,m


0
0
−i Qm+1,m+1 


0
0
− Qm+1,m+1  δm0 ,m+1
−i Qm,m+1 − Qm,m+1
0


1
−i
0
 −i −1 0  P m+1,m+2 δm0 ,m+2000In addition, it is shown that lower frequency modes
are more strongly excited by the increase of fast ion
temperature and the mode frequency having the highest growth rate gets close to TAE frequency. Thismode is considered as a kind of energetic particle
mode (EPM).

0.008

Pm,m∗

=

Qm,m∗

=

Rm

=

a

=

cm

=

∫
1
a
dv
vd 2 cm∗ fα0
4
1 − kkm vk /ω∫a1
dv
vd vk cm∗ fα0
2
1 − kkm vk /ω
∫1
vk 2 cm fα0
a dv
1 − kkm vk /ωIm(f) [MHz]

0.004

(3)

0.000

-0.004(4)-0.008
0.04

(5)

0.06

0.08

0.10

Re(f) [MHz]e2α,2ω 0 T αω∗m∂ ln fα0
iT α ∂ ln fα0 ∂T̂ −= −T α−ω∂εrωeα B ∂r ∂θ

Fig. 1 Contour of the wave amplitude in a complex
frequency plane with fast ion components 0.2%，
T= 100keV

In the case of steady state (Maxwellian distribution), Eq.(3) to (5) are easily calculated with plasma
dispersion relation. In the case that fast ion has
anisotropic distribution function or localized density
profile, it is necessary to integrate numerically incomplex velocity space.
In order to consider kinetic eﬀects such as finite
Larmor radius (FLR) eﬀects and coupling of Alfvénwaves with low frequency waves, the model based on
the gyrokinetic equation is required[2],[3].

0.008Im(f) [MHz]

0.004

0.000

-0.004-0.0080.04

3 Simulation results

0.06

0.08

0.10

Re(f) [MHz]Fig.1 and Fig.2 are examples of analyses by
TASK/WM describing excitation by fast ions in thecase of Maxwellian distribution function. These figures are contour of |E|2 in complex frequency space
and we can find eigenmode frequencies from them.

Fig. 2 Contour of the wave amplitude in a complex
frequency plane with fast ion components 0.2%，
T= 500keV
2

4

Future work

We are now developing a code to analyze AEs excited by fast ions which have arbitrary velocity distribution function. Then, it is possible to carry out
self-consistent analysis of AEs including the modification of velocity distribution function calculated by
the Fokker-Plank code.
REFERENCES[1] Y.M. Li, S.M. Mahajan and D.W. Ross, Phys.
Fluids 30, 1466 (1987)[2] L. Chen and F. Zonca, Nucl. Fusion 47, S727
(2007)[3] L. Chen and F. Zonca, EPL 96, 35001 (2011)[4] G. Y. Fu and J. W. Van Dam, Phys. Fluids B 1,2404 (1989)

Fig. 3 Comparison of the fast ion temperaturedependence of the growth rate between the
numerical results and these of the approximate
fomula.

3

Progress in applying gyrokinetic heat diﬀusivity model to transport
simulations for helical plasmas
S. Toda, M. Nunami, A. Ishizawa1 , S. Murakami2 , T. -H. Watanabe3 and H. Sugama
National Institute for Fusion Science, Oroshi-cho 322-6, Toki, Gifu, 509-5292, Japan
1 Graduate School of Energy Science, Kyoto Uni., Gokasho, Uji, Kyoto, 611-0011, Japan
2 Department of Nuclear Engineering, Kyoto Uni., Nishikyo, Kyoto, 615-8540, Japan
3 Department of Physics, Nagoya University, Furo-cho, Nagoya, Aichi, 464-8602, Japan

Turbulent transport is one of the most critical issues for plasma confinement in magnetic fusion
devices. Recently, a large number of the gyrokinetic simulations have been done in toroidal plasmas [1]. The GKV-X code solving the gyro-kinetic equation has been used to examine the ITG
mode and zonal flows in the LHD for studying the turbulent transport in helical plasmas [2]. The
reduced model of χi ∼ ρ2ti vti f (L, τ̃ZF )/R is taken [3] using the GKV-X code for the high-Ti LHD
discharge. Here, L is the mixing length estimate γ̃k̃y /k̃y2 integrated over the k̃y space, where γ̃k̃yis the normalized linear growth rate of the ITG mode and τ̃ZF is the normalized decay time ofzonal flows [4, 5]. However, it is costly to carry out linear calculations of the growth rate by the
gyrokinetic simulation at each time step of the dynamical transport code such as TASK3D [6],because the transport analysis of helical plasmas demands a high radial resolution so as to accurately evaluate the radial electric field and the field configuration. In this study, how to apply the
reduced model of the turbulent heat diﬀusivity for the ITG mode derived from the gyrokinetic
simulation to the transport code is shown with a low computational cost. Modeling of the term L
in the reduced model for the ITG mode is necessary to be involved with a parameter dependenceof the plasma instability in the dynamical transport code. The ion temperature gradient scale
length LTi (= −Ti /(∂Ti /∂r)) is chosen for the parameter to apply L to avoid the linear gyrokinetic
analysis at each time step of the transport simulation. The field configuration is fixed at the initialstate in the transport simulation. We have developed the formula for L in terms of LTi . The decay
time of zonal flows depends on the magnetic field configuration and is independent of the plasma
profile. The formula of the zonal flow decay time is needed to be calculated only at the initial state.
The calculation by substituting these formulae for L and τ̃ZF to the reduced model [3] reproduces
the results of the reduced model itself within allowable errors. The computational cost to obtainthe value of the turbulent ion heat diﬀusivity by this modeling at each time step of the transport
simulation is much smaller than that of the linear gyrokinetic simulation. The simulation results
for the high-Ti discharge of the shot number 88343 have been shown [7]. In this study, the modelingof the mixing length estimate term and the zonal flow decay time is done for the high-Ti dischargeof #109081. The value of the density gradient scale length in the case of #109081 is diﬀerent fromthat in the case of #88343. The density gradient scale length aﬀects the characteristic of the ITG
instability. This additional modeling is applied to the transport code and enables us to study thesimulation results with the experimental results in LHD.
Firstly, the linear analysis is done using the GKV-X code for the additional modeling of the turbulent ion heat diﬀusivity. The ITG instability is examined in the high-Ti LHD discharge #88343 [8].
On the other hand, the analysis of the plasma dynamics using TASK3D was done for #109081
in (LHD [9]. A reduced
model for the ITG turbulent heat diﬀusivity in terms of the functions)∫2= A1 Lα /(A2 + τ̃ZF /L1/2 ),L ≡ (γ̃k̃y /k̃y )dk̃y and τ̃ZF (= τZF /(R/vti )) was obtained as χi /χGBiis the gyro-Bohm diﬀusivity, γ̃ = γ/(vti /R) and k̃y = ky ρti . The numerical coeﬃcientswhere χGBi
are given in [3]. The characteristic length of the ion temperature gradient is considered to be theimportant parameter for the ITG instability. As the function of LTi , the parameter L is modeled

Figure 1: The radial dependence of (a) R/LTc and (b) a(ρ) is shown with filled circles for #109081
in LHD.
by(R
R
L = a(ρ)
−
LTi
LTc

)

,

(1)

where LTc is the normalized critical ion temperature gradient for the ITG instability. To examine the critical ion temperature gradient for the ITG mode, the dependence of L on R/LTi is
examined [10] with all plasma parameters fixed except the ion temperature gradient. The density and ion temperature profiles are obtained from the experimental results at t = 2.233s inthe high-Ti LHD discharge #88343. The critical ion temperature gradient for the ITG mode,
R/LTc and the slope a(ρ) in terms of R/LTi are obtained. The fitting polynomials of R/LTc
and a(ρ) are used as R/LTc = 4.0929 − 3.7681ρ + 19.712ρ2 + 11.087ρ3 − 14.272ρ4 and a(ρ) =0.38661 − 0.070919ρ + 0.2571ρ2 + 0.95949ρ3 − 0.92978ρ4 . The comparison between the right handside (a(ρ) (R/LTi − R/LTc )) and the left hand side (L) in equation (1) is examined with the root
mean square of ((a(ρ) (R/LTi − R/LTc ))/L − 1) given by σ = 0.13. The zonal flow decay timeτ̃ZF [5] is examined. The fitting function for the zonal flow decay time: τ̃ZF (f it) = 0.98565 −
0.65943ρ + 2.4471ρ2 + 3.2337ρ3 − 2.8382ρ4 is used throughout the transport simulation. When we
use the same modeling for the mixing length estimate and the zonal flow decay time at t = 4.2s for#109081 in LHD, the fitting polynomials of (a) the normalized ion temperature gradient R/LTc and(b) a(ρ) are obtained as R/LTc = 4.712−38.644ρ+512.56ρ2 −2481ρ3 +5991.3ρ4 −6810.3ρ5 +2911.3ρ6
and a(ρ) = 0.13049 + 6.3631ρ − 45.33ρ2 + 134.8ρ3 − 177.95ρ4 + 90.153ρ5 − 4.6622ρ6 in figure 1. The
fitting function for τ̃ZF : τ̃ZF (f it) = 0.34823 + 6.817ρ − 15.452ρ2 + 15.972ρ3 − 6.4126ρ4 is shown forthe case of #109081 in LHD. To obtain the value of the turbulent ion heat diﬀusivity at each time
step in the transport code, the value of the modeled turbulent ion heat diﬀusivity is calculated,
where the term a(ρ)(R/LTi − R/LTc ) is substituted for L and the fitting function τ̃ZF (f it) is usedfor τ̃ZF .
Now, the transport dynamics is examined using the modeled turbulent ion heat diﬀusivity, when
the TASK3D is performed. The radial profiles of the density and the electron temperature are
fixed. The dynamics of the radial Ti profile is simulated by solving the diﬀusion equation. The
profile of the radial electric field Er is derived from the ambipolar condition at the initial plasmastate. In #88343, three solutions of the ambipolar radial electric field are found in the radial
region 0.265 ≤ ρ ≤ 0.785. The profile of the radial electric field is dynamically fixed. The positive
Er is chosen. We study the simulation result of the Ti profile. The simulation results for theradial Ti profile show a good agreement with the experimental results. For modeling L in terms

of LTi , three runs at the diﬀerent values of LTi are needed at a radial point before the dynamicaltransport simulation. It takes about twenty minutes per one time of the program run. If the valueof L is calculated at each time step in the TASK3D, about one hundred thousand times of the
program run are necessary. Therefore, the transport simulation of an extremely low computational
cost can be achieved due to the additional modeling of the turbulent ion heat diﬀusivity basedon the gyrokinetic simulation. In #109081, one negative radial electric field is shown from theambipolar condition and is dynamically fixed. At the initial state, the ITG mode is stable in theregion 0.0 ≤ ρ ≤ 0.8. The simulation results using the diﬀusion equation for the ion temperatureshows the higher value than the experimental results in the core region. At the plasma center,the value of Ti of the simulation result is 4.2keV, which is higher than the value of Ti : 3.2keV forthe experimental results in LHD. The ITG mode becomes unstable in the region 0.18 ≤ ρ ≤ 0.42
at the stationary state of the simulation results. However, the neoclassical transport is dominantcompared with the turbulent transport driven by the ITG instability. In the case of #109081, thedensity profile is slightly peaky. The value of Ln becomes much smaller than the case in #88343,
where the density profile is flattened. At ρ = 0.5, the values of ηi (= Ln /LTi ) are −12.0 and 1.37 for#88343 and #109081 in LHD, respectively. In this study, we perform the gyrokinetic simulation
with the adiabatic electrons. It is not enough for showing the Ti profile close to the experimental
one only by the ITG stability analysis using the adiabatic electrons, especially for #109081 case.
The gyrokinetic simulation with the kinetic electrons shows the larger ion heat flux than that bysolving the gyrokinetic equation with the adiabatic electrons [11].
To reduce the simulation cost, how to apply the transport model based on a gyrokinetic simulation
to a dynamical transport code is shown within the accuracy of the gyrokinetic simulation. The
GKV-X code is performed for the additional modeling in terms of the ion temperature gradientscale length instead of the linear growth rate term for the ITG mode. The formula of the zonal
flow decay time in the reduced model is derived at the given magnetic field structure. The simulation results by applying the additional modeled turbulent ion heat diﬀusivity to the dynamicaltransport code are shown. The promising technique of a low computational cost for adapting the
gyrokinetic turbulent transport model to a transport code is proposed within allowable errors. The
dynamical transport simulation of the #109081 in LHD clarifies the importance of the gyrokinetic
simulation with the kinetic electrons. The construction of the reduced model for the heat diﬀusivity is needed by the electromagnetic gyrokinetic simulation with kinetic electrons. This is left for
future study.
The authors are grateful to the LHD Experimental Group for making the experimental data available. This work was partly supported by JSPS KAKENHI Grant Number 23561002, the NIFS
Collaboration Research program (Plasma Simulator), NIFS14KNST062 and the Collaborative Research Program of Research Institute for Applied Mechanics, Kyushu University, 27FP-11.
[1] Garbet X, Idomura Y, Villard L and Watanabe T -H 2010 Nucl. Fusion 50 0433002[2] Nunami M, Watanabe T -H, Sugama H and Tanaka K 2012 Phys. Plasmas 19 042504[3] Nunami M, Watanabe T -H and Sugama H 2013 Phys. Plasmas 20 092307[4] Sugama H and Watanabe T -H 2006 Phys. Plasmas 13 012501[5] Ferrando-Margalet S, Sugama H and Watanabe T -H 2007 Phys. Plasmas 14 122505[6] Yokoyama M, et al., 2012 Plasma Fusion Res. 7 2403011[7] Toda S, et al., 2014 Journal of Physics: Conference Series 561 012020[8] Tanaka K, et al. 2010 Plasma Fusion Res. 5 S2053[9] Murakami S, et al., 2015 Plasma Phys. Control. Fusion 57 054009[10] Nunami M, Watanabe T -H, Sugama H and Tanaka K 2011 Plasma Fusion Res. 6 001[11] Ishizawa A, et al., 2015 Nucl. Fusion 55 043024

13th Burning Plasma Simulation Initiative (BPSI) Meeting

Sensitivity and prediction of plasma profile inturbulent transport of helical plasmas
M. Nunami, M. Nakata, S. Toda, and H. Sugama
National Institute for Fusion Science, Toki 509-5292, Japan
Sensitivity to plasma temperature proﬁle in turbulent transport of helical plasmas is discussed by the local ﬂux-tube gyrokinetic simulations. Temperature proﬁles are reproduced
from the experimental observation within experimental error bars which can be regarded as
the standard deviation, and the allowable ranges of the ion temperature gradient within theerrors can be estimated. According to the allowable ion temperature gradient, the permissible ranges of the ion heat diﬀusivity are expected. Furthermore, if the heat diﬀusivities are
ﬁxed with the ﬂux-matching method, the radial proﬁles of the temperature gradients can bepredicted.

1 Introduction
Turbulent transport has been considered to be one of the most critical issues in the magneticallyconﬁned plasmas researches. The gyrokinetic approaches are powerful to analyze the transportphenomena, which are considered to be driven by the drift-wave plasma turbulence. To design
fusion reactors, especially, it is necessary to predict the transport ﬂuxes, the plasma temperature and density proﬁles. The turbulent transport ﬂuxes are quite sensitive to the proﬁles ofthe plasma temperature and density. Here, sensitivities of the transport ﬂuxes to the plasma
proﬁles are discussed within the experimental error ranges of the temperature proﬁles in helical plasmas. The allowable range of the temperature gradients is obtained from re-generated
radial functions for ion temperature proﬁles, which is obtained by random sampling from theredistributed temperature data within experimental error bars which can be regarded as thestandard deviations of the distributions.

2 Temperature gradient proﬁles within experimental errors
Radial proﬁles of plasma temperatures observed in the LHD experiments have the measurement
error bars at each radial position. The error bars are regarded as the width of the standard
deviation of the measurement data. Therefore, we can reproduce the temperature data by thenormal distributions with the standard deviations which correspond to the experimental error
author’s e-mail: nunami.masanori@nifs.ac.jp

13th Burning Plasma Simulation Initiative (BPSI) Meetingbars. If the normal distributions for the plasma temperatures are reproduced at each radial
position, we can give the ion temperature proﬁles by ﬁtting the randomly sampling points fromthe reproduced data with the function,
Ti (ρ) =k∑an ψ n (ρ) =

n=0

k
∑

bn ρ2n ,(1)

n=0

where the labeling index of the ﬂux surfaces, ρ ≡√ψ/ψa , is a dimensionless quantity. Hereψ represents the toroidal magnetic ﬂux, ψa is deﬁned at the last closed surface, an and bn are
ﬁtting coeﬃcients. If the ﬁtting functions for the ion temperature are obtained, the proﬁle of theradial gradient of the temperature with a certain allowable ranges according to the experimental
errors as shown in Fig. 1.

R0 / LTi

20.0

E 
P  
0%

+2

10.0

0.0
0.0

0.2

0.4

0.6

0.8

1.0

ρ
Fig. 1 Radial proﬁle of the gradients of the ion temperature with the allowable ranges according to theexperimental errors (red solid curves), and the predictions of the temperature gradients by the
ﬂux-matching method.

3 Permissible range of the transport coeﬃcients
We perform the gyrokinetic analyses for the ITG turbulent transport by GKV codes[1, 2, 3]
within the allowable range of the temperature gradients obtained in previous section. Figure 2
shows the ion heat diﬀusivities obtained by the ITG turbulence simulations and the LHD experiment. In outer radial region except for ρ > 0.8, the simulations within the experimental errorsof the ion temperature can cover the experimental diﬀusivities, since the experimental errors enhance the ranges of the simulation results. In the core region, on the other hand, the simulation
results have larger ambiguity rather than outer region, due to the fact that the experimental
errors are larger in the core region. Up-shifts of the critical gradients from linear critical values
are also observed and critical gradients go away from experimental values for more outer radial13th Burning Plasma Simulation Initiative (BPSI) Meeting
region. In the outer radial region, there exists large density gradients which may cause the TEM
instability, while the density proﬁle in core region is almost ﬂat.
On the other hand, if the heat diﬀusivities are ﬁxed to match the experimental observationsof the transport ﬂuxes, the radial proﬁles of the temperature gradients can be predicted byso-called ﬂux-matching method[4]. In ﬁg. 1, we show the predictions of the ion temperature
gradients, where the predictions agreed with the experimentally allowable ranges of the temperature gradients. Using the method, we may predict not only the temperature but also thedensity proﬁles against the error bars of the experimental density proﬁle including the diﬀerent
instability sources such as the TEM.

2

χi (m /s)

6.0

4.0

2.0

0.0
0.0

Experiment
Anomalous
Simulation w/ nominal grad.
Permissible range

0.2

0.4

0.6

0.8

1.0

ρ
Fig. 2 Ion heat diﬀusivities obtained by the experiment (solid curve) and the gyrokinetic simulations(blue symbols) with nominal temperature gradients. The simulations within the error bars oftemperature proﬁle enhance the ranges (red error bars).

Acknowledgment
This work is supported in part by the RIAM collaborative Research Program (27FP-12), theJapanese Ministry of Education, Culture, Sports, Science, and Technology, Grant (26820398),
NIFS collaborative Research Programs (KNTT031 and KNST080), and by use of Helios system
at International Fusion Energy Research Center (Project code: VLDGK ST and GKTBMDL).

[1] T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).
[2] M. Nunami, T.-H. Watanabe, and H. Sugama, Plasma Fusion Res. 5, 016 (2010).
[3] A. Ishizawa, et al., 25th IAEA Fusion Energy Conference, TH/P6-40 (2014).
[4] T. Görler, et al., Phys. Plasmas 21, 122307 (2014).

Multi-species ITG-TEM turbulence simulations in ITER
D-T-He plasma
Motoki NAKATA1 , Mitsuru HONDA2 , and Masanori NUNAMI1
National Institute for Fusion Science1 , Japan Atomic Energy Agency2
Multi-species ITG-TEM driven turbulent transport of D-T ions and He-ash in the burning
plasmas is explored by electromagnetic gyrokinetic Vlasov simulations with realistic ITERshape equilibria. Gyrokinetic-simulation-based quantitative evaluation of a steady burning
condition is realized by comprehensive nonlinear scans, and the strong impacts of the D-T
fuel ratio and He-ash on turbulent energy and particle fluxes are clarified. New findings inthis study, which are crucial for the burning plasma performance, are (i) imbalanced D-T
particle transport strongly influenced by He-ash accumulation, and (ii) identification of a
profile regime in which the burning condition with He-ash exhaust and fuel inward-pinch issatisfied.

1

Introduction

The first-principle based gyrokinetic simulation is widely recognized as a promising way to
predict turbulent transport in fusion plasmas, and nowadays the validation studies against
existing experiments become more active[1–3], where reasonable agreement in the ion andelectron heat fluxes between gyrokinetic simulations and experimental measurements has
been confirmed. Burning plasmas in ITER and DEMO are always composed of multiple
ion species such as fuel isotopes(D and T) and He-ash produced by the fusion reaction, andmore complex turbulent transport processes are expected in comparison to the single-ion
plasmas. Since simultaneous measurements of the kinetic profiles for all species are limited
even in experiments, systematic studies on the particle and heat transport are indispensable
to predict the confinement performance and to optimize the impurity exhausts and D/Tfueling.
In this study, the ITG-TEM driven turbulent transport in realistic ITER D-T-He plasmas is investigated by a gyrokinetic Vlasov simulation GKV[4,5] with multiple ions andfully kinetic electrons including their inter-species collisions[6,7], where a good prediction
capability has been confirmed against the actual JT-60U tokamak L-mode experiment[3].
In addition to the D-T fuel ratio dependence of the turbulent transport levels, we explore,
for the first time, an optimum profile regime, which satisfies the steady burning condition[8]with eﬃcient He-ash exhaust and fuel inward-pinch, from the gyrokinetic-simulation pointof view.

2

Impacts of D-T ratio and He-ash on turbulent transport

Nonlinear gyrokinetic Vlasov simulations with D, T, He, and real-mass kinetic electrons
are performed for an ITG-TEM unstable plasma with realistic ITER-shape[5], where thenormalized temperature and density gradient lengths at ρ = 0.5 are set as Rax /LTi = 6,
Rax /LTe = 8, and Rax /Ln = 2, respectively. The GKV simulations reveal diﬀerent saturationlevels and spatial structures of the turbulent fluctuations in D-T ions and He-ash [cf. Fig.
1(c)], and each transport level is evaluated by the time average in the quasi-steady phasesof the ITG-TEM turbulence. Figures 1(a) and 1(b) show the impact of the D-T fuel rationT /(nT + nD ) on the turbulent energy and particle transport levels of D and T in thecase with and without He-ash, where the ion fluxes are normalized by the electron one.
One finds that the D-T energy fluxes QD,T /Qe are reduced by the He-ash accumulation,
but the overall tendency of the linear D-T ratio dependence is unchanged [Fig. 1(a)], e.g.,
QD ≃ QT at the almost 50%-50% D-T ratio, regardless of the He-ash accumulation. On theother hand, in the particle fluxes ΓD,T /Γe , imbalanced D-T transport due to the nonlinear
D-T ratio dependence is revealed [Fig. 1(b)], where ΓD ̸= ΓT even at 50%-50% D-T ratio.
Moreover, the He-ash accumulation leads to not only the reduction in the particle flux levels,
but also stronger imbalance of the D-T particle transport.

3

GK-simulation-based evaluation of burning condition

In addition to the D-T ratio dependence of the turbulent transport, the present multiion species gyrokinetic simulations with He-ash can provide the quantitative evaluationof the steady burning condition which is never addressed by the single-ion approaches.
For sustainable fusion reactions with the eﬃcient He-ash exhaust(outward flux) and fuelpinch(inward flux), a steady burning condition is given by ΓHe > 0, ΓD,T < 0, and ηi Ti ΓHe >Qi /α∗ (or equivalently τHe < α∗ τE ), where α∗ ∼ 7 − 15 is a constant depending on thedivertor and wall conditions[8]. To this end, the density gradient(Rax /Ln ) scans of thesteady turbulent energy and particle transport levels are performed. As shown in Fig. 2(a),it is clearly identified that there exists a profile regime of Rax /Ln ≤ 1.27(hatched region)in which the above steady burning condition with α∗ = 7 is satisfied. The spectral analysisof the turbulent particle fluxes of He-ash and Tritium reveals that, in contrast to the weak
Rax /Ln dependence of the outward He-ash flux, the components of the inward Tritium flux
increase for the lower Rax /Ln , indicating significant oﬀ-diagonal eﬀects (∝ Rax /LTi ) on theparticle transport.

4

Summary

In this study, multi-species ITG-TEM driven turbulent transport of D-T ions and Heash in the burning plasmas is explored by electromagnetic gyrokinetic Vlasov simulations
with realistic ITER-shape equilibria. Gyrokinetic-simulation-based quantitative evaluationof a steady burning condition is realized by comprehensive nonlinear scans, and the strongimpacts of the D-T fuel ratio and He-ash on turbulent energy and particle fluxes are clarified.
These findings on the He-ash impact and on the profile dependence for the steady burning
condition can contribute to deeper understanding of the multi-species transport processes
in tokamak and helical burning plasmas, and encourage the quantitative predictions andoptimization for forthcoming ITER and DEMO.

Figure 1:

D-T ratio dependencies of (a) energy and (b) particle fluxes for D and T normalized
by the electron flux, where the cases with 0%(solid lines) and 10%(dashed lines) He-ash are plotted.
Imbalanced turbulent fluxes significantly appear in the D-T particle transport, rather than theenergy transport. (c)Turbulent fluctuations of D(blue), T(yellow), and He(red) in the multi-species
ITG-TEM turbulence.

Figure 2: (a) Density-gradient scan of turbulent particle and ion energy fluxes, where the steadyburning condition with α∗ = 7 is satisfied in Rax /Ln ≤ 1.27(hatched region). (b) Wavenumber
spectra of Tritium and He-ash fluxes at several density gradients, where the inward-flux components
appear in Tritium.

Reference[1] T. L. Rhodes, C. Holland et al., Nucl. Fusion 51 063022 (2011)[2] T. Golrer, A. E. White et al., Phys. Plasmas 21 122307 (2014)[3] M. Nakata, M. Honda et al., Proc. on the IAEA-FEC2014, THP7-38 (2014)[4] T. -H. Watanabe and H. Sugama, Nucl. Fusion 46 24 (2015)[5] M. Nakata, A. Matsuyama et al., Plasma Fusion Res. 9 1403029 (2014)[6] M. Nakata, M. Nunami et al., Comput. Phys. Commun., 197 61 (2015)[7] M. Nunami, M. Nakata et al., Plasma Fusion Res. 10 1403058 (2015)[8] D. Reiter, G. H. Wolf et al., Nucl. Fusion 30 2141 (1990)

トロイダルプラズマにおける乱流と計測に関するシミュレーション研究
Simulation study of turbulence and diagnostic in toroidal plasmas
N. Kasuya, K. Kawadu 1, M. Nunami 2 and M. Yagi 3
Research Institute for Applied Mechanics, Kyushu University, Kasuga, Fukuoka 816-8580,
Japan
1
Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, Kasuga,
Fukuoka 816-8580, Japan
2
National Institute for Fusion Science, Toki, Gifu 509-5292, Japan
3
Japan Atomic Energy Agency, Obuchi, Rokkasho-mura, Aomori 039-3212, Japan
E-mail: kasuya@riam.kyushu-u.ac.jp
Resent turbulence simulations in three-dimensional (3-D) magnetically confined plasmas
show the various aspects of plasma turbulence [1]. We have been carrying out numerical
diagnostics using 3-D simulation data of helical plasmas. Here we present the turbulence
analysis results (i) from simplified geometry for detailed nonlinear mechanism of heattransport, and (ii) from real 3-D geometry for comparison with experimental observations.
Competition of various temporal and spatial scales involved in the model, turbulent structures
revealed in the transient response, and their features through the observation are shown fromthe numerical diagnostics.
Many turbulence codes have been developed for the calculation in experimental magnetic
configurations. For the helical plasmas, using the GKV-X code, which is the gyrokinetic
Vlasov code with 3-D equilibrium, quantitative comparisons with experiments have been
carried out [2]. On the other hand, a reduced fluid code with simplified magnetic
configuration [3] is also useful for the fundamental understandings. It is important to integrate
the results from the various approaches. One of the ways for the integration is by analyzing
the various simulation data with the routines to give same diagnostics. Turbulence Diagnostic
Simulator (TDS) [4] is a combination of turbulence codes, measurement modules and analysis
routines, to carry out numerical experiments of plasma turbulence, which can be utilized as
the platform of the data analyses. The turbulence features revealed from the TDS are
explained in this paper.
First, the fundamental mechanism with simplified configuration, which has a circularpoloidal-cross-section, is described. Global nonlinear simulations of drift-interchange modes
in helical plasmas are carried out using a reduced MHD model [3]. The averaging with thehelical pitch is adopted, so the configuration can be treated as a quasi-2-D one with the helicalmagnetic curvature. The model is used as the fundamental one including the global andlocalized modes, nonlinear couplings with the Reynolds stress, and collisional transport
processes. The nonlinear saturated states are sustained with an applied heat source, and themode structures broadened in the radial direction and localized near their rational surfaces areboth excited. These fluctuations have the self-correlation time of 30, though those of higher
m modes are smaller than 10, where  is the ion acoustic time and m is the poloidal mode
number. A non-local effect has been studied in dynamical transport phenomena by thesimulation of the heat source modulation [3]. The nonlinear process plays the key role in theresponse, which takes a finite temporal duration ~ 30 for the energy redistribution. By
conditional averaging, the characteristic response can be extracted. Figure 1 shows the 2-D
pattern of the heat flux to propagate the modulation from the core to the edge. A spiral streamline appears, rotating to the electron diamagnetic direction with the ballooning feature. This
structure is not a linear mode structure, but a nonlinearly formed one, and is sustained for~100 that is longer than the turbulence decorrelation time. Concerning the electromagnetic
aspect, the other time scale appears. There exists
an oscillation with the period of ~ 160. The
model equation reduced with the ordering in thesimulation gives an Alfvénic oscillation of the m= 1 modes, whose nonlinear couplings are foundto give the oscillation of the mean components.
In addition, flattening by the magnetic curvaturein the pressure evolution stabilizes the mode
strongly. In this way, the reduced model canreveal the competition of various time-scales,
Fig.1: Conditional averaged 2-D patternwhen the MHD modes with a rather longof heat flux in the poloidal cross-section
wave-length play the main role in the saturation.
of the torus plasm
Next, the case with the 3-D equilibrium and shorter wave-length perturbations is explained.
Flux-tube simulations of Ion-Temperature-Gradient modes in the LHD configuration havebeen carried out with the GKV-X code. Quantitative comparison of the thermal diffusivity
coefficient has shown the good agreement with experiments [2], and other turbulent
characteristics, such as the perturbation spectrum, should be compared for the validation. We
develop an analysis routine for the GKV-X data taking into account of the line of sight of theexperimental diagnostics. The code uses the magnetic surface coordinate, so theinterpolation is necessary for evaluating thephysical quantities in the real space. The spatial
resolutions are limited by those of thesimulation, r ~ 0.4 [mm] and  ~ 50 [mm].
For the 2-D spectrum analysis of density
perturbation, 2-D profiles are obtained (Fig. 2).
Efficiency of the method to resolve the local
spectrum is tested. Assuming k// = 0, the typicalk direction (ky / kx) can be set, corresponding to
the magnetic field direction. The magnetic field
direction changes with the variation of themagnetic field in the vertical direction, which
can use for the reconstruction in helical plasmas.
The extracted component is dominant in thecertain vertical position, but the contributions
Fig. 2: Poloidal cross-section of
from the other positions are not negligible, due
perturbations and its local patterns at
to non-zero k//, so can use the estimation of thedifferent vertical positions along the line
experimental error bars.
of sight.
In this way, the numerical diagnostics have been carried out in the 3-D torus magnetic
configurations. Variety of the numerical diagnostics from several view-points can give
physical understanding and quantitative comparison of turbulent plasmas.

Acknowledgements
Authors acknowledge discussions with Professor K. Itoh, Professor S.-I. Itoh, Professor S.
Inagaki, Professor H. Sugama and Dr. M. Sasaki. This work is supported by the Grant-in-Aid
for Young Scientists (24760703), for Scientific Research (23244113) of JSPS, by thecollaboration program of NIFS (NIFS15KNST082, NIFS15KNTS040, NIFS13KOCT001)
and of RIAM of Kyushu University. Some numerical simulations were carried out on the
Plasma Simulator at National Institute for Fusion Science and the Helios Super Computer
System at Computational Simulation Center of International Fusion Energy Research Center.
References[1] X. Garbet, et al., Nucl. Fusion 50 (2010) 043002.
[2] M. Nunami, et al., Phys. Plasmas 19 (2012) 042504.
[3] N. Kasuya, et al., Phys. Plasmas 21 (2014) 110701.
[4] N. Kasuya, et al., Plasma Fusion Res. 8 (2013) 2403070.

Development of the multi-fluid transport codeon the tokamak flux coordinates
M. Honda1 , A. Fukuyama2
1 Japan Atomic Energy Agency, Naka, Japan
2 Graduate School of Engineering, Kyoto University, Kyoto, Japan

Introduction
We have developed a one-dimensional fluid-type transport code TASK/TX [1] in the cylindrical coordinates. The code is essentially based on a self-consistent two-fluid model, which
consists of two-fluid equations (conservation of mass, momentum and energy) plus Maxwell’sequations. It also involves the equations for beam ions and neutrals. It diﬀers from conventional∑
diﬀusive transport codes mainly in that: the quasi-neutrality condition ne = i Zi ni need not be∑
imposed and the ambipolar flux condition a Γa = 0 is intrinsically satisfied; the flux-gradientrelationship is not used for particle transport; the neoclassical features are self-consistently reproduced. Here we derive the governing equations of TASK/TX on the axisymmetric flux coordinates (ρ, θ, ζ) and then numerically implement them.
Maxwell’s equations
On the axisymmetric flux coordinates, the magnetic field is expressed as B = ∇ζ × ∇ψ + I∇ζ,
where I(ψ) = RBt . Maxwell’s equations consist of Gauss’s law (Poisson’s equation), Faraday’sA, it is found that the covariant comlaw and Ampère’s law. From Faraday’s law E = −∇Φ − Ȧ
ponents of the electric fields are strongly tied to (the temporal change in) the magnetic fluxes:√Eθ = E · g∇ζ × ∇ρ = −ψ̇t , Eζ = gζζ E · ∇ζ = ψ̇(= REt ). Ampère’s law relates (the spatial gradientof) the magnetic flux to the current. Taking the scalar product of Ampère’s law with ∇ζ and thesubsequent flux-surface average (FSA) yield (cf. FSA Grad-Shafranov equation)[ ⟨2 ⟩ ∂ψ ]1∂⟨ jζ ⟩∂ψ1 ∂ψ̇′ |∇ρ|
=
V
−
µ
≡ ψ̇,
0 −2 ,
2
′
−2
2
∂ρ
∂t
c ∂t V ⟨R ⟩ ∂ρ
R
⟨R ⟩

(1)

the latter of which is not only the definition of ψ̇ but also one of the governing equations that
constitutes Maxwell’s equations. This is because it is just alternative form of Faraday’s law.
Taking the scalar product of Ampère’s law with B and then subtracting (1) give the equations
for the toroidal flux in the form:⟨]2⟩ ∂ [⟨B j∥ ⟩ − I⟨ jζ ⟩1 ∂ψ̇t1∂ψt′ |∇ρ|
=
V
+
µ
,
0
⟨Bθ ⟩
c2 ∂t
R2 ∂ρ V ′ ⟨R−2 ⟩ ∂ρ∂ψt
≡ ψ̇t .
∂t(2)

When neglecting the displacement current term that is negligibly small, we readily find that theright-hand sides of (1) and (2) reproduce the expressions of the toroidal and parallel currents.
Finally, the Coulomb gauge allows us to simply write FSA Gauss’s law as follows:
[
]
1 ∂
1∑
′
2 ∂Φ
V
⟨|∇ρ|
⟩
=
−
ea na .
V ′ ∂ρ
∂ρ
ϵ0 a

(3)

Continuity equations
A flux-surface-averaged continuity equation is simply given by
1 ∂ ′
1 ∂
(V na ) + ′ [V ′ na ⟨(uua − u g ) · ∇ρ⟩] = S na .
′V ∂t
ρ V ∂ρ(4)

TASK/TX is distinguished from other transport codes in that na ⟨(uua − u g ) · ∇ρ⟩ is not approximated by the convection-diﬀusion relationship, but is treated as a dependent variable: The gridvelocity u g and the flux na ⟨uua · ∇ρ⟩ are both self-consistently calculated in the system.
Momentum equations
The parallel motion of plasma species regulates neoclassical properties in axisymmetric systems. The lowest order of the viscous term is the neoclassical viscosity. The viscous term andthe friction term can be calculated according to the moment approach [2]. Therefore, defining⟨Bq̂a∥ ⟩ = 2⟨Bqa∥ ⟩/(5pa ), the following equations should be solved:
ma na

∑
∑
∂⟨Bua∥ ⟩abab= − µ̂a1 (⟨Bua∥ ⟩ − BV1a ) − µ̂a2 (⟨Bq̂a∥ ⟩ − BV2a ) +ℓ11⟨Bua∥ ⟩ −ℓ12⟨Bua∥ ⟩∂tbb+ ea na ⟨BE∥ ⟩,

(5)
∑

∑

∂⟨Bq̂a∥ ⟩
5
ab
ab
ma na
= − µ̂a2 (⟨Bua∥ ⟩ − BV1a ) − µ̂a3 (⟨Bq̂a∥ ⟩ − BV2a ) −ℓ21⟨Bua∥ ⟩ +ℓ22
⟨Bq̂b∥ ⟩. (6)2∂tbbThe toroidal momentum equation is important for TASK/TX because it includes the v · ∇ρterm that provokes a particle flux as well as a j × B torque once losses of beam ions produce
the non-ambipolar flux. In conventional transport codes the particle transport coeﬃcients evaluated by external modules are directly substituted in the particle transport equation, whereas
in TASK/TX coupling of the equations self-consistently brings about particle transport throughthe continuity equation. While the first-order viscous stress term vanishes due to the CGL form,
the second-order term can be expressed as a combination of a convective momentum flux, a.k.a.
inward pinch, and a diﬀusive one plus a residual stress. We have to add a turbulent force FaQL
that drives a turbulence-induced quasilinear particle flux. The toroidal momentum equation isfinally given by, where ⟨La ⟩ ≡ ma na ⟨Ruaζ ⟩,[∂⟨Ruaζ ⟩
1 ∂ ′
1 ∂ ′
ρ
ρ
2
(V
⟨L
⟩)
=
−
V
⟨|∇ρ|⟩v
⟨L
⟩
+
(u
−
u
)⟨L
⟩
−
⟨|∇ρ|
⟩χ
m
n
a
aζ
a
a
aζ
a
a
a
g
V ′ ∂tV ′ ∂ρ
∂ρ
] ∑
∑
∂ψ
ρ
ab
ab I
+ ⟨Πres⟨Bq̂b∥ ⟩ + ea na ⟨REζ ⟩ + ea na ua . (7)ℓ11⟨Rubζ ⟩ −
ℓ12
a ⟩ +
2
∂ρ
⟨B
⟩
b
b
The radial momentum equation or the radial force balance equation is essentially equivalent
to the first-order flow within the flux surface. The leading order in δ is O(1) for pressure andLorentz force terms, and the other terms are practically ineﬀective. Thus, we obtain0=−

1 ∂pa 2 2
ea ∂Φ
ea
ea
(⟨B ⟩⟨R ⟩ − I 2 ) − na (⟨B2 ⟩⟨R2 ⟩ − I 2 ) + na I⟨Bua∥ ⟩ − na ⟨B2 ⟩⟨Ruaζ ⟩. (8)
ma ∂ψ
ma ∂ψ
ma
ma

Numerical results
For the sake of simplicity, the large aspect ratio plasma with circular cross section is assumed.
JT-60U-like plasma parameters are given in the following simulation. Figure 1 shows the radial
profiles of (left) the diﬀerence in the charge density and the resultant Er , indicating that automatically the quasi-neutrality is well satisfied and the ambipolar flux condition is satisfied aswell. A break of the local charge neutrality in the core region, (ne − ni )/ne , is of the order of
10−7 , which is always neglected in conventional transport codes, but this tiny discrepancy is
indispensable to give rise to the radial electric field Er through Eq. (3). The right figure shows
the comparison of the parallel current profiles estimated by TASK/TX and three neoclassicaltransport model/modules. The breakdown of the current in TASK/TX cannot be known, but wefind that the parallel current calculated by TASK/TX well agrees with that by other models.
This work was supported by JSPS KAKENHI (No. 25820442).
References[1] M. Honda and A. Fukuyama, J. Comput. Phys. 227, 2808 (2008).
[2] S. P. Hirshman and D. J. Sigmar, Nucl. Fusion 21, (1981) 1079.
1.5

Er
Neutrality

1.0

Er [kV/m]

4
0.5

2
0

0

-2

-0.5

-4
-6
-8
0

0.2

0.4

ρ

0.6

0.8

1

TASK/TXNCLASS
SauterMI

5j// [MA/m2]

6

Neutrality [1013Cm-3]

8

4
3
2

-1.0

1

-1.5

0
0.0

0.2

0.4

ρ

0.6

0.8

1.0

Figure 1: Radial profiles of (left) Er and quasi-neutrality and (right) the parallel current.

Bifurcation of the interchange mode growth rate androtation frequency due to the perpendicular heat
conductivity
T. Nicolas1 , K. Ichiguchi1,2 , Y. Todo1,2 , M. Sato1 , Y. Suzuki1,2 , and
A. Ishizawa1,2
1National Institute for Fusion Science, 322-6 Oroshi-cho, Toki 509-5292, JapanThe Graduate University of Advanced Study, Sokendai, Toki 509-5292, Japan2

The linear stability against interchange modes of heliotrons like the Large Helical Device has not been fully understood yet. The interchange mode is the dominant instability
in nearly current-free 3D configurations. Recently, experiments in the so-called inwardshifted configuration have shown the machine can be operated safely up to hβi ∼ 5%
without major MHD event [1]. Here hβi is the volume averaged ratio between kinetic andmagnetic pressure. The rotation of the mode, which seems to play a role in the stability [2],is also not understood. Namely, the modes are ob×10served to rotate in the electron direction [3], whereas
3
Bifurcation
extended MHD theory without heat conductivity
but including electron and ion diamagnetic effects
2
predicts rotation in the ion direction for Mercier
1
unstable modes and almost no rotation for Mercier
stable modes [4].
0
In this contribution, we present a new effect on−1
Growth Ratethe linear Mercier unstable interchange growth rate
Frequency
χccaused by the perpendicular heat conductivity χ⊥ .
−2
10
10
When the heat conductivity is introduced in a 3χ⊥field reduced model of Strauss type equations [5], Figure 1: Frequency and growth
it is found that a second mode with an interchange rate of the first two most unstable
parity and no sign change in the potential wave func- modes
tion appears, in addition to the usual most unstableinterchange mode. The second mode is more stable, so its growth rate is lower. As can beseen in Fig 1, for small values of the heat conductivity, the two modes have zero frequency
and converge toward each other when χ⊥ is increased. When χ⊥ is equal to a critical valueχc , the growth rates of the two most unstable modes become equal. For χ⊥ > χc , the twomodes have same growth rate and opposite frequency, as shown in Fig. 1. Therefore there
is a degeneracy, since two modes have the same growth rate, but opposite frequencies. If
the plasma is in this regime, then the observation of the linear phase of the interchange
mode corresponds to standing waves. The direction of rotation of the mode cannot be
identified.
This degeneracy is removed by the diamagnetic effects. In the domain χ⊥ > χc , thediamagnetic flows have a very different effect compared to what happens for ideal modes.
For ideal modes, the diamagnetic effects are such that the growth rate and frequency ofGrowth rate, Frequency−2−5

Corresponding author: Timothee Nicolas timothee.nicolas@nifs.ac.jp

−4

−2

4.0 ×10a)3.5

ωi⋆ = 0γ1,1 (Positive branch), ωi⋆ = 1.8 × 10−2

γ1,2 (Negative branch), ωi⋆ = 1.8 × 10−2

Growth rate γ

3.0
2.5
2.0
1.5
1.0
0.5
0.0
0.06

b)

Frequency ωr

0.04
0.02
0
-0.02
-0.04
-0.06
-0.08 −610ωi⋆ = 0
ω1,1 (Positive branch), ωi⋆ = 1.8 × 10−2

ω1,2 (Negative branch), ωi⋆ = 1.8 × 10−2

10−5

10−4

Heat conductivity χ⊥

Figure 2: Frequency and growth rate of the first two most unstable modes with diamagnetic effects ωI? = 1.8 × 10−2 and no viscosity
the interchange mode are given [6] byr?
 γ = γI 1 − ωi
2γI ,
1

 ω = ω?
2 i

(1)where ωi? and γI are the ion diamagnetic frequency and the ideal growth rate respectively.
This law is approximately valid for small values of the heat conductivity, resistivity andviscosity. However, in the domain χ⊥ > χc , this law is not valid anymore. There are
two cases. If there is no viscosity, then the most unstable mode is destabilized and thesecond most unstable mode is stabilized. The frequency of the most unstable mode isin the ion direction, so the mode rotates in the ion direction. This situation is shown in
Fig. 2. In the second case, if there is a viscosity large enough to affect the ideal growthrate, then the most unstable mode is stabilized whereas the most stable mode, rotating in
the electron direction, is destabilized. As a result, the mode rotation, which is still given

−2Growth rate γ

3.0 ×10a)2.5

ωi⋆ = 0γ1,1 (Positive branch), ωi⋆ = 1.8 × 10−2

γ1,2 (Negative branch), ωi⋆ = 1.8 × 10−2

2.0
1.5
1.0
0.5

0.0
0.06 b)Frequency ωr

0.04
0.02
0
-0.02
-0.04

ωi⋆ = 0-0.06

ω1,2 (Negative branch), ωi⋆ = 1.8 × 10−2

10−6

ω1,1 (Positive branch), ωi⋆ = 1.8 × 10−2

10−5

10−4

Heat conductivity χ⊥

10−3

Figure 3: Frequency and growth rate of the first two most unstable modes with diamagnetic effects ωI? = 1.8 × 10−2 and viscosity ν = 10−5 .
by the dominant mode, is in the electron direction, as shown in Fig. 3.
This effect, which sheds a new light on the stability of the interchange mode in stellarators, is also found analytically and confirmed with 3D non-linear MHD simulations
utilizing the MIPS code [8, 9]. In the MIPS simulations, the linear phase of an ideallyunstable interchange mode is simulated in presence of some resistivity, viscosity and heat
conductivity. For large values of β, the relative effect of dissipation is small and a result
similar to the ideal case is found. For smaller values of β, closer to marginal stability but
still in the ideally unstable domain, it is found that the diamagnetic effects are destabilizing for small values of ωi? , and that the mode rotates in the electron direction, rather
than in the ion direction.
These results are detailed in references [8] and [9].

References[1] H. Yamada et al. Nucl. Fusion 51(9), 4021 (2011).
[2] Y. Takemura, et al. Nucl. Fusion 52(10), 102001 (2012).
[3] Y. Takemura, et al. Plasma and Fusion Research 8(0), 1402123 (2013).
[4] S. Gupta, et al. Phys. Plasmas 9(8), 3395 (2002).
[5] H. Strauss. Plasma Physics 22(7), 733 (1980).
[6] J. Connor, et al. Nucl. Fusion 24(8), 1023 (1984).
[7] Y. Todo, et al. Plasma and Fusion Research 5, S2062 (2010).
[8] T. Nicolas et al. Nucl. Fusion 56(2), 026008 (2015).
[9] T. Nicolas, et al. Plasma and Fusion Research (10) (2015).

Fokker-Planck Simulation for Runaway Electron Generation
including the Hot-tail Eﬀect
H. NUGA1 , A. MATSUYAMA2 , M. YAGI2 , A. FUKUYAMA1
1

Department of Nuclear Engineering, Kyoto University
2
Japan Atomic Energy Agency

Introduction In large size tokamak, such as ITER, disruptions have a potential to generate therunaway electron current which can reach to a few mega-amperes. Because in ITER, mitigateddisruptions by means of the massive gas injection (MGI) tends to yield the shorter thermal quench
time than in unmitigated disruptions, so-called hot-tail eﬀect[1, 2] might enhance the “primary”
runaway electron (RE) generation rate, which requires full kinetic calculation of the RE generation
process.
The aim of this paper is to compare the RE generation, especially RE current density profile,including and excluding the hot-tail eﬀect in ITER. For this purpose, the Fokker-Planck code,
named TASK/FP[3], is applied to simulations of the evolution of the relativistic electron momentum
distribution function, the RE current, and the induced electric field in a self-consistent manner.
The code treats the electron distribution function in two-dimensional momentum space for subrelativistic region (less than 0.5 MeV). In the present calculation, though the treatment of relativistic
REs is simplified, the thermal and sub-relativistic electron distribution function is analyzed with agood resolution, being useful for detailed calculation of the primary RE generation rate.
The RE generation in a high-current tokamak like ITER is expected to be dominated by the“secondary” generation. Nevertheless, our simulation shows that the hot-tail formation plays animportant role in determining the RE current profile after the current quench. Since even a small
increment of the seed current density is amplified by the avalanche eﬀect, the secondary RE current
density profile is modified by the hot-tail eﬀect.
Models In our present research, we have developed a Fokker-Planck simulation code TASK/FP
in order to calculate the time evolution of the relativistic momentum distribution function and thatof induced toroidal electric field self-consistently. It solves the Fokker-Planck equation:)(↔
∂f= −∇ · SC,E = −∇ · D C · ∇ f + FC,E f,(1)∂twhere ∇ is a derivative operator in momentum space and f is a relativistic electron momentum
distribution function. In this paper, the weak relativistic isotropic background collision term[4]↔

with relativistic Maxwellian is adopted for D C and FC .
The induced electric field E obeys following equations:∫dnrp=∇ · SC,E dp
(2)
dt
(
)−1/2
√
4π(Ze f f + 1)2
Ē − 1dnrs
πφ= nr
1 − Ē +≡ S s (nr , Ē), (3)dtτr ln Λ 3(Ze f f + 5)
3φ(Ze f f + 5)(Ē 2 + 4/φ2 − 1)1

()
1 ∂ ∂E∂r= µ0 jr ∂r ∂r∂t
j = σ sp + ecnr

(4)
(5)

where σ sp is a Spitzer conductivity, nrp and nrs are the primary and secondary RE densities, Ē ≡E/EC , and the critical electric field EC = ne q3e ln Λ/4πϵ02 me c2 , respectively. The total RE density
nr is evaluated from the sum of the primary (nrp ) and secondary (nrs ) REs. The generation rate of
the primary RE, which is generated by Dreicer mechanism, is defined in eq. (2) as electrons go
out from the momentum calculation domain: 0 < p < pmax and p2max /m ∼ 1MeV. On the other
hand, the generation rate of the secondary RE, which is generated by avalanche mechanism[5], isexpressed as the function of nr and Ē = E/EC in eq. (3), where φ = 1 − 1.46ε1/2 + 1.72ε, ε = r/Rand τr = 4πϵ02 m2 c3 /ne q4 ln Λ. Equation (4) expresses the electric field diﬀusion and this equation
is closed by the Ohm’s low (eq. (5)). Here, it is assumed that the current density consists of two
components: the bulk current , σ sp E, and the RE current, ecnr . Here all of REs are assumed to↔

travel with the velocity of light. The conductivity, σ sp , and the collision coeﬃcients ,D C and FC ,require the background plasma temperature, therefore, the evolution of the plasma temperature isgiven as a function of time. Since the behavior of the plasma during thermal quench is too complex
to describe accurately in RE generation time scale, here, the decay of the plasma temperature ismodeled with the time constant of the thermal quench, τq :()
T (t, ρ)=(T 0 (ρ) − T f (ρ)) exp −t/τq + T f (ρ),(6)
T 0 (ρ)=(T 0 (0) − T 0 (1))(1 − ρ2 )2 + T 0 (1),T f (ρ)=T f (0)(1 − 0.9ρ2 ),

(7)
(8)

where ρ denotes the normalized minor radius, T 0 and T f are the initial and post-quench temperature,and τq is the time constant of the thermal quench time, respectively. Electrons and ions are assumed
to have the same temperature.

current [MA]Results In order to investigate the eﬀect of the hot-tail, comparisons have been made for the
results including and excluding the hot-tail eﬀect. The results excluding the eﬀect are obtained by
the simulation which use the theoretical primary RE generation rate model[6] instead of eq. (2).
Under the condition of the steady electric field and constant plasma temperature, the values of the
primary RE generation rate calculated by these two models have good agreement.
We employed the ITER like plasma parameters, such as major radius R = 6.2m, minor radiusτs (2vth)=1.15
τs (3vth)= 3.88
a = 2.0m, initial plasma current I p = 15MA, ini14
12
total RE current
tial temperature T 0 (0) = 22.7keV, T 0 (1) = 2.27keV
10
post-quench temperature T f (0) = 10eV, and electron
8
density ne (0) = 1 × 1020 m−3 . The eﬀective charge is
Solids : including the hot-tail6
Dotted: excluding the hot-tail4
assumed to be Ze f f = 1.
primary RE current
2
Figure 1 shows the relation of RE current includ0
1
2
3
4
5
τq [ms]ing and excluding the hot-tail eﬀect to the thermal
quench time τq . Red and green curves denote thetotal and primary RE currents and solid and dotted Figure 1: Generated primary and total RE
currents against the thermal quench time.
2

2

current dens. [MA/m ]2current dens. [MA/m ]curves denote RE currents including and excluding the hot-tail eﬀect. Dashed vertical lines show
2 2 3
4
the e-e slowing down time, τees (v) = 2πϵ0 m v /ne q ln Λ, calculated with the velocity 2vth and 3vth ,where vth is the thermal velocity for the pre-quench temperature. From this figure, it is found that
if the thermal quench time is shorter than τees (2 − 3vth ), the hot-tail eﬀect enhances the total RE
generation.
Figure 2 shows the primary and secondary RE
current density profiles for τq = 1.5ms and 4.0ms.
(a) 6For τq = 4.0ms, which is slower than τees (3vth ), proτq =1.5ms
5files between including and excluding the hot-taildotted: excluding hot-tail4
eﬀect have good agreement. In the contrary, forsolid: including hot-tail3τq = 1.5ms, the hot-tail eﬀect changes RE current
primary2
density profiles. The eﬀect enhances the primary RE
secondary1current density in the whole region. On the other0
hand, the secondary RE current density decreases in
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1ρ
inner region (ρ < 0.4) and increases in outer region(b) 7(ρ > 0.4). The decrease of the secondary RE in innerτq =4.0ms
6
region comes from the reason that the increment of5
the primary RE current in inner region reduces the4secondary3
electric field. The increase of the secondary RE in2
outer region comes from the fact that the seed REs
primary
1in the region increase slightly owing to the hot-tail0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
eﬀect. The increment of the seed RE does not aﬀect
ρ
the electric field because it is negligible in the viewof the poloidal flux. As a consequence, the hot-tail Figure 2: Primary and secondary RE current
eﬀect make the RE current density profile broader density profile for (a) τq = 1.5ms and (b)even in the secondary RE dominant case.
4.0ms.
Summary The present work has investigated thecontribution of the hot-tail eﬀect to the RE generation in ITER. Despite the fact that the RE generation is dominated by the secondary current in ahigh-current tokamak like ITER, our results shows that the RE current density profile is modified
by the hot-tail eﬀect. Since the small increment of the primary RE current is amplified by theavalanche eﬀect, the final RE profile changes associated with the hot-tail generation.

[1] H. M. Smith and E. Verwichte, Phys. Plasmas, 15, 072502, (2008).
[2] T. Fehér, et. al., Plasma Phys. Control. Fusion, 53, 035014, (2011).
[3] H. Nuga, A. Matsuyama, et at., Plasma and Fusion Research, 10, 1203006, (2015).
[4] C. F. F. Karney, Computer Physics Report, 4, 183, (1986).
[5] M. N. Rosenbluth and S. V. Putvinski, Nuclear Fusion, 37, 1355, (1997).
[6] J. W. Conner and R. J. Hastie, Nuclear Fusion, 15, 415, (1975).

3

統合コード TASK による炉心プラズマ立ち上げシミュレーション

Burning plasma start-up simulation by the TASK code
永井一輝，福山淳，奴賀秀男
Kazuki Nagai, Atsushi Fukuyama, Hideo Nuga
京都大学大学院工学研究科
Graduate School of Engineering, Kyoto University
The start up of burning plasma has to be carried out within a limited external heating
power. Therefore, it is important to analyze the behavior of fast fuel ions quantitatively andself-consistently. To simulate the start up of ITER burning plasma, a tool to analyze thetime evolution of momentum distribution functions of fuel ions has been developed by using
the core plasma integrated code TASK. In this report, we introduce the Legendre expansion
in calculating the fusion reaction rate in order to reduce the computation time.
1 Introduction
In order to describe the behavior of burning plasmas, it is required to calculate the enhancement of the fusion reaction rate by the fast fuel ions produced by neutral beam injection orion-cyclotron wave heating. High resolution in momentum space is required for correctly
analyzing the wide energy range from 1keV to 3.5MeV. In order to improve the momentum
resolution and the computational performance, TASK/FP is parallelized with respect to theminor radius, the momentum and the particle species. Since the quadruple integral overpa , pb , θa , θb in the calculation of the fusion reaction rate requires a long computation timeand a large memory, the reduction of necessary computational resource is required.
2 Calculation of fusion reaction rate
In order to calculate the fusion reaction rate, it is necesarry to integrate over the momentum
space. For species a and b, we employ the momentum pa and pb , the magnitude of momentum
pa and pb , the pitch angle θa and θb , the momentum distribution function fa (pa , θa ) andfb (pb , θb ), the fusion reaction cross section σT and the particle relative velocity v 0 . The fusion
reaction rate R is calculated as∫∫R=fa (pa )fb (pb )σT v 0 dpa dpb
∫∫∫∫
2
= 4π
fa (pa , θa )fb (pb , θb )hσT v 0 i(pa , pb , θa , θb )p2a p2b sin θa sin θb dpa dpb dθa dθb , (1)where hσT v 0 i(pa , pb , θa , θb ) is σT v 0 averaged over the gyro phase φa and φb . The calculation of(1) takes a long computation time due to the quadruple integration. In order to reduce thecomputation time, the momentum distribution function and hσT v 0 i were expanded using the
Legendre polynomials Pl (cos θ) such asf (p, θ) =∞∑fl (p)Pl (cos θ)

(2)

l=0

and
hσT (E)v 0 i(pa , pb , θa , θb ) =

∞ ∑
∞
∑

hσT (E)v 0 il0 k0 (pa , pb )Pl0 (cos θa )Pk0 (cos θb )l0 =0 k0 =0

1(3)

wherefl (p) =2l + 12∫
f (p, θ)Pl (cos θ) sin θdθ

(4)

and
hσT (E)v 0 il0 k0 (pa , pb )∫∫2l0 + 1 2k 0 + 1=
hσT (E)v 0 i(pa , pb , θa , θb )Pl0 (cos θa )Pk0 (cos θb ) sin θa sin θb dθa dθb . (5)22Using the orthogonality of the Legendre polynomials, the final form of fusion reaction rate isexpressed as
∫∫
∑∑ 2
2
2
fal (pa )fbk (pb )hσT v 0 ilk (pa , pb )p2a p2b dpa dpb .
R = 4π(6)2l + 1 2k + 1lk

By reducing the number of Legendre polynomial, compared with number of pitch angle mesh,
we can reduce the computational resource required.
3 Source and sink in the fusion reaction
The fusion products are assumed to be generated with the fusion reaction rate as isotropic
in velocity space and mono-energy neglecting the kinetic energy of fuel ion. On the other
hand, the sink term Sa (pa , θa ) of fuel ions a is given by∫∫
Sa (pa , θa ) = −2πhσT (E)v 0 i(pa , θa , pb , θb )fa (pa , θa )fb (pb , θb )p2b sin θb dpb dθb .
(7)
Using the Legendre polynomials, it can be expressed as∑( 2 )∫∑Pl (cos θa )fa (pa , cos θa )Sa (pa , θa ) = −2πp2b hσT (E)v 0 il,k (pa , pb )fk (pb )dpb .
2k + 1lk(8)Since the double integral over pb and θb is reduced to one-dimensional pb integral, this calculation also requires less computational resource.
4 Benchmark test
We have compared the calculation of the fusion reaction rate with the original quadruple
integral model and the present Legendre expansion model. Two patterns of benchmark test
were carried out. One test was with the Maxwellian distribution function and the otherincludes neutral beam injection.
Table 1: The comparison of computation time

NTHMAX
initial setup time[s]time evolution : time[s]DD fusion reaction rate[×10−5 ]Without Legendre expansion

With Legendre expansion

32
0.790
12.2006.089232
0.750
10.5906.094164
2.870
25.5106.0855128
11.660120.550
6.084664
2.800
19.6806.0867128
10.72062.1706.0849In the test with Maxwellian distribution function, the initial parameters at the center are
electron density 1.0 × 1020 m−3 , deuterium ion density 0.5 × 1020 m−3 , tritium ion density
0.5 × 1020 m−3 and temperature 20keV. Table 1 shows computation time without and with2

Legendre expansion. NTHMAX is number of the pitch angle division. According to Table 1,
when NTHMAX is large, the eﬀect of speedup in calculation is observed.
Next, 1MeV deuterium NBI with beam power of 20MW is assumed. The initial parametersat the center are electron density 1.0×1020 m−3 , deuterium ion density 0.45×1020 m−3 , tritium
ion density 0.45 × 1020 m−3 , α particle density 0.05 × 1020 m−3 and temperature 20keV. Fig.1
shows the example of momentum distribution function with NBI. Table 2 shows computational
PPERP

FMIN =-1.2000E+01 FMAX = -1.2009E+00

12

8

4

0

-12

-8

-4

0

4

8

12

PPARA

Fig. 1: Momentum distribution function

time without Legendre expansion, with Legendre expansion (maximum order of Legendre
expansion L is 0) and with Legendre expansion (maximum order of Legendre expansion Lis 4). According to Table 2, the computational time using Legendre expansion becomes lessthan half of the time without Legendre expansion for NTHMAX= 64. DT fusion reaction
rate is calculated with suﬃcient accuracy.
Table 2: The comparison of calculation speed with the NBI
Without Legendre
expansion

With Legendre
expansion (L=0)

With Legendre
expansion (L=4)

8.660
138.9108.58588.290
52.5308.587512.44053.0908.5875initial setup time [s]time evolution: time [s]DT fusion reaction rate [×10−3 ]Therefore, by using the Legendre expansion, the fusion reaction rate is calculated much
faster.
5 Summary and remaining issue
In order to analyze the time evolution in burning plasma start-up, it is required to describe
time-evolution starting from around 2keV to around 20keV in temperature in the presence ofadditional heating (NBI, ICRF) and fusion reaction. In this situation, the kinetic transport
analysis with the momentum distribution functions requires high resolution in momentum
space in order to describe collisional relaxation between electrons and ions and between bulk
and tail accurately. For such a large scale analysis, it is important to improve the computational performance of calculating the fusion reaction rate in addition to the non-linear
collisional terms. In this study, we have shown that the calculation of source and sink terms
in fusion reaction rate can be accelerated by the Legendre expansion.
The remaining issue is to use the TASK/FP analysis for the startup of ITER burning plasma
and to compare with those of the analysis using 1D diﬀusive transport code TASK/TR.
3

Hydrogen Isotope trapped in tungsten vacancy
Kazuhito Ohsawa1, Yuji Hatano2, Masatake Yamaguchi3

1Research Institute for Applied Mechanics, Kyushu University
2Hydrogen Isotope Research Center, University of Toyama

3Center for Computational Science and e-system, Japan Atomic Energy Agency

1. Introduction
Divertor armor tiles installed in fusion reactors are exposed to the most intense
plasma particle irradiation. So, it is planning that the divertor armor tiles are coveredby tungsten (W) or its alloys because they are excellent properties as plasma facing
materials (PFMs). For example, they have low hydrogen (H) solubility, high thermal
conductivity, and high melting point. However, it is concerned that a large amount of H
isotopes, deuterium (D) and tritium (T), are retained in the W materials in theirradiation circumstance. In particular, tritium retention of the PFMs is very serious
safety concerns because T is a radioisotope whose physical half-life is 12 years.
Usually, H isotopes are not retained in bulk W because H dissolution in W is
considerable endothermic reaction. However, it is expected that vacancy type latticedefects are nucleated by the plasma particle irradiation and H isotopes can beaccommodated in the defects. In the present work, we focus on the H isotopes trapped in
W mono-vacancy. The binding energy between H isotopes and W mono-vacancy isestimated by first-principle calculations on the basis of density functional theory. The
zero point energy of H isotope depends on the mass number. The binding energy oflighter H isotope to a mono-vacancy is larger than that of heavier one by taking into
account the zero point energy. We calculate the amount of vacancy-hydrogen (V-H)
cluster to perform a quantitative analysis to estimate H isotopes trapped in the V-Hclusters in terms of a thermodynamic model assuming thermal equilibrium.
2. Methodology
First-principle calculations are performed using Vienna ab initio simulation package(VASP). Binding energy eka of hydrogen isotopea

H to a W mono-vacancy where aindicates mass number of hydrogen isotope (a=1,2,3) is defined , as follows.

[

]

(

[

]

)

eka = E[Wn-1 V ] − E Wn-1 V a H k − Z[ a H V k ] − k E[Wn ] − E Wn a H T − Z[ a H T ](1)

The function E is the cohesive energy of the supercell and Z is zero point energy ofhydrogen isotopes H . V represents a W mono-vacancy. For example, Wn-1V Hkaa

indicates a supercell composed of n-1 W, a vacancy V, and k aH. Z[aHVk] is the total zero
point energy of k aH trapped in the vacancy. WnaHT means a solute aH in bulk W at T site
and Z[aHT] is the zero point energy of the aHT. We use a large simulation cell, 6 × 6 × 6
bcc lattice i.e. n=432 to calculate vibrational entropy, as mentioned later.

V-H cluster concentrations are calculated in an equilibrium thermodynamic model.
Helmholtz free energy is(2)

F = U − TS
The internal energy is expressed as

U = Nev − ∑ nij eij ,

(3)

ij

where ev is empty vacancy formation energy. nij is the number of V-H clusters
composed of two types of hydrogen isotopes, VaHibHj. The binding energy eij is assumed to
be a weighted average, where eka and ekb are defined in eq. (1).

()

eij = ieia+ j + jeib+ j / (i + j )(4)

N is the total number of the V-H clusters. Entropy is composed of two terms,
configurational and vibrational entropy.
(5)

S = S con + S vib

The configurational entropy is derived from Boltzmann equation.

S con = k B ln Ω ,(6)where Ω is the number of possible arrangement of V-H clusters and solute H isotopes
distributed over the T site in bcc W lattice. The vibrational entropy is described in thelattice vibration ratio in high temperature approximation

S vib = k B ∑ nk ln PkkPk = ∏ν i /ν ki′ ,(7)iwhere ν i and ν ki′ are vibration mode of perfect lattice and i-th nearest neighbor atomin the vicinity of V-H cluster VHk..
We investigate H isotope trapping in two kinds of model. (Model I) The total numberof vacancies, N, in bulk W is assumed to be constant. Therefore, the vacancies are
introduced by plasma particle collision. (Model II) The vacancy formation energy is
considerably reduced by H presence, i.e. superabundant vacancy formation induced by
H isotopes is taken into account. The number of each V-H cluster nij is determined byminimization of the free energy in eq. (2) in the thermodynamic model.
3. Results
Figure 1 indicates binding energies of H isotopes eka to a W vacancy in eq. (1). The
formation energy for an empty vacancy ev is estimated to be 3.028 eV. As mentioned
before, binding energy of lighter H isotopes is larger than that of heavier one. According
to Fig.1, vacancy formation energy is quite reduced by hydrogen trapping. If the numberof H atoms is more than three, the V-H cluster formation energy become negative, thatis, the V-H clusters are spontaneously nucleated in the present thermodynamic model.

Figure 2 shows the total amount of H
isotopes trapped in W vacancies in thecase of model I, where the total number ofvacancies is assumed to be constant. As an
example, the vacancy concentration, i.e.,
number ratio of vacancy and W atom, isassumed to be 1 × 10 −5 . Each H isotope
initially induced in the bulk W is 5 × 10 −4 .
Figure 3 indicates the concentrations of
V-H cluster as a function of bulk H
concentration and temperature in the caseFig. 1: Binding energy of H isotopes(H, D, T) to a W mono-vacancy.

of model II. In the limit of low H concentration, the V-H cluster concentration at a given
temperature is usual thermal vacancy one. However, the V-H clusters rapidly increase
when the bulk H concentration exceeds a threshold value. It is possible to beaccommodated more than six H atoms in a W vacancy. However, the majority of V-Hcluster is VH6 in a wide range of temperature and bulk H concentration.
4. Discussion
Figure 2 suggests that the amount of lighter H isotope trapped in W vacancies ismore than that of heavier one, which reflects the binding energy difference, as shown in
Fig. 1. This effect will be useful to remove or reduce tritium retained in the FPMs
installed in fusion reactors.

Fig. 2: Total amount of H isotopes trapped Fig. 3: Dependence of V-H cluster
in W vacancies in finite temperature. Red concentration on bulk H concentration
and blue lines correspond to H-D and D-T in finite temperatures.
combination, respectively.
References
K. Ohsawa, F. Nakamori, Y. Hatano, M. Yamaguchi, J. Nucl. Mat. 458 (2015) 187.

Eigenmode analysis of geodesic acoustic modes driven by fast ions in the limitof high safety factor
M. Sasaki1,2, N. Kasuya1,2, K. Hallatschek2,3, K. Itoh2,4, Y. Kosuga2,5, M. Lesur6 and S.-I. Itoh1,21Research Institute for Applied Mechanics, Kyushu University, Kasuga 816-8580, Japan
2
Research Center for Plasma Turbulence, Kyushu University, Kasuga 816-8580, Japan
3
Centre for Interdisciplinary Plasma Science, Max-Planck Institut fur Plasmaphysik, EURATOM-IPP
Association, D-85748 Garching, Germany
4
National Institute for Fusion Science, Toki 509-5292, Japan
5
Institute for Advanced Study, Kyushu University, Fukuoka 812-8581, Japan
6
Lorraine University, Institut Jean Lamour 54-506, France
Here we report a new mechanism of excitation of geodesic acoustic modes (GAMs) in the presence of fastions in the limit of high safety factor. Resonance to the magnetic drift of fast ions gives an unstable branchof GAMs driven by fast ions. One of the prominent characteristics of the mode is its spatial structure,
which has a steep poloidal gradient. The associated large wavenumber can cause effective ion heating
mediated by the energy transfer from the GAMs.

Importance of geodesic acoustic modes driven
by energetic particles (EPs), E-GAMs, [1] hasbeen recognized in the study of magnetic
confined plasmas. A gyrokinetic simulation has
shown that the E-GAMs couple with turbulence
so that it may be possible to control turbulent
transport by EPs [2]. GAMs transfer their energy
to bulk ions through dissipation such as Landau
damping, and the energy channeling from EPs tobulk ions via GAMs has been proposed as GAM
channeling [3]. Actually, E-GAM have been
observed in experiments [4], and several branchesof GAMs have been found [5]. Identification ofeach branch of E-GAMs is required.
In this study, the eigenmode analysis of
E-GAMs is performed based on the gyrokinetic
equation in the limit of high safety factor. We
show there are several branches of GAMs in thepresence of EPs. In accordance with the energy ofEPs, the resonance between particles motion and
the GAMs changes so that a new branch appears.
Here, we report a new branch of E-GAM.
The plasma configuration under consideration isa high aspect ratio tokamak with a circular crosssection. The linearized gyrokinetic equation for the
non-adiabatic response of the ion velocity
distribution to a GAM oscillation is solved in thelimit of high safety factor. The response function of
the EPs has two kinds of resonances; (i): theresonance due to the transit frequency, and that due
to the drift frequency. In the limit of high safety
factor, the resonance due to the drift frequency
becomes dominant, and a new branch with afrequency close to the drift frequency of the fastions appears. In this case, the poloidal side band
coupling is strong. The poloidal mode expansion

cannot be used, so that the eigenfunction is treatedas a function of the poloidal angle. The poloidal
eigenfunction of the unstable branch, which isrelated to the electrostatic potential fluctuation, isshown in Fig. 1. The real part of the eigenfunction
has bumps where the resonance of the magnetic
drift with the wave is strong. The imaginary part ofthe eigenfunction becomes important near the
resonance positions, indicating the particles andassociated energy transfer. The effective
wavenumber parallel to the magnetic field, which isevaluated from the poloidal angle derivative,
increases due to the bumps, and the effect of ion
heating through Landau damping becomes large.

RealImaginaryFig.1. Poloidal eigenfunction of the new branch[1] G. Y. Fu, Phys. Rev. Lett., 101, 185002 (2008).
[2] D. Zarzoso, et. al., Phys. Rev. Lett., 110, 125002(2013).
[3] M. Sasaki, et. al., Plasma Phys. Control. Fusion, 53,
085017 (2011).
[4] R. Nazikian, et. al., Phys. Rev. Lett., 101, 185001(2008).
[5] T. Ido, et. al., IAEA Fusion Energy Conference,
PD/P8-16 (2012).

Analysis of triplet term in parallel momentum flux
Y. Kosuga[1,2] , S.-I. Itoh[2,3] , P.H. Diamond[4] , K. Itoh[5,3][1] Institute for Advanced Study, Kyushu University, Fukuoka, Japan[2] Research Institute for Applied Mechanics, Kyushu University, Fukuoka, Japan[3] Research Center for Plasma Turbulence, Kyushu University, Fukuoka, Japan[4] University of California, San Diego, USA[5] National Institute for Fusion Science, Gifu, Japan
Transport of the momentum along the magnetic field (or shortly, parallel momentum
transport) is an important issue for magnetic fusion. In fusion plasmas, these parallel
flows (or approximately toroidal flows) can stabilize harmful MHD instability and control
transport by contributing to the E × B shear. While parallel flows are externally driven
by NBI, recent research has revealed the intrinsic property of magnetized plasmas todrive parallel flows intrinsically[1,2]. Namely, parallel flows are generated in the absenceof appearant momentum source, such as in RF heated plasmas. In these cases, it isoften thought that the flows are driven by turbulence. Direct measurement of parallel
momentum flux has been performed, and the role of turbulence in driving flows have been
confirmed[3].
A recent development of the problem of parallel momentum transport is the observation of the triplet correlation in the total momentum flux. The total momentum flux isgiven byΠxz = hnihṽx ṽz i + hvz ihṽx ñi + hñṽx ṽz i(1)Here, x is the radial direction and z denotes the direction of the magnetic field. The
first term is Reynolds stress, the second term is convective term, and the last term is thenonlinear flux. Recent measurement on the linear device PANTA confirm the contribution
from the triplet term[4]. The triplet term can be a major part of the total momentum
flux, as reported from the TORPEX experiment[5]. Then it is an important issue toanalyze this term and clarify the parameter dependence, in order to asses it impact on
different plasmas.
From modeling perspective, the triplet flux requires the closure calculation. In a recent
work[6], the triplet was calculated for the case of strong turbulence, where blobs are amajor player in the turbulent fluctuation. On the other hand, turbulence often contains
waves, as indicated by the sharp peaks in the frequency spectrum. In this case, we needto address the role of waves in driving the triplet flux[7]. In this case, the necessary phase
is for producing the triplet flux is set by the wave resonance, rather than by the nonlinear
decorrelation. As a result, since drift waves are dispersive, we would expect that the classof triads that satisfy the frequency resonance condition is very limited. We can identify
triads that contribute to the nonlinear flux dominantly. We later discuss this in moredetail.
In order to calculate the triplet flux, we use Hasegawa-Mima equation with parallel
flow coupling. Within the model, we focus on the wave momentum density,2hq̃ṽk i ∼ kk (1 + k⊥)|φ̂k |2 /ωk ≡ Pkk(2)

Here we used normalized quantities. The space is normalized by ρs , the time is normalizedby ωci , and the potential is normalized by e/Te . q = (1 − ∇2⊥ φ). Note that Ek =

2)|φ̂k |2 is the energy density of drift waves. So Pkk = Ek /(ωk /kk ) is a momentum(1 + k⊥density in the parallel direction. This may be understood from the particle mechanics,
where momentum p is obtained by dividing energy by velocity, p = E/v. Note that thesign of the Pkk is set by the symmetry breaking in the parallel direction of drift waves. Thiscan be provided from several processes, including the intensity gradient, radial electricfield shear, etc.
The wave momentum density can be spatially transported. One way is due to thewave propagation. In this case, the flux is given by vgr Pkk , where vgr is the group velocity
of drift waves. Typically vgr > 0 for drift waves, due to the out-going wave boundary
condition. The flux of wave momentum by the drift wave propagation exert residual stressand drive parallel flows.
In addition to the linear propagation, the momentum density can be spatially scattered by turbulence and can spread. This effect may be captured by the nonlinear, triplet
contribution hṽx q̃ v˜k i. A closure calculation is necessary to calculate this term. By employing the two scale direct interaction approximation, we can calculate the nonlinear flux
as:X
hṽx q̃ṽk i =(−Dk1 ∂x Pk1 k + Vk1 Pk1 k )(3)
k1

Here Dk1 is the scale dependent diffusivity and Vk1 is the scale dependent convective
velocity. They are defined asX
Dk1 =θk1 ,k2 |φ̂k2 |2 F1 (k1 , k2 )
(4)
k2

Vk1 =

X

θk1 ,k2 ∂x |φ̂k2 |2 F2 (k1 , k2 )

(5)

k2

Here θk1 ,k2 ∼= πδ(ωk1 +k2 − ωk1 − ωk2 ) is the triad interaction time. F1 and F2 are coupling
coefficients between k1 and k2 . As an estimate, if we use gyro-Bohm like diffusivity andthe convective velocity, we find thatρs(6)
hṽx q̃ṽk i ∼ −v∗ kk Nk
LI
Here LI is the intensity gradient and the sign convention is such that when the fluctuation
intensity increases toward the edge (LI > 0) the flux is toward the core.
From the approximated form, we find that the flux becomes important around region
with steep intensity gradient. For example, the nonlinear flux is important in the region
with transport barriers, since across the barrier the fluctuation intensity varies dramatically. So the triplet flux is important to understand the generation of flows in suchregions.
We acknowledge useful discussion with Drs. S. Inagaki, T. Kobayashi, L. Wang, andZ. Guo. This work was supported by Grants-in-Aid for Scientific Research of JSPF ofJapan (23244113, 15H02155, 15K17799), Kyushu University Interdisciplinary Programsin Education and Projects in Research Development (TT26705).

References[1] P.H. Diamond, et al., Nucl. Fusion 53 104019 (2013)

[2] K. Ida and J.E. Rice, et al., Nucl. Fusion 54 045001 (2014)[3] Y. Xu, et al., Nucl. Fusion 53 072001 (2013)[4] S. Inagaki, et al., submitted to Sci. Rep. (2016)[5] B. Labit, et al., Phys. Plasmas 18 032308 (2011)[6] L. Wang, et al., Phys. Plasmas 22 052302 (2015)[7] Y. Kosuga, et al., in preparation (2016)

Transport analysis of impurity LHD peripheral regions during neon gas-puff
Gakushi Kawamura, Masahiro Kobayashi, Shuyu Dai, Kiyofumi Mukai,
Ryuichi Sano, Shwetang N. Pandya1, Byron Peterson, Yuhe Feng2,
and the LHD Experiment Group
National Institute for Fusion Science, 322-6 Oroshi-cho Toki, Gifu 509-5292, Japan
1
Institute for Plasma Research, India
2

Max-Planck-Institut für Plasmaphysik, Greifswald Germany

Transport analysis of neon seeded by gas puff was carried out. LHD
discharge #117478 was simulated by EMC3-EIRENE code and radial
distribution of plasma measured by Thomson scattering was reproduced.
The results suggest that neon accumulates more than carbon because it has
deeper penetration of neutral particles and accesses to the region where
thermal force is dominant. The transport of neon ions along a magnetic field
line was strongly affected by the force balance between friction and thermal
force. The friction force determined by the plasma flow speed can beenhanced by locally strong of magnetic flux density similar to a magnetic
nozzle. The simulation suggests that helical ripple can enhance impurity
screening.
Heat load on divertor plates is one of critical issues in a fusion device and a future reactor.
Methods of power reduction such and advanced divertor concepts, impurity seeding,detachment are proposed and actively studied. Impurity gas puff is one of possible methods andhas been investigated experimentally in LHD [1]. Direct measurement of impurity transport in
the plasma is difficult and thus transport modeling with a simulation code is necessary tounderstand it. We have developed transport model of plasma, impurity and neutral particles in
the peripheral plasma by using EMC3-EIRENE [2,3] and applied to LHD [4-6]. The code solves
fluid equation with Monte-Carlo technique and has the capability of resolving perpendiculartransport across a magnetic field line in 3D space. We modeled a 1/20 toroidal block of the LHD
plasma, which consists of plasma and vacuum regions and they cover every cross section of thevacuum vessel. The plasma mesh system of EMC3 covers the peripheral plasma which has longconnection length, roughly >10m. The core region is excluded from the simulation and modeled
as boundary conditions at the last closed flux surface (LCFS). The magnetic axis is located at
Rax=3.6m and the vacuum field, i.e., zero beta value, is employed to make the calculation mesh.
We analyzed a LHD discharge with neon gas puff by EMC3-EIRENE. Input parameters werebased on the discharge #117478 (17th cycle, 2013) and determined to reproduce radial
distribution of plasma (see Fig. 1): input power Pin=11.6MW, electron density at core-SOL
boundary ne=3.3×1019/m3, perpendicular transport coefficients D=Dimp=1.0 m2/s and χe=χi=0.4
m2/s. The duration time of the gas puff is approximately 100ms and small increase of electron
density, 10~20%, and large increase of radiation, approximately twice, were observed. The
plasma become steady state after short transition time and the enhanced radiation continued
until the end of the discharge. That fact indicates that some amount of neon gas was confined inthe device and recycled through ionization and surface recombination in steady state. Therefore
we modeled the neon as no gas puff but surface source on the divertor plates and no wall loss in
EMC3-EIRENE code. The amount of neon is determined to have the same amount of total
radiation as that measured by bolometer, 2.5MW. We simulate carbon and neon, which have
0.24MW and 1.4MW radiation, respectively, and the remaining 0.86 MW radiation is fromhydrogen.

Fig. 1: Comparison of electron density and temperature between Thomson scattering (red dots)
and EMC3-EIRENE (blue lines). The green lines represent the distribution of connection length.
We compared distribution of density and radiation of carbon and neon and found that muchneon accumulate in the ergodic region than carbon, see Fig. 2(a). The physical reason of it canbe understood from neutral transport and force acting on ions. Neon has higher ionizationenergy and deeper penetration into the plasma. The penetration limit is drawn in Fig 2 (a) as
white lines. Carbon has high density outside the limit, e.g. left side of the white line, and neonhas mainly inside the limit. This limit is determined by the distribution of electron temperature.
Figure 2(b) shows distribution of force balance on impurity ions. Positive value (red) means
friction force dominant and negative value (blue) means thermal force dominant. Friction force
pushes ions toward downstream and thermal force is opposite. Two black lines are thepenetration limit of neutral particles. Neon has deeper penetration and can access to the thermal
force dominant region and therefore neon is transported toward upstream more easily thancarbon. That is the reason of larger accumulation in the ergodic region.
(a)

(b)

Fig. 2: distributions of (a) C
and Ne density and (b) force
balance on impurity ions. The
white and black lines in thebothfigures
represent
penetration limit of neutral
impurity particles.

Helical plasma has complicated structure of magnetic field lines and impurity transport is also
complicated in space. Figure 3 shows density distribution of neon ions along a magnetic fieldline. The field tracing was carried out from the density peak in Fig. 2(a) toward the both
directions, e.g. positive and negative. Thermal force is dominant in field line except near thedivertor plates, where fast plasma flow exists because of collisional presheath. The neon ionizedin the divertor plasma experiences successive ionization to higher states and transported toward
the upstream. In the middle of the field line, flow speed of impurity becomes zero and causes alarge peak. Ions flow from the both side and escape from the field line due to perpendiculartransport. Another field line has large area where friction force is dominant. Such area plays arole of barrier of impurity transport, see Fig. 4. The plasma flow speed is essential factor for
friction force and affected by distribution of magnetic flux density. Helical plasma hasperiodical change of the magnetic flux density, e.g. helical ripple. We confirm that the ripple
plays a role of magnetic nozzle and causes locally enhanced flow speed. That causes locallylarge friction force and pushes impurity ions toward the downstream. It can contribute toimpurity screening.

Fig. 3: distributions of neon
density with three chargestate groups (green, blue andmagenta lines) and force
balance (red line) along amagnetic field line.

Fig. 4: distributions of neon density of all
charge states (blue line), force balance (redline) and magnetic flux density (greenline). The distribution of magnetic flux
density and force balance has a strong
correlation.

References[1] K. Mukai et al., Nucl. Fusion 55 (2015) 083016[2] Y. Feng et al., Contrib. Plasma Phys. 44 (2004) 57.
[3] D. Reiter et al., Nucl. Fusion 47 (2005) 172.
[4] G. Kawamura, et al., Contrib. Plasma Phys. 54 (2014) 437.
[5] M. Kobayashi, et al., Fusion Sci. Technol. 58 (2010) 220.
[6] M. Kobayashi et al., Nucl. Fusion 53 (2013) 033011.

Bohm Criterion and Virtual Divertor Model for SOL-Divertor Simulation
T. Takizuka1, S. Togo2, S. Azuma3, A. Fukuyama3
1

2

Graduate School of Engineering, Osaka University
Graduate School of Frontier Sciences, University of Tokyo
3
Graduate School of Engineering, Kyoto University

Bohm criterion, Vf ≥ Cs (Vf : plasma flow speed at the sheath entrance, Cs : sound speed), hasbeen derived from the sheath formation condition. When the Bohm criterion is analyzed for
two-ion-species plasma, two branches of the criterion are derived. Numerical experiments with
PARASOL, however, give always the upper branch but never the lower branch. Introducing alinear stability analysis, it is found that the upper branch is stable but the lower branch unstable.
Also found that the Bohm criterion is not directly related with the sheath formation condition.
Generalized fluid model with anisotropic ion temperature for SOL-divertor simulation does not
require a down-stream boundary condition on Vf at the divertor plate. For this sake, a virtual
divertor model has been developed. Simulations with a new code demonstrate that Vf at thedivertor plate exceeds Cs automatically without using the boundary condition of Bohm criterion.

1. Introduction
Bohm criterion is well known as a boundary condition for the plasma flow speed Vf at the sheathentrance in front of a plasma facing material. A simple criterion, Vf ≥ Cs (Cs: sound speed), is given
from a condition that the solution to Poisson’s equation for the sheath potential is not an oscillating
profile but a monotonic profile [1]. From the plasma fluid equations, a linear response of the electron
density ne to the sheath potential φ is derived as δne/n0 = e δφ /Te ≡ Φ and that of ion density ni as δni/n0= Φ {Te /(miVf2 – γaTi)}, where γa is the adiabatic index. The Poisson’s equation is then described by
d2Φ /dx2 = κ2 Φ = (1/λD2){1 – Te /(miVf2 – γaTi)} Φ, where λD is the Debye length. For the monotonic
profile of Φ ~ exp (κ x) with κ2 ≥ 0, the condition miVf2 ≥ Te + γaTi is required. This elementary theory
on the sheath formation has been believed as a base of the Bohm criterion.
Applying the above analytical method, the Bohm criterion for two-ion-species plasma was investigated[2]. The analysis gives two curves like a shifted hyperbola in a (V1, V2) diagram. On an upper branch,
V1 > Cs1 and V2 > Cs2, while on a lower branch V1 < Cs1 and V2 < Cs2, where Vj and Csj are the flow
speed at the sheath entrance and a specific sound speed for species-j ion (j = 1, 2), respectively.
Analytical results were compared with PARASOL simulation results [2]. The PARASOL simulation
showed that V1 and V2 are always located near the upper branch. The study in Ref. [2] is briefly
reviewed in the next section.
In order to solve the reason why the lower-branch solution is absent, we have introduced a linear
stability analysis to the sheath theory. We found that the upper branch is stable but the lower branch
unstable [3,4]. This is because the upper branch is realized but the lower branch does never appear in
the PARASOL simulation. This stability analysis in Ref. [3,4] is shown in Section 3.
Recently we have developed a new code for SOL-divertor simulation with a generalized plasma
fluid model taking account of anisotropic ion temperatures [5,6]. The momentum transport equation is
described by a spatially 1st-order differential equation without the 2nd-order viscosity term. Thisequation does not especially require a boundary condition at the divertor plate. In order to simulate theexistence of the divertor plate, we have introduced a virtual divertor (VD) model [5,6]. Simulation
results with our new code demonstrated that the flow speed, Vf, at the divertor plate exceeds the soundspeed, Cs, automatically. The criterion Vf ≥ Cs is required for the temporal-stability condition at theplasma boundary. When Vf < Cs , Vf is temporally increased up to Cs until the stability condition issatisfied. The topics of the VD model and Bohm criterion are presented in Section 4.

2. Bohm criterion for two-ion-species plasma
The Bohm criterion for two-ion-species plasma was investigated [2]. Basic equations in the analysis
are; momentum balance equation, ∂neTe/∂x = e ∂φ /∂x, for electron, and equation of continuity, ∂njVj /∂x = 0,
and momentum balance equation, ∂/∂x (njmjVj2 + njTj) = – e ∂φ /∂x, for species-j ion. Linear responses ofne, n1 and n2 to φ are derived as δne/n0 = Φ, δn1/n0 = F1Φ
× {Te /(m1V12 – γaT1)} and δn2/n0 = F2Φ {Te /(m2V22 – γaT2)},
where Fj = n0j /n0 is the density fraction of species-j ion
and F1 + F2 = 1. Poisson’s equation is given by d2Φ /dx2 =κ2 Φ = (1/λD2){1 – Te /(m1V12 – γaT1) – Te /(m2V22 – γaT2)} Φ.
The monotonic sheath of Φ ~ exp (κ x) can be formed when
κ2 ≥ 0. This condition is rewritten as(V12 – Cs12) (V22 – Cs22) ≥ F1 F2 Te2 /m1 m2 ,(1)

where Cs12 = (Fj Te + γaTj)/mj is the specific sound speed forspecies-j ion. Condition (1) is satisfied inside blue regions
Fig. 1 Regions satisfying Condition (1)in the (V1, V2) diagram shown in Fig. 1. On an upper branch,in (V1, V2) diagram for flow speeds at theV1 > Cs1 and V2 > Cs2, while on a lower branch V1 < Cs1
sheath entrance for two-ion-species plasma.
and V2 < Cs2. This analytical result was compared with the
1D PARASOL simulation results for various mass ratio m2/m1, collisionality L// /lmfp etc. Figure 2
shows results of a numerical experiment using PARASOL code [7]; m2/m1, = 4 and L// /lmfp = 0.01. The
flow speeds, V1 (red line) and V2 (blue line) at the sheath entrance are measured. It was found from the
PARASOL simulation that V1 and V2 are located near the upper branch as plotted in Fig. 3 by a cross (×),

Fig. 2 Profiles of two-ion-species SOL plasma in a 1D PARASOL simulation, (a) potential φ and (b) flow
velocity V// parallel to magnetic field. The hot-particle source is given in the central SOL region, and thesheath is formed in front of the divertor plate.

Fig. 3 Flow speeds at the sheath entrance, V1 and V2, measured in PARASOL simulations with m2/m1 =4 for (a) rare-collision case and (b) collisional case. Simulation results (×) are located near the upper
branch. Blue dashed line is V2 = (m1/m2)1/2 V1, and orange dashed line is V2 = V1.

and the lower-branch result was never realized. For the rare-collision case (a) L// /lmfp =0.01, the (V1, V2)
have a relation m1V12 = m2V22 as plotted by a blue dashed line. This is due to the same amount of energy
gain from the electric field in the SOL plasma. On the other hand for the collisional case (b) L// /lmfp =100, (V1, V2) tend to have the same speed as plotted by an orange dashed line. This is simply due to thecollisional relaxation between V1 and V2. In either case, (V1, V2) are located near the upper branch, V1 >Cs1 and V2 > Cs2. In the next section, we study the reason why the lower-branch solution is absent.
3. Stability analyses for the Bohm criterion [3,4]We consider the time-dependent perturbation of the form Φ ~ exp (γ t + κ x). Including the growth rateγ, linear responses of ne, n1 and n2 to Φ are written similar to those for Section 2; δne/n0 = Φ, δn1/n0 =F1Φ {Te /(m1V1*2 – γaT1)} and δn2/n0 = F2Φ {Te /(m2V2*2 – γaT2)}, where V1* = V1 + γ /κ and V2* = V2 + γ /κ.
The condition for the monotonic sheath formation is obtained similar to Condition (1):(V1*2 – Cs12) (V2*2 – Cs22) ≥ F1 F2 Te2 /m1 m2 .

(2)When (V1, V2) are located inside the upper blue region in Fig. 4 (a) (assuming V > 0 and κ > 0), Condition(2) can be satisfied not only for the null growth rate γ = 0 ( (V1*, V2*) = (V1, V2) stay inside the upper
blue region), but also the negative γ (shifting downward by (γ /κ, γ /κ), (V1*, V2*) stay inside the lower
blue region). This means the upper solution of (V1, V2) is stable. On the contrary, for the case of (V1,
V2) inside the lower blue region in Fig. 4 (b), Condition (2) is satisfied for the null growth rate γ = 0( (V1*, V2*) = (V1, V2) stay inside the lower blue region), and also for the positive γ (shifting upward by(γ /κ, γ /κ), (V1*, V2*) stay inside the upper blue region). This means the lower solution of (V1, V2) is
unstable. Accordingly the stable upper branch is realized but the unstable lower branch does neverappear in the PARASOL simulation.

Fig. 4 Stable upper-branchsolution (a) and unstable
lower-branch solution (b)of (V1, V2).

We can understand that the Bohm criterion is a result of the stability condition of plasma fluids.
Now we examine whether the Poisson’s equation is essential or not for the Bohm criterion. For
simplicity we treat a single-species-ion plasma. In the quasi-neutral plasma (λD2 → 0), a relation δne =δni holds, and the following equation is obtained:( κ V f + γ ) 2 – κ 2 C s2 = 0 ,(3)where Cs2 = (Te + γaTi)/mi. At the outlet of SOL plasma flow, the perturbation is increased monotonicallytowards the outlet, i.e., κ Vf > 0. For the stable perturbation (γ ≤ 0), the Bohm criterion, Vf2 ≥ Cs2, isrequired. Note that Eq. (3) is a dispersion relation for a sound wave with Doppler shift, (ω – k Vf) )2 = k2 Cs2,
when (γ, κ) are replaced with (– i ω, i k).
In this section, it is clarified that the Bohm criterion is fundamentally realized from the stability
condition of the plasma fluid flowing out to the outlet. Also proved that the Bohm criterion is notdirectly connected with the sheath formation. This finding is a backbone of the SOL-divertor
simulation with a virtual divertor model (see next section), which does not require the Bohm criterion
as a boundary condition at the divertor plate.

4. Virtual divertor model and Bohm criterion
We have developed a new code for 1D SOL-divertor plasmas with a generalized fluid model taking
account of anisotropic ion temperatures, Ti// and Ti⊥ (parallel // and perpendicular ⊥ to the magnetic field)[5,6]. The parallel momentum transport equation is described by a spatially 1st-order differential equation
without the 2nd-order viscosity term. This equation does not especially require a boundary condition at
the divertor plate. In order to simulate the existence of the divertor plate, we have introduced a virtual
divertor (VD) model [5,6]. The VD region is set beyond the divertor plate position (blue vertical line)
as shown in Fig. 5 (a). Plasma fluid equations are solved
continuously from the SOL-divertor region to this VD
region, where artificial sinks for particle, momentum,
and energy are added. By coupling the anisotropic-iontemperature formulation and the VD model, simulation
results demonstrated that the flow speed, Vf, at the divertor
plate exceeds automatically the sound speed, Cs (or Mach
number M ≡ Vf/Cs ≥ 1). The criterion M ≥ 1 is required
for the temporal-stability condition of the monotonicprofile plasma formation at the plasma boundary. Strong
sinks in the VD region generate the monotonic change inthe plasma profile at the VD entrance. When Vf < Cs, Vf
temporally grows up to Cs until the stability condition issatisfied as shown in Fig. 5 (b).
Using this new 1D code, we are studying various
physics of SOL-divertor plasmas in tokamaks, such as
Fig. 5 (a) Profile of Mach number M in a
the detached divertor physics [8]. Application of the new
SOL-divertor plasma parallel to magnetic
field. (b) Evolution of M at divertor platemodels to comprehensive divertor simulation codes is(s = 2m) and M x at x-point (s = 1m).
the future work.
Acknowledgements The authors thank Prof. Y. Ogawa (U. Tokyo), Dr. K. Ibano (Osaka U.), Dr. K.
Hoshino (JAEA) and Dr. M. Nakamura (JAEA) for the collaboration on virtual divertor modeling. One ofauthors, TT, is grateful to Prof. Y. Ueda of Osaka University for his continuous encouragement.

References[1] K.-U. Riemann, J. Phys. D: Appl. Phys. 24 (1991) 492.
[2] S. Azuma, A. Fukuyama, T. Takizuka, Contrib. Plasma Phys. 52 (2012) 512.
[3] S. Azuma, “Particle simulation of peripheral plasmas in tokamaks”, Master thesis, Graduate School ofEngineering, Kyoto University (2014).
[4] T. Takizuka, S. Togo, S. Azuma, A. Fukuyama, “Bohm criterion and virtual divertor model for SOLdivertor simulation”, presented at 15th Int. Workshop on Plasma Edge Theory in Fusion Devices (PET-15)(9-11 Sep 2015, Nara, Japan).
[5] S. Togo, T. Takizuka, M. Nakamura, K. Hoshino, Y. Ogawa, J. Nucl. Mater. 463 (2015) 502.
[6] S. Togo, T. Takizuka, M. Nakamura, K. Hoshino, K. Ibano, T. L. Lang, Y. Ogawa, “Self-consistenttreatment of the sheath boundary conditions by introducing anisotropic ion temperatures and virtual
divertor model”, to be published in J. Comput. Phys.
[7] T. Takizuka, Plasma Sci. Technol. 13 (2011) 316.
[8] S. Togo, T. Takizuka, M. Nakamura, K. Hoshino, K. Ibano, T.L. Lang, Y. Ogawa, “Simulation studyof detached plasmas by using one-dimensional SOL-divertor fluid code with virtual divertor model”,
presented at PET-15; to be published in Contrib. Plasma Phys.

Numerical analysis of the eﬀects of radial plasma diﬀusion on the
detached divertor plasmas
Satoshi Togo, Tomonori Takizukaa) , Makoto Nakamurab) , Kazuo Hoshinob) , Kenzo Ibanoa) , Yuichi
Ogawa
Graduate School of Frontier Sciences, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa 277-8568, Japana)
Graduate school of Engineering, Osaka University, 2-1 Yamadaoka, Suita 565-0871, Japanb)
Japan Atomic Energy Agency, 2-166 Omotedate, Obuchi-aza, O-aza, Rokkasho 039-3212, Japan
Radial plasma diﬀusion is considered to be important as a candidate to reduce the particle flux to the divertor
plate in addition to the volume recombination. We investigated the eﬀects of radial plasma diﬀusion on the
detached divertor plasmas by using a one-dimensional (1D) scrape-oﬀ-layer (SOL)-divertor (DIV) plasma fluid
model incorporating anisotropic ion temperature. A virtual divertor (VD) model is implemented in our code inorder to reproduce the eﬀects of divertor plate and sheath. The roll-over of the particle flux to the divertor plate
becomes smoother with radial plasma diﬀusion.

1

Introduction

Development of accurate simulation models of scrape-oﬀlayer (SOL)-divertor (DIV) plasma is crucial for the handling of the divertor heat load in the future tokamaks. SOLDIV code packages such as SOLPS [1], SONIC [2, 3] andUEDGE [4] have been developed and widely used in designing the divertor and predicting the heat load. It is reported, however, that the simulation results of these code
packages do not satisfactorily agree with the experimental
results [5, 6].
Radial plasma diﬀusion is considered to be important as a candidate to reduce the particle flux to the divertor plate in addition to the volume recombination [7].
Moreover, it was shown by a Particle-in-Cell (PIC) codewith binary collision model that supersonic flows could
be driven by the eﬀect of radial plasma diﬀusion which
was considered to become strong in the partially detacheddivertor plasmas with short radial decay length and help
the plasma transit from attached to detached regime moresmoothly [8]. A one-dimensional (1D) SOL-DIV plasma
fluid simulation code with inviscid approximation showedthat supersonic flow was observed in the low DIV temperature conditions and that the particle flux to the divertor plate became lower than that predicted with two-point
model with subsonic flow [9]. In a simulation with a code
package, on the other hand, the Mach number at the sheathentrance Mt = 1 which was given as a boundary condition was changed into Mt > 1 in order to drive supersonicflow artificially. The particle flux, however, did not change
largely [10]. We suspected that this diﬀerence was brought
about by the existence of parallel viscosity term in the parallel momentum transport equation.
We investigated the eﬀects of radial plasma diﬀusionon the detached divertor plasmas by using a 1D SOL-DIV
author’s e-mail: togo@ppl.k.u-tokyo.ac.jp

plasma fluid model incorporating anisotropic ion temperature T i,k , T i,⊥ [11–13]. Because the parallel viscosity termin the parallel momentum transport equation is an approximated form of the ion stress tensor πi ≡ 2n(T i,k − T i,⊥ )/3 inthe assumption of πi << nT i ≡ n(T i,k + 2T i,⊥ )/3, the sum ofnT i and πi is replaced by nT i,k and the parallel momentum
transport equation becomes the first-order diﬀerential [15].
Therefore, the flow velocity is determined only by the upstream condition and the boundary condition for Mt is nolonger necessary. Instead, we implemented a virtual divertor (VD) model in our code in order to reproduce theeﬀects of divertor plate and sheath [11–14].
In this paper, the numerical model is explained in Sec.
2 and the results are shown in Sec. 3. In Sec.4, summary
and future works are described.

2

Numerical Model

2.1

Plasma fluid model

The 1D fluid equations with anisotropic ion temperature
T i,k , T i,⊥ in the system of the homogeneous magnetic field
are as follows [11–13];∂∂n(nV) = S ,
+
∂t ∂s
)
∂mi nV∂ (+mi nV 2 + nT i,k + nT e = Mm ,
∂t
∂s
(
)
1
∂ 1
mi nV 2 + nT i,k
∂t 2
2
(
)
∂ 1
3
+
mi nV 3 + nT i,k V + qi,k
∂s 2
2
(
)
(
)
n T i,⊥ − T i,kme n T e − T i,k+= Qi,k +
τrlx
mi
τe
∂nT e
−V
,
∂s

(1)
(2)

(3)

∂nT i,⊥)∂ (+nT i,⊥ V + qi,⊥
∂t
∂s(
)
(
)
n T i,⊥ − T i,k
2me n T e − T i,⊥= Qi,⊥ −
+
,
τrlx
mi
τe
(
)
(
)
∂ 3
∂ 5
nT e +
nT e V + qe
∂t 2
∂s 2
3me n (T i − T e )∂nT e= Qe +
+V
.
mi
τe
∂s

2.3
(4)

The source terms for the plasma and neutrals are as follows;
S = S core + nn νiz − nνrc + S ⊥ + S VD ,inn+ mi ninn
n,recy Vn,recy sin θ (νiz + νcx )1
5
mi nV 3 + nT i V + qi = γi nT i V,
2
2

(6)

5
nT e V + qe = γe nT e V,
2

(7)

2
2
Qi,k = Qcore
i,k + mi nn,recy Vn,recy (νiz + νcx )3
1+ nn,diﬀ T i (νiz + νcx ) + Q⊥i,k
(14)
(2
)
1
1
− mi V 2 + T i,k (nn νcx + nνrc ) + QVD
i,k ,
2
2
Qi,⊥ = Qcore
i,⊥
(
)
4
2
+ mi nn,recy Vn,recy+ nn,diﬀ T i (νiz + νcx ) (15)3
− T i,⊥ (nn νcx + nνrc ) + Q⊥i,⊥ + QVDi,⊥ ,3− iz nn νiz − nT e νrc − Lz rimp n2
Qe = Qcoree2+ Q⊥e + QVD
,
e
(
)
out
out
S n,recy = −nn,recy νcx + νiz + νL,recy

Neutral fluid modelout,VD+ S n,recy
,

(
)
inn
S n,recy= −ninnn,recy νcx + νiz + νL,recy

Neutrals are divided into recycling neutrals generated onouter and inner divertor plate and diﬀusion ones generated
inn
in the plasma whose densities are denoted by noutn,recy , nn,recy
and nn,diﬀ , respectively. This model is similar to the firstflight corrected diﬀusion (FFCD) model in Ref. [18]. The
equations are as follows [12];∂t)∂ ( outout
out+nn,recy Vn,recy= S n,recy,∂x)∂ ( inn
inn
inn
nn,recy Vn,recy= S n,recy
,
∂t
∂x
(
)
∂nn,diﬀ∂nn,diﬀ∂+−Dn= S n,diﬀ .
∂t
∂x
∂x
∂ninn
n,recy

+

(8)

(9)
(10)

Here, the coordinate x is in the poloidal direction andx = s sin θ. We set the pitch angle to be( θ = 6◦ . ) The
out
innvelocities of the recycling neutrals Vn,recy= −Vn,recy
andthe diﬀusivity of the diﬀusion neutrals Dn are the same as
Ref. [12]. Assuming that the diﬀusion neutrals impinging the divertor plates do not reflect and become recycling
neutrals, we impose the following boundary condition at
the divertor plate [12];−Dn∂nn,diﬀ= Cn nn,diﬀ v0 .
∂x

(11)

Here, the coeﬃcient Cn is given from kinetic considerations.

(13)

⊥
VD
− mi V (nn νcx + nνrc ) + Mm
+ Mm,

where qi = qi,k + qi,⊥ , and γi and γe are the sheath heattransmission factors for ion and electron, respectively.

∂noutn,recy(12)out
Mm = mi noutn,recy Vn,recy sin θ (νiz + νcx )(5)Here, we assume that the inner/outer SOL-divertor system is symmetry. The definitions of the notations are
the same as those in Refs. [11–13]. Here, we set theheat flux limiters to be αi,k = αi,⊥ = 0.5 and αe = 0.2.
The relaxation time of ion temperature is assumed to beτrlx =√2.5τi [16, 17]. The sound speed is defined as()
cs =T e + T i,k /mi . The boundary conditions for theheat fluxes of ion and electron are prescribed as follows;

2.2

Source terms

inn,VD+ S n,recy,()
S n,diﬀ = −nn,diﬀ νiz + νL,diﬀ + nn,recy νcxinn,VD+ nνcx + S n,recy
.

(16)

(17)

(18)

(19)

Here, the superscript “⊥” represents the sinks due to radialplasma diﬀusion which are explained in the next section.
We introduce the simple radiative cooling model with radiation rate of carbon Lz [19]. The impurity ratio to theplasma density is denoted by rimp . The definitions of othernotations are the same as Ref. [12]. For the energy loss due
to ionization and excitation, iz = 30 eV was used.

2.4

Radial plasma diﬀusion

Similar to the particle and heat fluxes from the core plasma,
we treat the eﬀect of radial plasma diﬀusion with sinkterms in DIV as follows;nS⊥ = − ⊥,
τ

(20)

mi nV
,
τ⊥
(
)
1 1
⊥
2
⊥ 1
Qi,k = − ⊥ mi nV + gi,k nT i,k ,
τ 2
2
⊥
Mm
=−

Q⊥i,⊥ = −g⊥i,⊥ nT i,⊥
τ⊥

,

(21)
(22)
(23)

2.5
Mt

M

2.0
1.5

MX

1.0
0.5
0.8

0.9
(RΓ/Rp)(Tt/TX)1/2

1

Fig. 2 Mach number at the sheath entrance Mt (circles) and at
the X-point MX (squares) as functions of C. The theoretical values are shown by the dashed lines.
Fig. 1 Schematic of the SOL-DIV 1D system with a VD model.

(()∂nn,diﬀ
nn,diﬀnout,VD= η VD H (−V) + ηn VD H
S n,recy, (31)
∂x
τ
τn

)

1
3
Q⊥e = − ⊥ g⊥e nT e ,
τ
2

(24)

nn,diﬀninn,VDS n,recy= η VD H (V) + ηn VD Hττnwhere the time constant of the radial plasma diﬀusion isestimated in terms of the input radial decay length d1 andparticle diﬀusion coeﬃcient D⊥ as follows;
τ⊥ =

d12
D⊥

.

(25)

The eﬀect of non-zero radial heat diﬀusion coeﬃcients,
χi,⊥ and χe,⊥ , can be considered by setting the values of
the coeﬃcients g⊥ to be larger than 1. In this paper, weused g⊥i,k = g⊥i,⊥ = g⊥e = 3 assuming that χi,⊥ = χe,⊥ = 2D⊥ .
In this study, we used D⊥ = 0.5 m2 /s.

2.5

Virtual divertor model

Instead of imposing boundary conditions, such as Mt = 1,
sheath heat transmission factors for plasma and recycling
flux for neutrals, we use a virtual divertor (VD) model[11–14]. We introduce an artificial region (VD region) beyond the sheath entrance, as shown in Fig. 1, and set artificial sinks for plasma particle, momentum, energy, diﬀusion neutrals and an artificial source for recycling neutrals
there as follows;n
S VD = − VD ,
τ

(26)

VD
Mm
=−

(27)

(
)
mi nV
∂
VD ∂V
m
nD
,
+
i
m
∂s
∂s
τVD
(
)
1 1
2
VD 1
QVD
=
−
m
nV
+
g
nT
i
i,k ,
i,k
i,k
2
τVD 2
QVD
i,⊥ = −gVDi,⊥ nT i,⊥

,
τVD
(
)
1
VD 3
QVD
=
−
g
nT
e ,
e
τVD e 2

(28)

(29)
(30)

VD
S n,diﬀ=−

nn,diﬀ
.
τVD
n

(

)
∂nn,diﬀ, (32)−∂x(33)

Here, the definitions of the notations are the same as Ref.
[12]. We used τVD = 5 µs and τVD= 1 µs. The energynsink coeﬃcients were set to be gVD=
3, gVDi,⊥ = 1.2 andi,kVDge = 2.5. Strictly speaking, we need to control the values
VD
of τVDin order to satisfy the boundary conditions,n and gEqs. (6), (7) and (11). Nevertheless, we used fixed values
for them in this study because the ranges of output valuesof γi , γe and Cn were suﬃciently narrow. The recycling
rates were η = ηn = 1.

3

Results

3.1

Supersonic flows in no-neutral conditionsFirst, we conducted simulations for diﬀerent values ofd1 /dDIV in no-neutral conditions in order to study the supersonic flows. Figure 2 shows the Mach number at the
sheath entrance Mt and at the X-point MX as functions of
C which is defined by the following equation;
RΓ
C≡
Rp

√

Tt
.
TX

(34)

Here, RΓ and R p represent the particle and momentum flux
amplification factors, respectively, and T ≡ T i,k + T e . The
theoretical values shown by the dashed lines are obtainedin the same way as Refs. [8, 13]. By setting finite values ofd1 /dDIV , supersonic flows are driven and they agree wellwith the theoretical values.

10smoother with radial plasma diﬀusion.
For the present study, the homogeneous magnetic field
and the symmetric inner/outer SOL-DIV systems are assumed for simplicity. In order to simulate the inhomogeneous magnetic field and the asymmetric SOL-DIV systems, the eﬀects of the mirror force and the plasma current should be considered. For the transient simulations,
we need to control the values of gVD and τVD
n in order tosatisfy the prescribed boundary conditions. These are our
future works.

(a)8

d1/d = 0.5

6

d1/d = 1RΓ

d1/d = ∞d1/d = 0.3

4
2
0
1.0

References(b)0.8
d1/d = 0.3
Rp

0.6
d1/d = 0.5

0.4
0.2

d1/d = ∞d1/d = 1

0.0
0

1

2
3
4
22
Γsep (10 /s)

5

6

Fig. 3 Particle and momentum flux amplification factors R as
functions of Γsep for d1 /dDIV = ∞, 1, 0.5 and 0.3.

3.2

Eﬀect of radial plasma diﬀusion on detached plasmas

We investigated the eﬀect of the radial plasma diﬀusion on
the relations between the flux amplification factors R andthe particle flux from the core plasma Γsep . Here, the geometrical parameters are the same as Ref. [12] and the heat
flux from the core plasma is set to be Psep = 4 MW. The
impurity ratio is rimp = 1 %. The result is shown in Fig.
3. When d1 /dDIV = 0.5, RΓ starts to decrease (roll-over)with smoothly decreasing R p , which is qualitatively similar to the experimental observations. In these simulations,supersonic flows are not observed.

4

Summary and Future Works

Radial plasma diﬀusion is considered to be important as a
candidate to reduce the particle flux to the divertor plate inaddition to the volume recombination. We investigated theeﬀects of radial plasma diﬀusion on the detached divertor plasmas by using a 1D SOL-DIV plasma fluid model
incorporating anisotropic ion temperature [11–13]. A VD
model is implemented in our code in order to reproduce theeﬀects of divertor plate and sheath [14]. The supersonicflows are driven by the eﬀect of the radial plasma diﬀusion which agrees well with the theoretical values. The
roll-over of the particle flux to the divertor plate becomes[1] R. Schneider et al., Contrib. Plasma Phys. 46 (2006)
3.
[2] H. Kawashima et al., Plasma Fusion Res. 1 (2006)
031.
[3] K. Shimizu et al., Nucl. Fusion 49 (2009) 065028.
[4] T. D. Rognlien et al., J. Nucl. Mater. 196-198 (1992)
347.
[5] A. V. Chankin et al., J. Nucl. Mater. 390-391 (2009)319.
[6] M. Wischmeier et al., J. Nucl. Mater. 415 (2011)
S523.
[7] K. Hoshino et al., J. Nucl. Mater. 463 (2015) 573.
[8] T. Takizuka et al., J. Nucl. Mater. 290-293 (2001)
753.
[9] O. Marchuk and M. Z. Tokar, Contrib. Plasma Phys.
48 (2008) 164.
[10] K. Hoshino et al., Plasma Fusion Res. Ser. 9 (2010)592.
[11] S. Togo et al., J. Nucl. Mater. 463 (2015) 502.
[12] S. Togo et al., Contrib. Plasma Phys. to be published.
[13] S. Togo et al., J. Comput. Phys. 310 (2016) 109.
[14] T. Takizuka et al., 13th Burning Plasma Simulation
Initiative (BPSI) Meeting (2016).
[15] R. V. Bravenec et al., Phys. Fluids 25 (1982) 608.
[16] E. Zawaideh et al., Phys. Fluids 29 (1986) 463.
[17] W. Fundamenski, Plasma Phys. Control. Fusion 47(2005) R163.
[18] A. K. Prinja et al., J. Nucl. Mater. 196-198 (1992)
340.
[19] D. E. Post, J. Nucl. Mater. 220-222 (1995) 143.

Predictive simulation of fusion plasma by data assimilation
Takahisa Ozeki
Fusion Research and Development Directorate
Japan Atomic Energy Agency
Introduction and motivation: Experiments for fusion energy have been carried out in many
tokamak facilities in the world, and experimental data related to the fusion plasma isaccumulated by these facilities. In ITER remote experimentation center (REC) based on the
broader approach (BA) activity of the joint program of Japan and Europe (EU), the storage
system is being prepared as one of the important function, which can store/access the hugeamount of fusion relevant experiment data. For the efficient use of the experiment data, theintroduction of statistics is expected, such as Bayesian analysis for the data analysis todetermination of the unknown parameter of Modified Rutherford Equation for thestabilization of NTM, etc., and Big data analyses using analysis engine, visualization, datacollection/mining, work flow, etc. For the simulation, data assimilation to the simulation canbe expected, that is the prediction of the fusion plasma by the theoretical simulation withexperiment data. Here, predictive simulation of fusion plasma by the data assimilation is
discussed.
In DEMO reactor, where possible diagnostics are very limited, the simulation is strongly
expected for the guide of the operation. The integrated simulation of the fusion plasma isuseful tool for the guide of the plasma operation by the integration of the plasma model. In
addition to the development of the plasma mode for the integration, the efficient use of thereal experimental data with the simulation is useful method to navigate the DEMO operation.
The assimilation of the experiment data can be obtained by adapting the experimentallyobserved data on numerical simulation of the fusion plasma. Other methods are to set the
optimal initial conditions by assimilation of past data or to embed and to adapt the observed
value data to the physical model of the simulation. It has already been performed in the fieldof Meteorology, such as weather forecasts1-2), but its application in the nuclear fusion field
has not been enough yet.
Methods of data assimilation: The simple methods of the data assimilation on the simulation
are the direct insertion that is to replace simulation results by observed data or the Nudging
that is to slightly bring simulation results to observed values. These methods are simple but
need many observed data point, otherwise the assimilation is not enough for accuracy. For
more accurate assimilation, methods of the Optimal Interpolation and 3-dimensionalvariational method are proposed as the 3-dimensional data assimilation, where the generalized
form can be obtained by using the Bayes’ theorem. For the method including the past data,the 4-dimensional data assimilation can be applied, such as Sequential data assimilation of
Kalman filter or Ensemble Kalman filter, and 4-dimensional variational method. Thesemethods can be understand as the minimization method of the difference of the value between
the simulation data and the observed data
Formulation for the data assimilation: The linear system is discussed here as the first step.
The discretized numerical simulation model can be described by xt=Ft xt-1 to estimate the stateof xn . The system model and the observed model can be written as follows respectively:
xt = Ft xt−1 + vt , yt = H t xt + wt
where, xt is variable vector of the system model at the time t , yt is variable vector of theobserved model at the time t , vt is uncertainty by the modelling and the boundary condition x0,

wt is observed noise. Dimension of xt is larger than that of yt. Therefore, Ht is written like as,
In this case, the dimension of the observed point, yt, is three that is smaller than the dimension" 1 0 0$
H =$ 0 1 0$ 0 0 1#

⋅ ⋅ 0 %'
⋅ ⋅ 0 '
⋅ ⋅ 0 '&of the simulation data, xt.
The most plausible estimated value, xt, est can be obtained in the case of maximizing p(xt | yt),which is written asp(yt | xt )p(xt ),p(xt | yt ) =∫ p(yt | xt )p(xt )dxt
the conditional probability that can be obtained using the Bayes’ theorem. Here, p(xt) is theprobability of “xt at the point of no observation would be around this value”. p (xt) and p(yt |xt) are assumed to follow a multi-dimensional (k order) normal distribution, such as follows:" 1%1p(xt ) =exp $ − (xt − xt,b )T V −1 (xt − xt,b )'k&( 2π ) | V | # 2

p(yt | xt ) =" 1%
exp $ − (yt − H t xt )T R −1 (yt − H t xt )'&( 2π ) | R | # 2
1

kusing V, R is the variance-covariance matrix (positive value symmetric matrix). Xt,b is themean vector. From Bayes’ theorem, for minimizing the following function11(xt − xt,b )T V −1 (xt − xt,b ) + (yt − H t xt )T R −1 (yt − H t xt ) ,22
xt becomes the most plausible estimated value , xt, est, as follows:
x t, est = xt, b +VH tT (H tVH tT + R)−1 (yt − H t xt, b )This means maximizing the probability of p(xt | yt). For the sequential data assimilation using
Kalman filter, the conditional probability p(xt | yt) is obtained by using the average data andthe variance-covariance matrix of xt-1 on previous step. Values of the next step, t|t-1, areobtained from value of the previous step, t-1. The averaged value of p(xt), xt|t-1 andvariance-covariance matrixVt|t-1 are obtained from the system model as follows:
xt|t−1 = Ft xt−1|t−1 , Vt|t−1 = FtVt−1|t−1FtT + Q t
The filtered value of xt|t-1 and p(xt | yt) are obtained by the filtering of . xt|t-1 and Vt|t-1 .
xt|t = xt|t−1 +Vt|t−1H tT (H tVt|t−1H tT + R t )−1 (yt − H t xt|t−1 )Vt|t = "# I −Vt|t−1H tT (H tVt|t−1H tT + R t )−1 H t $%Vt|t−1
Implementation of Kalman filter to the fusion plasma: The heat and particle transport of thefusion plasma using the well known transport equation of the ion and electron is proposed as
the implementation of the data assimilation simulation.
3∂1∂e+ Peq − Prad + Pne(neTe ) = − (rQe ) + POH + Pheat
2 ∂t
r ∂r
3∂
1∂
i
− Peq − Pcx + Pni(niTi ) = − (rQi ) + Pheat
2 ∂t
r ∂r
∂ Te 3∂T 3∂nhere, Qe = −ne xe+ Te Γ e , Qi = −ni xi i + Ti Γ i and Γ i = Γ e = −D i − niVin
∂r 2
∂r 2
∂r

At the first step, it is assumed that these transport equations can be linearized as the systemmodel,
xt=Ft xt-1 + vt .
3 ∂x
1 ∂
n =−
(rQ) + ⋅⋅⋅⋅⋅ , it can be imagine that xt can be the
2 ∂t
r ∂r
electron temperature, Te, or the ion temperature, Ti, at the data on the plasma minor radius. vtof the system model is produced by the uncertainty due to the transport model, such as χe , χi.
Variance-covariance matrix V of xt for the various uncertainty of N is

From the form of the equation,

1 N (i)(i)
Vt|t−1 =xt|t−1 − x̂t|t−1 ) ( xt|t−1− x̂t|t−1 )∑(N −1 i=1T.

1 N (i)∑ xt|t−1 For the observed data, yt, the diagnostic
N i=1
data from Thomson scattering, ECE, etc. can be applied. The observation noise can beobtained from the error and noise of the diagnostic system.
Here, the average value x̂t|t−1 are x̂t|t−1 =For the non-linear system model, Ensemble Kalman filter can be applied, where the systemmodel will be xt = ft (xt−1, vt ) . The procedure of the calculation is same manner in that of theN

Kalman filter. The system noise of N can be obtained {vt(i) } . The estimated value xt|t-1 is
i=1

calculated for each i as x

(i)
t|t−1

= ft

(i)

(x

(i)
t−1|t−1(i)t,v) . The observed noise is made from theNdiagnostic data for each i as {wt(i) } . Using the variance-covariance matrix V and R, xt|t-1 willi=1

be filtered, namely modified.
(i)(i)
xt|t(i) = xt|t−1+ K t ( yt − H t xt|t−1)

Here, Kt is Kalman gain. These steps will proceed for each i in parallel, namely theseprocesses are ensemble.
In Summary, the possibility of the implementation of the data assimilation to the simulation of
the fusion plasma is discussed. Considering the difficulty with the treatment of the observed
data, it seems to be easiest way of the application to the transport equations, however it is still
not clear in the validity of the assumption on application to the Kalman filter, the handle thesystem noise, and so on. It is necessary to investigate of the other filtering models and toproceed the implementation to the actual fusion plasma model.
Acknowledgement:
I would like to thank for the discussion with Drs. O. Naito and S. Ide of JAEA.
References
1) K.Nakamura, G.Ueno, T.Higuch, Tokei-Suri, Vol.53, 2, 211-229 (2005)
2) K.Yumimoto, Nagare Vol. 28, 37-44 (2009)

Development of PHZ pinch modelwith plasma rotation effects on particle orbit
T. Fujita, Y. Shimizu
Graduate School of Engineering, Nagoya University
Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan

The inward pinch model for high-Z impurity ions in a toroidally rotating plasma, the PHZ pinchmodel, has been improved to take account of the velocity pitch angle distribution and effects of thecentrifugal force and the Coriolis force. It was found that the PHZ pinch velocity is large forparticles with larger orbit width. The PHZ pinch velocity averaged over the pitch angle becomeslarge with increase in the plasma rotation as in the original model.

1. Introduction
Related to the observation that accumulation oftungsten (W) was enhanced with increase in the
counter toroidal rotation in JT-60U H-mode
plasmas [1], two pinch models (PHZ pinch and Er
pinch) for high Z impurity ions caused by thetoroidal rotation and the radial electric field were
proposed [2]. The PHZ pinch is caused by change
in the ion charge number Z along its drift orbit. The
Er pinch is caused by collision in the deformed drift
orbit due to the radial electric field Er. These twopinch models have been introduced into TOTAL
and the results of the JT-60U experiments were
simulated. Since the assumption employed fordevelopment of the Er pinch velocity formula was
not satisfied, only the PHZ pinch was evaluated.
Dependence of W accumulation on the plasma
rotation was seen, but was weaker than theexperimental observation [3].
In the original PHZ pinch model, the situationthat the toroidal rotation speed of impurity ions,which is equal to that of plasma (bulk ions), ismuch larger than their thermal speed is assumed
and then only passing particles with velocity
constant during its orbit circulation is considered.
In the experiment, however, concentration of highZ impurity on the low-field-side (LFS) was
observed in rotating plasmas [4]. In the frame
rotating with the same speed as the plasma, theorigin of this phenomenon can be regarded as thecentrifugal force acting on the particles. In thisarticle, we report on the formulation for adding
this effect on the original PHZ model and show
some initial results obtained by the new formula.

2. PHZ for arbitrary pitch angles
In the original PHZ model, the following
assumptions were made.
(i) The particles have monotonic energy and its

velocity component perpendicular to themagnetic field 𝑣! is negligibly small compared
to its velocity component parallel to themagnetic field 𝑣∥ . So, conservation of magnetic
moment µ does not cause any decrease in 𝑣∥ on
the high-field side (HFS) and hence no trapped
particles exist.
(ii) The change in the particle energy caused by its
radial drift with presence of the radial electricfield Er is not taken into account while change inthe poloidal velocity due to the ExB drift istaken into account.
Concentration of heavy ions on the LFS
indicates existence of trapped heavy ions. So we
calculate particle guiding center orbit with arbitrary
pitch angles, namely including both passing andtrapped particles, and evaluate the PHZ pinch
velocity generated in the orbit
It is assumed that the W ions rotate with thesame speed as the bulk plasma (bulk ions). So, their
speed in the rotating frame is only due to thermal
motions and is isotropic. So, first we extend the
PHZ model for arbitrary particle pitch angles,
before including plasma rotation effects, namely thecentrifugal force and the Coriolis force.
A new code for calculating the PHZ pinch was
produced. We plan to introduce it into TOTAL codein future. The guiding center orbit is numericallycalculated with the 4th order Runge-Kutta method
in a circular cross-section equilibrium for a particle
with monotonic energy and with an arbitrary pitch
angle. In this article, 𝜉   =    𝑣∥   /𝑣 is used for the
pitch angle parameter. First, we calculate the orbit
with constant charge Z of ions. Then, we calculate
variation of Z along the orbit, which is caused by
the temperature variation along the drift orbit andevaluate the PHZ pinch velocity 𝑉!"# by
integrating 𝑣! during orbit circulation. This isbasically the same procedure as employed in theoriginal PHZ model, though it is now evaluatednumerically while it was evaluated analytically inthe original model [2].
1 !
𝑉!"# =
𝑣 𝑑𝑡
𝑇 ! !
𝑚!
1
𝑣! =
𝑣∥! + 𝑣!! sin 𝜃
𝑍𝑒𝑅! 𝐵!
2
Here T is the period for completing the closed
orbit in the poloidal cross-section, 𝑚! and Z are
the mass and charge number of impurity ions, 𝑅!
is the plasma major radius, 𝐵! is the magnetic
field at 𝑅 = 𝑅! , and 𝜃 is the poloidal angle.
We calculated the orbit in the laboratory frame,
namely without the centrifugal force or the
Coriolis force, but with arbitrary 𝜉 , and thencalculated 𝑉!"# . We compared 𝑉!"# for theparticle with 𝜉   = 1 and verified good agreementwith the original model [5]. Next the orbits and thePHZ pinch velocity for arbitrary pitch angle
particles were calculated. An example of the resultis shown in Fig. 1. Figure 1 (a) shows the half widthof the orbit ∆  =    𝑟!"# − 𝑟!"# /2 , where 𝑟!"#and 𝑟!"# are the maximum and minimum valuesof the minor radius in the orbit circulation. Here theboundary between trapped and passing particles isgiven by 𝜉 = 0.61. The ∆ is larger for trapped

particles except for deeply trapped ones, or near the
point 𝜉   = 0, than for passing particles. The 𝑉!"#is large for the particles with large orbit width, suchas barely trapped ones, as shown in Fig. 1 (b). The
PHZ pinch velocity for particles with a pitch angle
distribution would be evaluated by averaging 𝑉!"#given in Fig. 1 (b) with weight proportional to its
velocity distribution function, which is denoted by
𝑉!"# here. The 𝑉!"# becomes larger withincrease in the particle velocity or energy.
3. Centrifugal force and Coriolis force effecton PHZThe equation of motion for the guiding centersin the rotating frame is given in [6,7]. Here therotation speed of the frame is taken to be the sameas the plasma (bulk ion) rotation speed, assumingthat the rotation speed of the impurity ions is thesame as that of bulk ions. Furthermore, it isassumed that Er = 0 in the rotating frame. Note that
this assumption does not mean that Er = 0 in thelaboratory frame. The Er in the laboratory frame
𝐸!!"# and in the rotation frame 𝐸!!"# are related by
𝐸!!"# = 𝐸!!"# + 𝑅𝐵! Ωwhere   Ω is the angular frequency of the rotation.
Then, 𝐸!!"# > 0   in co-rotating plasma and𝐸!!"# < 0 in counter-rotating plasma when
𝐸!!"# = 0. This is usually the case in the tokamak
experiment.
With presence of the centrifugal force, the sumof the potential energy derived by the centrifugal
force and the kinetic energy is conserved [7]:1
𝑚 𝑣 ! − 𝑅 ! 𝛺∗!   ,2 !
𝑍𝑚! 𝑇!
Ω∗ = Ω 1 −𝑚! 𝑇! + 𝑇!
ℇ=

(1)
!/!

  ,

where 𝑚! denotes the mass of bulk ion. The Ω∗ isthe effective rotation angular frequency, or
modification of Ω considering the electric fieldalong the major radius direction generated by
non-uniform density distribution of bulk ions on theflux surface due to the centrifugal force.
The drift velocity is given by as follows, taking
into account effects of the centrifugal force and the
Coriolis force [6]:Fig. 1 (a) The half width of the guiding center
orbit Δ and (b) the PHZ pinch velocity 𝑉!"# asa function of the pitch angle parameter 𝜉 forthree values of particle velocity, 0.5×10! m/s(red solid line), 1×10! m/s (blue dotted line)
and 1.5×10! m/s (green dot-dash line).

𝑣! =𝑚
𝑍𝑒𝑅! 𝐵!

𝑣∥ ! +𝑣! !
+ 𝛺∗ ! 𝑅 ! + 2𝑣∥ 𝛺𝑅   .
2(2)

The velocity parallel to the magnetic field isdetermined by conservation of ℇ in eq. (1) in
addition to the conservation of magnetic moment µ.

done by the radial electric field on particles moving
along the minor radius changes the kinetic energyof particles as done in eq. (1). The vertical drift
velocity given by eq. (2) agrees with the usual drift
velocity if the effect of the electric field along the
major radius direction, generated by the
non-uniform distribution of bulk ions on the flux
surface due to the centrifugal force, is taken into
account.

Fig. 2. (a) The PHZ pinch velocity VPHZ as afunction of pitch angle parameter ξ, for toroidal
plasma rotation velocities of 0 m/s (red solidline), −0.3×10! m/s, −0.5×10! m/s,   −0.7×
10! m/s and  −1×10! m/s (purple dotted line).
(b) The PHZ pinch velocity averaged over the
pitch angle <VPHZ> (red solid line) and thefraction of trapped particles ftrap (blue dottedline) as a function of the toroidal plasma rotation
velocity. The symbols "A" to "E" denote thevalues observed in the JT-60U experiment andused in the simulation in [3].
The analysis was done for the W ions with thethermal speed (4.58×10! m/s) at 𝑟 𝑎 = 0.76 in
JT-60U plasmas shown in [1]. The PHZ pinch
velocity is shown in Fig. 2 (a) as a function of 𝜉for counter-rotating plasmas. The 𝑉!"# becomeslarger with increase in the counter toroidal
rotation speed. The fraction of trapped particles
also becomes larger with increase in the countertoroidal rotation speed as shown by the dottedline in Fig. 2 (b), which is caused by the effect of
the centrifugal force on ℇ and then on 𝑣∥ shownin eq. (1). The PHZ pinch velocity averaged over𝜉 becomes larger with increase in the countertoroidal rotation speed as shown in Fig. 2 (b). The
𝑉!"# changes by a factor of 5 within the rangeof the plasma rotation observed in the JT-60U
experiments as shown by symbols A to E in Fig. 2(b).
It should be noted that the equation of motion in
the laboratory frame would also lead to the same
results as those in this study where the equation of
motion in the rotating frame is employed. The work

4. Summary
The PHZ pinch model has been improved byintroducing the distribution of the velocity pitch
angle and effects of the bulk plasma rotation on theparticle orbit (the centrifugal force and the Coriolis
force). It was found that the PHZ pinch velocity islarge for particles with larger orbit width such as
the barely trapped particles. The PHZ pinch
velocity averaged over the pitch angle becomeslarger with increase in the plasma rotation velocity
as observed in the original PHZ model. We plan toextend the revised PHZ pinch model to include theradial electric field and then introduce it into the
TOTAL code to simulate transport of the high Z
impurity in rotating plasmas.
Acknowledgments
The authors are grateful to Drs. N. Hayashi, K.
Hoshino, T. Nakano and M. Honda for fruitful
discussion on the PHZ pinch model and on the
JT-60U experiment data. This work was carried out
under National Centralized Tokamak Collaboration
Research Program in JAEA. This work was
supported by JSPS KAKENHI Grand Number
25420895.
References[1] T. Nakano, et al., J. Nucl. Mater. 415,
S327-S333 (2011).
[2] K. Hoshino, et al., Nucl. Fusion 51, 083027(2011)[3] Y. Shimizu, et al., Plasma Fusion Res. 10,3403062 (2015).
[4] C. Angioni, et al., Nucl. Fusion 54, 083028(2014).
[5] Y. Shimizu, et al., submitted to Plasma Fusion
Res.
[6] K. L. Wong, et al., Phys. Fluids B 1 545 (1989).
[7] K. G. McClements, et al., Plasma Phys. Control.
Fusion 51, 115009 (2009).