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

第 14 回核燃焼プラズマ統合コード研究会
14th Burning Plasma Simulation Initiative (BPSI) Meeting(Ver.2)
日時：2016 年 12 月 7 日(水)－8 日(木)場所：九州大学筑紫キャンパス 応用力学研究所 2 階大会議室(18 min talk+7 min discuss or 10 min talk+5 min discuss)
12 月 7 日(水)9:00 – 12:45 核融合エネルギーフォーラムサブクラスター会合
12:45 – 13:45 昼休み
13:45 – 13:55 はじめに
村上 (京大)(座長：村上)13:55 – 14:20 講演 1-1
藤田 (名大)Analysis of axisymmetric equilibrium including magnetic field generated
by energetic particles
14:20 – 14:45 講演 1-2
本多 (量研機構)
Predictions of toroidal rotation and torque sources arising in
non-axisymmetric perturbed magnetic fields in tokamaks
14:45 – 15:10 講演 1-3
福山 (京大)Status of integrated modeling activity for ITER
15:10 – 15:35 講演 1-4
林 (量研機構)
Integrated modeling of JT-60SA high-beta steady-state scenario
15:35 – 15:50 休憩(座長：藤田)15:50 – 16:15 講演 1-5
T. S. Hahm (Seoul Univ.)
Role of trapped electron dynamics in ITB formation
沼波 (核融合研)
16:15 – 16:40 講演 1-6
Turbulent transport of multi-ion-species plasmas in helical systems
16:40 – 17:05 講演 1-7
糟谷 (九大)Turbulence diagnostic simulation for comparison with experiments
17:05 – 17:30 講演 1-8
奴賀 (核融合研)
Fokker-Planck simulations of runaway electron generation in tokamak
disruption1

17:30

散会

19:00 – 21:00 懇親会 (馳走 やまとや 博多駅東店にて)

12 月 8 日(木)9:00 – 9:05 事務連絡(座長：登田)9:05 – 9:30 講演 2-1
佐々木 (九大)Status and future plans of turbulence simulation on linear plasmas
9:30 – 9:55 講演 2-2
小菅 (九大)Cross scale energy transfer in parallel flow shear driven turbulence
9:55 – 10:10 講演 2-3
松井 (九大)Nonlinear simulation on profile formation with drift-wave instability inlinear devices
10:10 – 10:25 講演 2-4
阿部 (九大)Turbulence simulation of structural formation considering neutral particle
profiles in linear devices
10:25 – 10:40 休憩(座長：福山)10:40 – 11:05 講演 2-5
滝塚 (阪大)Kinetic simulation for edge plasma and detached divertor
11:05 – 11:30 講演 2-6
東郷 (筑波大)Research on behavior of divertor particle flux with a one-dimensional
numerical model and an analytical one
11:30 – 11:55 講演 2-7
大澤 (九大)Vacancy cluster growth in tungsten induced by hydrogen
11:55 – 13:00 昼休み(座長：糟谷)13:00 – 13:25 講演 3-1
矢木 (量研機構)
Simulation study on internal transport barrier formation using gryofluidmodel
13:25 – 13:50 講演 3-2
瀬戸 (量研機構)
Linear analysis of non-ideal ballooning mode instability with real electron
inertia2

13:50 – 14:05 講演 3-3
池田 (京大)PIC simulation of EC wave propagation and absorption in plasmas
14:05 – 14:20 講演 3-4
永井 (京大)Start-up simulation of burning plasma by integrated modeling code TASK
14:20 – 14:35 休憩(座長：矢木)
14:35 – 15:00 講演 3-5
登田 (核融合研)
Construction of reduced transport model by gyro- kinetic simulation withkinetic electrons in helical plasmas
15:00 – 15:25 講演 3-6
村上 (京大)Prediction of neutron production in the deuterium experiment of LHD
15:25 – 15:40 講演 3-7
前田 (京大)Heat and particle transport modeling of LHD plasma by integrated
transport simulation
15:40 – 16:00 まとめ糟谷 (九大)16:00 散会3

Analysis of axisymmetric equilibrium
including magnetic field generated by energetic particles
T. Fujita
Graduate School of Engineering, Nagoya University
Furo-cho, Chikusa-ku, Nagoya 464-8603, Japan

A module for calculating an axisymmetric equilibrium considering current by fast particles withorbits shifted from the flux surface has been developed in the TOTAL code. Initial results indicate
shift of the magnetic axis owing to the fast ion current.

1. Introduction
MHD equilibrium calculation is one of essential parts in an integrated transport code. In usual, the current
density profile is evaluated in the form of the flux surface average of inner product of the current density 𝒋and the magnetic field 𝑩, 𝒋 ∙ 𝑩 to be substituted to the Grad-Shafranov equation [1, 2].
∆∗ 𝜓! =

𝜕 ! 𝜓! 1 𝜕𝜓! 𝜕 ! 𝜓!
−+
= −𝜇! 𝑅𝑗!"
𝜕𝑅 !
𝑅 𝜕𝑅
𝜕𝑍 !

𝐵!!
𝑑𝑝
𝐹 𝑑𝐹
𝑑𝑝
𝑗!" 𝑅, 𝜓 = 𝑅
+
=𝑅
1− !
𝑑𝜓 𝜇! 𝑅 𝑑𝜓
𝑑𝜓
𝐵

(1)

+

𝒋∙𝑩
𝐵!
𝐵!

𝒋 ∙ 𝑩 = 𝒋𝑶𝑯 ∙ 𝑩 + 𝒋𝑩𝑺 ∙ 𝑩 + 𝒋𝑵𝑩 ∙ 𝑩 + 𝒋𝑹𝑭 ∙ 𝑩 + ⋯

(2)
(3)

where 𝑅, 𝜙, 𝑍 is the cylindrical coordinate system, 𝜓! is the ploidal flux function by the 𝜙 componentof the plasma current density 𝑗!" , 𝜓 = 𝜓! + 𝜓!" is the total ploidal flux function, 𝜓!" is the poloidal flux
function by the coil current, 𝑝 is the total plasma pressure including the fast particle pressure and𝐹 𝜓 = 𝑅𝐵! . The 𝒋𝑶𝑯 , 𝒋𝑩𝑺 , 𝒋𝑵𝑩 and 𝒋𝑹𝑭 are the inductive current density, bootstrap current density,
NB-driven current density and RF-driven current density, respectively.
This treatment seems to be incorrect for the current generated by energetic or fast particles whose
orbits are not on the flux surfaces. In current hole plasmas, the shift of ion orbit from the flux surface is close
to the order of the plasma minor radius [3]. In EC heated low-Ip plasmas, the electron orbit expands evenoutside of the last closed flux surface LCFS [4]. A self-consistent equilibrium including the fast particles was
obtained before [5]. The model, however, focused on the runaway electrons and then the gyro-motion ofenergetic particles was not included.
In this article, we report on the formulation for finding an axisymmetric self-consistent equilibrium
including the fast particle current and show some initial results.

2. Basic equations
Force balance equations in thermal or bulk plasma fluid and fast particle fluid are given by
𝟎 = 𝒋𝒑 ×𝑩 − 𝛁𝑝!!

(4)

𝑛! 𝒖𝒇 ∙ ∇ 𝑚! 𝒖𝒇 = 𝒋𝒇 ×𝑩 − 𝛁∥ 𝑝𝐟∥ − 𝛁! 𝑝𝐟!

(5)

where 𝒋𝒑 and 𝑝!! are the current density and the pressure of the thermal plasma fluid, 𝒖𝒇 is the
fast-particle fluid velocity and 𝑝𝐟∥ and 𝑝𝐟! are components parallel and perpendicular to 𝑩 of the pressureof the fast-particle fluid. In [4], monotonic energy beam without thermal pressure is assumed in the second
equation. Then 𝒖𝒇 can be replaced with the particle velocity 𝒗𝒇 and eq. (5) becomes𝑛! 𝒗𝒇 ∙ ∇ 𝑚! 𝒗𝒇 = 𝒋𝒇 ×𝑩.

(6)

From eq. (6) a partial differential equation is obtained for Uf, the stream function of 𝒗𝒇 , which is solved
simultaneously with the Grad-Shafranov equation.
In this study, we do not use the force balance in the fast particle fluid but trace the guiding-center orbit offast particles. The slowing down and the pitch angle scattering by Coulomb collision are neglected. The
guiding center orbit is calculated by using conservation of the 𝜙 component of canonical angularmomentum 𝑃! given by eq. (7), where 𝑚 and 𝑞! are the mass and charge of the particle and 𝑃!" is thevalue of 𝑃! at the birth point of the particle. The toroidal component of the particle velocity 𝑣! is given byeq. (8), where 𝑣∥ is the velocity parallel to 𝑩, 𝑣!" is the 𝜙 component of the drift velocity 𝒗𝒅 given byeq. (9). The 𝑣∥ and 𝑣! (the gyro-motion velocity perpendicular to 𝑩) are given by eqs. (10) and (11) byconservation of magnetic momentum and energy, where 𝑣 is the magnitude of the particle velocity, 𝐵! isthe magntiude of 𝑩 at the birth point and 𝜉   is the ratio of 𝑣∥ to 𝑣 at the birth point.
𝑃! = 𝑚𝑅𝑣! + 𝑞! 𝜓 = 𝑐𝑜𝑛𝑠𝑡. = 𝑃!"
𝑣! =
𝒗𝒅 =

𝐵!
𝐵

𝑣∥ + 𝑣!"

(7)
(8)

𝑚
𝑣!! 𝑩×𝛁𝐵
𝑣∥! +
𝑞! 𝐵
2
𝐵!

(9)

𝐵
1 − 𝜉!
𝐵!

(10)𝑣∥ = 𝑣 1 −𝑣!! = 𝑣 !

𝐵
1 − 𝜉!
𝐵!

(11)

Note that 𝜓 in eq. (7) contains the flux function by the fast particle current, 𝜓! . Namely𝜓 = 𝜓! +    𝜓! + 𝜓!"(12)

In this study it is assumed that no space electrostatic potential is generated by the fast particle even when
their orbit expands outside LCFS. The potential would be canceled by existence of SOL plasma.
We fix the birth point of fast particles and trace the contour of 𝑃! . The rectangular grind in the 𝑅𝑍
plane used for solving the Grad-Shafranov equation is used for tracing 𝑃! Once the 𝑃! contour is obtained,
the distribution of the toroidal current by the fast particle is calculated along the orbit if the total number ofparticles is given. In this study, the total toroidal current by fast particles 𝐼! is given instead of the total
number of particles.
The Grad-Shafranov equation from eq. (4) is the same as eqs. (1) and (2), except that 𝑝 is replaced with
pth in eq. (2). In an integrated transport code the parallel current or 𝒋 ∙ 𝑩 is obtained by the transportmodule. In this study, the toroidal current density is given by the following equation instead of eq. (2).
𝑗!" 𝑅, 𝜓 = 𝑅

𝑑𝑝!!
𝐹 𝑑𝐹
𝑅
𝑅!
+
= 𝑐!
+ 1 − 𝑐!
𝑔 𝜓𝑑𝜓
𝜇! 𝑅 𝑑𝜓𝑅!
𝑅(2')

where 𝑐! is a constant, 𝑅! is the plasma major radius, 𝑔 𝜓 is a function of 𝜓, and 𝜓 is the normalized
poloidal flux function such that 𝜓 = 0 at the magnetic axis and 𝜓 = 1  on the LCFS. The function 𝑔 𝜓 isdetermined such that the total thermal plasma current and the safety factor 𝑞 at the axis become the given
values.
A self-consistent solution of the Grad-Shafranov equation (1) and (2') and the orbit of fast particles (7)together with eq. (12) is obtained by iterating the following two steps. The first step is solving theGrad-Shafranov equation with the fixed fast particle current and then renewing 𝜓! and 𝜓!" . The second
step is following the fast particle orbit with the given 𝜓 and then renewing 𝜓! . This process is iterated until
𝜓! and 𝜓! are converged. This procedure has been introduced into the initial equilibrium construction partin the integrated transport code TOTAL.

Fig. 1 Comparison of tokamak equilibria with and without fast ion current; (a) with bulk plasma currentof 1 MA and no fast ion current (b) with bulk plasma current of 0.6 MA and fast ion current of 0.4 MA.
It is known that a return current by electrons is generated when a toroidal current is driven by fast ionscirculating the torus. This effect is not yet included in the present analysis.
3. Results
The equilibrium including the fast particle current has been calculated as shown in Fig. 1. Here, theequilibrium without the fast particle current and that with the fast particle current are compared. The total
toroidal current 𝐼! , sum of the bulk plasma current and the fast particle current is fixed to 1 MA. The shapeof the LCFS is also fixed. As the fast particle, deuterium ions with 1 MeV energy and with 𝜉 = 1 are
generated at 𝑅 = 3 m and 𝑍 = 0. The fast particle current 𝐼! was 0 MA or negligibly small in Fig. 1 (a)and was 0.4 MA in Fig. 1 (b). It is found that the radial width of the fast ion orbit shown by dotted lines issmaller and the magnetic axis is shifted toward the high-field side in Fig. 1 (b) compared to Fig. 1 (a).
The 𝑞 profiles are shown in Fig. 2. A nearly flat 𝑞 profile with weak magnetic shear was employed. A
localized negative shear region appears around 𝑟/𝑎   =   0.5 due to the localized fast ion current. Change inthe 𝑞 profile has possibility for affecting the shift of the magnetic axis shown in Fig. 1. It is needed toseparate the effect of 𝑞 profile change and that of current
by fast particles with orbits shifted from the flux surface.
The convergence became slower with increasing 𝐼! 𝐼! .
Tracing the particle orbit with fixed 𝜓 would be one of thecauses since it does not hold consistency between theparticle orbit and 𝜓! . Correction on this point would help toimprove the convergence.
Acknowledgments
This work was supported by JSPS KAKENHI Grand
Number 25420895.
References[1] S. Tokuda, et al., J. Phys. Soc. Japan 58, 871 (1989).
[2] K. Tani, et al., J. Comp. Phys. 98, 332 (1992).
[3] T. Fujita, Nucl. Fusion 50, 113001 (2010).
[4] T. Maekawa et al. Nucl. Fusion 45, 1439 (2005).
[5] Z. Yoshida, Nucl. Fusion 30, 317 (1990).

Fig. 2. The radial profile of safety factor 𝑞
of equilibria shown in Fig. 1. The dottedline denotes the case without the fast ion
current and the solid line denotes the casewith the fast ion current.

Predictions of toroidal rotation and torque sources arising in
non-axisymmetric perturbed magnetic fields in tokamaks
M. Honda, *S. Satake, *Y. Suzuki, K. Shinohara, M. Yoshida, §N. Aiba, J. Shiraishi, N.
Hayashi, G. Matsunaga, *M. Nakata, §A. Matsuyama and S. Ide
National Institutes for Quantum and Radiological Science and Technology, Naka,
Ibaraki 311-0193 Japan*National Institute for Fusion Science, Toki, Gifu 509-5292 Japan§

National Institutes for Quantum and Radiological Science and Technology, Rokkasho,
Aomori 039-3212 Japan

Integrated modeling on toroidal rotation for developing operation
scenarios
Toroidal rotation has attracted much attention in terms of not only the physics itself
but also its influence on confinement and transport [1]. Accurate predictions of toroidal
rotation are indispensable for those of plasma confinement. Recent advancement of thetoroidal rotation modeling in the integrated suite of codes, TOPICS, with the 3D
neoclassical transport solver, FORTEC-3D, and the 3D nonlinear equilibrium code,
VMEC, enables us to perform analyses of toroidal momentum transport and the radial
electric field 𝐸𝐸" , and predictions of toroidal rotation with torque sources. The NTV stems

from the non-axisymmetric perturbed magnetic fields in tokamaks and is an intrinsic

torque source, which is expected to be a major momentum source in ITER and a DEMO
reactor. Another intrinsic torque is the one due to the residual stress, which can spin up
rotation from rest. The residual stress generally depends upon the temperature gradient[2]; hence, it is imperative to include the reliable residual stress model for developing
high-performance operation scenarios with transport barriers.

The intrinsic torque due to perturbed magnetic fields
Thus far we have dealt with the n=18 symmetric non-resonant perturbed magnetic

field produced by the toroidal field (TF)
coils and the ferromagnetic inserts (FIs)[3,4], where n is the toroidal mode
number. However, ITER in reality havemore complex perturbed magnetic fieldcreated by the heating and diagnostic
ports and the TBMs. In particular, theFigure 1: (a) Vacuum perturbed magnetic field at(R,Z)=(8.186,0.476) versus the toroidal angle. (b)Displacement of the LCFS of the equilibrium.

TBMs give rise to large ripple amplitude of about 1% locally in the toroidal direction,
which potentially enhances the loss of fast ions and transport of bulk ions. Eachphenomenon is directly linked to the neutral beam injection (NBI) and NTV torques,
respectively. The quantitative investigation of the impacts of the realistic perturbed
magnetic fields on the estimate of the torque and toroidal rotation is of quite importance.
Effects of the non-axisymmetric perturbed magnetic fields on confinement of fast ionsproduced by NBIs and fusion reactions and on torque sources for the ITER 15MA DT Lmode scenario are investigated using the magnetic field produced by TF coils, FIs, ports
and TBMs. The plasma response is considered in the equilibrium with the displacementof the last closed flux surface (LCFS) up to 1% of the minor radius. Figure 1 shows the𝛿𝛿𝛿𝛿 amplitude and the plasma displacement 𝛿𝛿𝑅𝑅, where one can see the existence of thethree pairs of TBMs and the slight amplitude enhancement due to the ports. For reference,the ‘TF+FI’ case without effects of the TBMs and ports is also superimposed. The OFMC
analysis shows that fast ions produced by 1MeV D-NBIs are never lost. Also, alphaparticle confinement was investigated, indicatingthat 0.007% alphas merely ran away especially onto
the divertor plates and no deterioration of alphaconfinement due to TBMs was observed. Theseenergetic particles are mostly passing ones and areinsensitive to 𝛿𝛿𝛿𝛿. These results are the beneficial

findings that the TBMs and ports are unlikely todeteriorate fast-ion confinement and yield lossrelated torques. In contrast, thermal ions are

Figure 2: Profiles of torque density
due to NBI, NTV and alphas.

significantly affected, as seen from the torque profiles in Fig. 2. Since the NTV increases
towards the edge as the volume does, its torque becomes -23.5Nm, which is quite larger
than that for the ‘TF+FI’ case, and is even comparable to the NBI torque of 30.3Nm.
Applying the NTV decelerates core toroidal rotation from 180km/s to 145km/s.
Simulations elucidate the essential importance of the NTV regarding rotation predictions
in ITER.

The intrinsic torque due to the residual stress
The major unknown in the context of the rotation
modeling is the residual stress. Calibrating thetheory-based model [2] that can capture theexperimental tendency is a practical way to attainrotation predictive capability for high-performance
plasmas. A numerical factor C whose sign depends
upon the dominant instability is added to the modeland is examined with JT-60U data. For these shots
possible torque sources have in advance been

Figure 3: Improvement of the VϕC
reproducibility using the residual
stress model with C=0.3.

estimated, and the linear stability has been computed by the gyrokinetic code GKV. In
Fig. 3, it is found that reproducibility of the carbon toroidal rotation 𝑉𝑉'( is to large extent
improved, from dashed line to solid one, by the residual stress model with C=0.3 for thisH-mode shot, where the ion temperature gradient (ITG) mode predominates in the edge
region. It suggests a positive C, consistent with the C sign adopted here. The range of thisfactor is found to be -1≤C≤1 based on the several shots, a fact which can limit the rangeof predictions of toroidal rotation.

References[1] Mantica P et al., Phys. Rev. Lett. 107, 135004 (2011).
[2] Kosuga Y et al., Phys. Plasmas 17, 102313 (2010).
[3] Honda M et al., Nucl. Fusion 54, 114005 (2014).
[4] Honda M et al., Nucl. Fusion 55, 073033 (2015).

Integrated modeling of JT-60SA high-beta steady-state scenario
N. Hayashi, K. Hoshino, K. Shimizu, M. Honda, S. Ide
QST
The JT-60SA [1] mission is to contribute to early realization of fusion energy by supporting the
exploitation of ITER and by complementing ITER in resolving key issues for DEMO reactors.
Especially, JT-60SA explores the high-beta and steady-state operation regime to complement ITER.
For the purpose, JT-60SA is equipped with flexibly applicable actuators of NB and EC. The
scenario will use lower plasma current than the maximum one (5.5 MA), the plasma density needs
to be lower and thus the reduction of divertor heat load is a critical issue. Impurity seeding to thedivertor region can reduce the heat load below a preferable level of 10 MW/m2, however, theaccumulation of seeded impurity in the core may reduce the core confinement by enhancing theradiation and dilution.
In this paper, we predict the JT-60SA plasma by integrated modeling code TOPICS [2] which
includes a divertor integrated code SONIC and a core impurity transport code IMPACT. For theprediction, we validate typical anomalous heat transport models, which are the major uncertainty inthe prediction, for JT-60U and JET experiments with internal transport barriers in the Japan-EU
collaborative activity. In the models, CDBM agrees withthe experiments better than the other models, orunderestimates temperatures. Thus, we consider that
CDBM can be used for the conservative prediction.
TOPICS with CDBM is used for the prediction without
impurity seeding. A stationary state is obtained at the flat
top phase of plasma current of Ip=2.3 MA with Bt=1.7 T.
Various states with high-βN (>3.5) and nearly full CD
condition can be obtained by using various sets ofactuators. By using the total power of 24 MW, a plasma
with βN=3.9 and nearly full CD condition is obtained. The
scenario requires the low electron density at the separatrix(ne,sep < 1.7 × 1019 m-3), leading to the high heat load on
divertor plates. SONIC simulations are carried out to find asolution satisfying the conditions by Ar seeding. With the

FIG. Profiles of (a) Ar densities withcharge states from 15+ to 18+, (b)neoclassical convective velocities, (c)electron temperature and radiation
power density in a case only with
Carbon radiation (broken lines) and acase with seeded Ar (solid lines).

Ar puff amount of 0.17 Pa m3/s, a solution with ne,sep = 1.6 × 1019 m-3 and qdiv,out = 9.7 MW/m2 isobtained. TOPICS with IMPACT can clarify the accumulation of seeded Ar in the core and its
effect on the plasma performance. Figure shows profiles of Ar densities, neoclassical convective
velocities, electron temperature and radiation power density, where Ar densities at the edge
evaluated by SONIC is used as boundary conditions in IMPACT to obtain the maximum
accumulation. Considering the maximum accumulation, impurity anomalous diffusivities are set tothe minimum level of neoclassical ones. The inward neoclassical pinch due to the bulk density
gradient exists in the core and the outward convection due to the temperature gradient in thepedestal. This neoclassical convection results in the accumulation of Ar16-18+ in the core. The Arradiation power increases by 3.8 MW and thus the temperature decreases in the core. The Ar
accumulation reduces βN to 3.6 so that the full CD condition is no longer satisfied. However, it canbe recovered by adding the NB power to supplement the radiation increase.
Acknowledgments
This work was supported by a Grant-in-Aid for Scientific Research (C) (No 26420862) from theJapan Society for the Promotion of Science and was carried out using the HELIOS supercomputerat International Fusion Energy Research Centre, Aomori, Japan, under the Broader Approachcollaboration between Euratom and Japan.
References[1] JT-60SA Research Plan - Research Objectives and Strategy Version 3.3 2016, March,	 http://www.jt60sa.org/pdfs/JT-60SA_Res_Plan.pdf[2] Hayashi N. et al 2010 Phys. Plasmas 17 056112

実験との比較のための乱流計測シミュレーション
Turbulence diagnostic simulation for comparison with experiments
N. Kasuya and M. Nunami 1
Research Institute for Applied Mechanics, Kyushu University, Kasuga, Fukuoka 816-8580,
Japan
1
National Institute for Fusion Science, Toki, Gifu 509-5292, Japan
E-mail: kasuya@riam.kyushu-u.ac.jp
Improvement of performance of helical plasmas is remarkable in these days. Consideration of
the 3-dimensional (3-D) magnetic configurations is necessary for understanding their
turbulent transport, but it is not simple to identify the observation region in the actual 3-D
configurations, so numerical simulations can be used for the identification. Numericalsimulations can be utilized for the synthetic analysis to simulate the experimental
measurement [1,2] and the validation for the simulation result [3,4]. Many turbulencesimulation codes have been developed for calculations in experimental magnetic
configurations. For a helical plasma, using the GKV-X, which is a gyrokinetic Vlasov codewith 3-D equilibrium, quantitative comparisons with experiments have been carried out [4].
Turbulence Diagnostic Simulator (TDS) [5] is a combination of turbulence codes,
measurement modules and analysis routines, to carry out numerical experiments of plasma
turbulence, which can be utilized as a platform of the data analysis. A numerical diagnostic tosimulate an experimental measurement is carried out using gyrokinetic simulations of
Ion-Temperature-Gradient (ITG) modes in the LHD configuration. The effect of line
integration along the line of the sight (LS) is considered as in the Phase Contrast Imaging(PCI) measurement, and the spatial resolution for reconstruction of the local spectrum isevaluated.
Flux-tube simulations of ITG modes in the LHD configuration have been carried out with the
GKV-X code [4]. The GKV-X is a nonlinear gyrokinetic Vlasov flux-tube code to perform
turbulent transport simulations. The field aligned coordinate in the real space is used for
calculations, and geometrical information of the flux surfaces is constructed by using VMEC
code [6]. The adiabatic condition is used, so the electron density and the electrostatic potential
is connected with the Boltzmann relation. Quantitative comparison of the thermal diffusivity
coefficient has shown the good agreement with experiments [4], and other turbulent
characteristics, such as a fluctuation spectrum, should be compared for validation.
We develop an analysis routine for the GKV-X data taking into account of a LS ofexperimental diagnostics, such as PCI. The PCI is used for the measurement of density
fluctuation with high temporal and wavenumber k resolutions [7]. The signals are given as theintegral of the density fluctuation along the LS. The GKV-X code uses the magnetic surface
coordinate, so interpolation is necessary for evaluating the physical quantities in the real space.
It must be noted that the obtained signal includes components perpendicular to the LS, so is
usually combination of the radial and poloidal components, depending on the selection of theLS. To obtain k spectrum a finite region around the LS are picked out, and the quantities in it
are evaluated with a finite spatial resolution, corresponding to the experimental diagnostic.
The width and mesh size of the region determine the resolvable k range (resolution in the kspace).

The GKV-X is the flux-tube code, so a set of the data is localized near one flux surface with afinite radial width. Figure 1 shows three sets of the density fluctuation data with different flux
tubes ( = 0.43, 0.65, 0.83). The simulation parameters are written in Ref. [4]. Figure 1 also
shows the LS for the numerical diagnostic, which viewsnear the magnetic axis and is perpendicular to the
mid-plane. The density patterns are calculated, and the kspectrum is obtained with Fourier transform with the
Hanning window function. The line integration andtemporal average gives the 2-D k spectrum from eachflux-tube data. Here, the main components in the x and y
directions are the toroidal and poloidal components,
respectively. As the observation, all the density along theLS should be integrated, so every spectrum data issummed up with the weight of the mean density at theposition. ITG modes have longer wavelengths in the
direction of the magnetic field, so there are characteristic
FIG 1. Fluctuation profiles fromkx - ky relation in accordance with the magnetic field
GKV-X gyrokinetic simulation
direction. Note that components with finite k along theat  = 0.43, 0.65, 0.83. The LS
LS become weak by the line integration, which reducesis also indicated by the graythe contained components in the obtained signal.
bold line.
As the obtained signal from PCI is the line integrated one, so there is a problem to resolve thelocal values from the signal. The pitch angle of the magnetic field is used to help theidentification [7]. The method to resolve the local
spectrum from the integrated signal (2-D PCI) is
applied using the simulation data. Assuming k// = 0,
the typical k direction can be set to be perpendicular
to the magnetic field direction, which changes inaccordance with the vertical position in helicalplasmas. A vertical profile of the |k| spectrum canbe obtained by extracting the perpendicularcomponents from the integrated signal (Fig. 2). The
target fluctuations exist at the limited radial range,
but the reconstructed profile spreads in the wider
region. Larger part of the leakage comes fromdiscreteness of the data with the finite spatial
resolution. The discrete Fourier transform cannot
avoid the leakage, so the resolution is limited to themesh size in the k space. The variation of themagnetic field direction in the region to calculate
the spectrum gives other cause for the errors, whose
FIG 2. (a) Vertical profile of the kcontribution is comparable to the former one.
spectrum decomposed by using theBroadened fluctuation spectrum in local positions ismagnetic shear method. (b) Spectra
also one of the causes for the deterioration. For theat z, where the fluctuation data exists,
decomposition using the magnetic shear, the mesh
are also plotted.

is sparser with smaller k, so the spatial resolution is worse with smaller k: z ~ 0.2 at ki = 0.1,
z ~ 0.1 at ki = 0.4 in the case at z = 0.6 in Fig. 2, where i is the ion Larmor radius. The
local spectrum must be discussed with careful consideration of the spatial resolution.
The analysis routine can give a fluctuation pattern at an arbitrary position, and comparison
between the local fluctuation and the decomposition of the integrated signal gives
understanding of the observed signals. The ITG modes have a characteristic wavelength andfrequency, and difference in the spectrum can be distinguished at different radial positions,considering the spatial resolution. The propagation direction can be also deduced from thisspectrum. This is the rather simple case to observe the k components, but other LSs givecombination of the kr and k components, which makes the spectrum complex forinterpretation.
In this way, the numerical diagnostics have been carried out in the 3-D torus magnetic
configurations. For the 2-D PCI, one way to improve the spatial resolution is to apply theother conversion method, such as the maximum entropy method [7]. Considering the roughspatial resolution, quantities distant enough from each other can be identical, so thediagnostics will be applied for evaluation of the up-down asymmetry and detection of
long-range correlation in future. Variety of the numerical diagnostics from several view-points
can give physical understanding and quantitative comparison of turbulent plasmas.
Acknowledgements
Authors acknowledge discussions with Prof. K. Itoh, Prof. S.-I. Itoh, Prof. S. Inagaki, Prof. K.
Tanaka, Prof. K. Ida, Dr. S. Toda, and Mr. K. Kawadu. This work is supported by JSPS
KAKENHI Grant Number JP24760703 and JP16K06938, by the collaboration program of
NIFS (NIFS15KNST082, NIFS16KNXN323, NIFS13KOCT001) and of RIAM of Kyushu
University.
References[1] LIN, L., et al., Phys. Plasmas 16 (2009) 012502.
[2] HOLLAND, C. et al., Phys. Plasmas 16 (2009) 052301.
[3] RHODES, T.L., et al., Nucl. Fusion 51 (2011) 063022.
[4] NUNAMI, M., et al., Phys. Plasmas 19 (2012) 042504.
[5] KASUYA, N., et al., Plasma Fusion Res. 8 (2013) 2403070.
[6] SUZUKI, Y., et al., Nucl. Fusion 46 (2006) L19.
[7] TANAKA, K., et al., Rev. Sci. Instrum. 79 (2008) 10E702.

Fokker-Planck simulations of runaway electron generation
in tokamak disruption
H. NUGA1 , M. YAGI2 , A. FUKUYAMA3
1

2

National Institute for Fusion Science
National Institutes for Quantum and Radiological Science and Technology
3
Department of Nuclear Engineering, Kyoto University

Introduction Disruption is one of the most serious events in tokamak discharges, since it induces
huge electromagnetic force to the device and generates high-energy runaway electrons (REs) which
may cause the damage of the plasma facing components[1]. Especially in large size tokamak, suchas ITER, since there is a possibility that the RE current can reach to a few mega-amperes andcause intolerable heat load to the first wall, the damages should be avoided or suppressed. To avoid
the disruption and related damages, the fast plasma shutdown by means of Massive Gas Injection(MGI) has been proposed and demonstrated in several devices. It is expected that the high electron
density plasma achieved by MGI may reduce the RE generation owing to its high collisionality.
MGI for mitigating disruption, however, tends to yield the shorter thermal quench duration andmay enhance the primary RE generation rate through the so-called “hot-tail eﬀect”[2, 3].
Aim of this paper is the investigation of the RE current generation in the MGI mitigatingplasma including hot-tail (H-T) eﬀect. To include the H-T eﬀect, the primary RE generation rate isevaluated by the electron momentum distribution function obtained by using Fokker-Planck code
TASK/FP. Evolutions of RE generation and induced electric field are calculated self-consistently.
Additionally, the drop of the plasma temperature and the increase of the electron density and theeﬀective charge owing to MGI are implemented simply.
Models In the present research, we use Fokker-Planck code TASK/FP to obtain the primary RE
generation rate and the induced toroidal electric field self-consistently. The detailed expression of
the models of the primary and secondary RE generation rates, the induced electric field, and thethermal quench are displayed in [4].
Additionally, the evolutions of electron and ions densities and the eﬀective charge assuming
MGI are modeled. Usually, noble gases, such as neon or argon, are used for MGI. In the lowtemperature plasma (∼ tens of eV), however, most of impurity ions are not fully ionized. Since
the issue of ionization state of impurities is also complex, to simplify our model, we introducefpost-MGI eﬀective charge Zeﬀ
and virtual impurity ion species, which has no electron-bound state.
fThey satisfy the following relations: nef = n0e + Zi nif , Zeﬀ= (nD + Zi2 nif )/nef , where superscript 0
and f denote pre- and post-MGI density or eﬀfective charge, respectively. We assume that electron
density evolves as follows:ne (ρ) = (nef (ρ) − n0e (ρ))= nef (ρ)

t∆tMGI(∆tMGI ≤ t),1+ n0e (ρ)(t ≤ ∆tMGI )

(1)
(2)

(a)

R = 6.2m, a = 2.0m, b = 2.4m
15 MA
T 0 (0) = 20 keV, T 0 (1) = 2 keV
T f (0) = 10 eV
τq = 1 ms∆tMGI = 5ms

RE current [MA](b)n0e (0) = 1.0 × 1020 m−3n0e (1) = 1.0 × 1019 m−3
f0
Zeﬀ= 1, Zeﬀ=3

Current [MA]solid: including Hot-taildotted: excluding Hot-tail9
RE+ohmic current

3

RE current

0

0
02

10

-02

10

primary RE current
-06

10

-10

2

4

6

20 8
nfe [10 m ]1012mary, and secondary RE current for τq = 1 ms,
f∆tMGI = 5ms, and Zeﬀ= 3.

(a) 15

6

total RE current
secondary RE current

4

Figure 2: nef (0) dependence of the total, pri-

Deuteron12

solid: including Hot-taildotted: excluding Hot-tail8

10

Current dens. [MA/m2 ]Radii
Initial current
Initial temperature
Post-quench temp.
Thermal quench time
MGI duration
Pre-quench density
on axison edgeeﬀective charge
Ion species

12

RE current [MA]Figure 1: Plasma parameters8

pre-quench ohmic current
solid: RE current with hot-taildotted: RE current without hot-tail

6
4
2
0
0

0.2

0.4

ρ

0.6

0.8

1

20
16

10

(d)
ρ=0.0

Figure 4: Current density profile for nef (0) = 12 ×secondary gen. rate

1020 m−3 .

12

10

12

10
primary gen. rate

810

primary RE generation rate [m-3s-1]

RE generation rate [m-3s-1]

10

20

10

(e)

16

10

ρ=0.3

secondary gen. rate

12

10

primary gen. rate

8

10

0

5

10

15
time[ms]

20

25

30

blue solid: Δt=5ms
green dotted: Δt=3ms

8

10

4

10

0

10
12
10
8

10

4

10

ρ=0.3

0

10

Figure 3: Evolutions of (a): Total and RE current,

ρ=0.0

0

2

4 time[ms] 6810(b): primary and secondary RE generation rates onaxis and (c) on ρ = 0.3.
Figure 5: Primary RE generation rate on axis andρ = 0.3 for ∆tMGI = 5 ms and 3 ms.

where ∆tMGI is the duration in which the electron density increases. Here, impurity transport isomitted for simplicity.
Results In the following calculation, we chose parameters tabulated in Table 1.
Figure 2 shows the nef (0) dependence of the (a) total, secondary, and (b) primary RE current. It
is found that the density increase achieves to suppress the RE generation as expected. In low density
region (nef (0) = 2 × 1020 m−3 ), most of pre-quench ohmic current is converted into RE current (total
RE current: IRE = 12.7 MA). The primary RE current excluding H-T eﬀect is Iprim = 3.96 MA andthat including H-T eﬀect increases to Iprim = 4.73 MA.
Subsequently, we focus on high electron density case (nef (0) = 1.2 × 1021 m−3 ). Although thehigh electron density succeeds to suppress the RE generation, the presence of the H-T eﬀect makes
the total RE current twice in the region (IRE = 1.09 and 2.12 MA for excluding and including H-T
eﬀect cases). In this case, H-T eﬀect enhances the primary RE current about four times (Iprim = 3.3
mA to 14 mA), though the magnitude of the increment (∼ 11 mA) is very small compared to total
2

RE current. This small increment of the primary RE current is multiplied by the avalanche eﬀect
to 1.03 MA and total RE current becomes twice. Figures 3 show the evolutions of (a) total and RE
current, and (b) primary and secondary RE generation rates on axis and (c) ρ = 0.3. From fig. 3 (b)and (c), it is found that the H-T eﬀect makes a additional peak of the primary RE generation rate(t ∼ 5 ms). This is because, H-T electrons can be REs easily rather than bulk electrons even withthe weak electric field due to the low collisionality. Owing to the presence of earlier seed REs, thesecondary RE generation is also triggered earlier. On ρ = 0.3, since the primary RE generation rate
at first additional peak (t ∼ 5 ms) is suﬃciently greater than that at subsequent peak, the secondary
RE current density also becomes greater. Figure 4 shows the initial ohmic and RE current density
profiles for nef (0) = 1.2 × 1021 m−3 case. The RE current density profiles have a peak on axis rather
than the initial profile. Furthermore, the H-T eﬀect makes the peak lower and broadens the RE
current profile.
This result shows that the presence of the first additional peak of (t ∼ 5 ms) multiplies the total
RE current. Since ∆tMGI = 5 ms is chosen in this result, the RE generation from H-T electrons
begins before the electron density reaches to post-MGI density nef . Therefore, ∆tMGI = 3 ms ischosen to investigate the behavior of the H-T RE peak.
Figure 5 shows the primary RE generation rate for ∆tMGI = 5 ms and 3 ms cases. It is foundthat the faster increase of the electron density achieves to suppress the H-T RE generation (t ∼ 5ms) both on axis and ρ = 0.3. Consequently total RE current is also suppressed from 2.12MA(∆tMGI = 5ms) to 1.10 MA (∆tMGI = 3ms), which nearly equals to the result excluding H-T eﬀect(1.09MA, ∆tMGI = 3ms).
Summary The present simulation study makes clear that the H-T eﬀect is important for the RE
generation with short thermal quench even in high electron density plasma. We have confirmedthat the RE generation is suppressed with the increase of the electron density. If the H-T eﬀect
is not considered, the RE generation is reduced to ∼ 1 MA in the high density region (nef (0) =1.2 × 1021 m−3 ). On the other hand, if the H-T eﬀect is included, the eﬀect multiplies the primary
RE current slightly (∼ 10 mA) and, consequently, total RE current is doubled due to avalanche
eﬀect. This multiplication comes from the fact that H-T electrons can be REs earlier than bulk
electrons owing to its low collisionality. If the electron density increases up insuﬃciently until
the H-T RE generation starts, the H-T eﬀect multiplies the primary RE current. Conversely, it
is also confirmed that suﬃciently fast increase of the electron density may suppress the H-T RE
generation. Because of these reasons, the reliable estimation of the thermal quench, increase of theelectron density, and impurity transport are inevitable.

[1] T.C. Hender et al., Nuclear Fusion, 47, S128 (2007)[2] H. M. Smith and E. Verwichte, Phys. Plasmas, 15, 072502, (2008).
[3] T. Fehér, et. al., Plasma Phys. Control. Fusion, 53, 035014, (2011).
[4] H. Nuga, et at., Phys. Plasmas, 23, 062506, (2016).

3

Status and future plans of turbulence simulation on
linear plasmas
M. Sasaki1, 2
1Research Institute for Applied Mechanics, Kyushu University, Kasuga 816-8580, Japan

2 Research Center for Plasma Turbulence, Kyushu University, Kasuga 816-8580, Japan

1

Introduction

Importance of the parallel flows has been recognized in magnetized plasmas. Intrinsic toroidal
rotation in toroidal plasmas and the flows around cusp region of the earth’s magnetic field [1, 2].
In the case when the shear of the parallel flow becomes large, the flow shear instability, which iscalled D’Angelo mode, occurs [3, 4]. In a linearly magnetized plasma device, PANTA, this shear
flow instability was observed, and the density peaking by the D’Angelo mode was observed.
The relation between the particle and momentum transport was investigated [5], and the spatial
structure of the D’Angelo mode was identified [6]. In this way, the experimental studies havebeen developed. The turbulence simulation on the D’Angelo mode is required to accelarate thefundamental nonlinear processes of the shear flow instability.

2

Model

In this study, we develop a turbulence simulation code to simulate the D’Angelo mode. Basedon the Numerical Linear Device (NLD) code [7], the ion parallel flow is included. The model is areduced fluid model to calculate the time evolution of the density, electrostatic potential and theparallel ion flow. This mode is an extension of Hasegawa-Wakatani model with coupling of theion parallel flow, which can simulate the resistive drift wave and the D’Angelo mode. The meandensity and the ion parallel flow are introduced by giving the particle and parallel momentum
sources, respectively.

3

Results

The bifurcation of the turbulence states between the drift wave dominant state and the D’Angelo
mode dominant state by changing the strength of the density source with the shot-by-shotmanner, where the parallel momentum source is kept constant. The mode spectrums in eachturbulence state are shown in Fig. 1. The change between the two states is abrupt, which has
not been observed in experiments. The drift wave dominant case is obtained in the case that thedensity source is large. The symmetry parallel mode spectrum is obtained in this case, which isconsistent with the previous studies as in [7]. The drift waves forms the streamer in this case.
When the density source is small and the density gradient becomes small, the D’Angelo mode is1

excited by the parallel flow shear. In this case, the asymmetry of the parallel mode spectrum isobtained. The higher harmonics are excited and the nonlinear structure with (m, n) = (0, n) was
obtained. The spectrum asymmetry and the spatial structure of the eigenfunction are consistent
with the theoretical prediction [4] and the experiments [6]. The roles of each turbulence on thedensity and the parallel flow are also studied by evaluating the density and parallel momentum
flux. In the case of the drift wave dominant state, the density fluctuation flattens the density
profile, and the parallel momentum flux is relatively weak due to the fact that the parallel mode
spectrum is almost symmetry. In the D’Angelo mode dominant state, the fluctuation with theparallel mode spectrum asymmetry effectively drives the parallel momentum flux, which works
to flatten the profile of the parallel flow. For the density flux, we obtain the inward particle
flux, which steepens the density profile. These characteristics are found to be similar to theexperiment [5].
10

−4

10

r/a=0.2
8

8

−5

6

−5

m

m

6
−6

4

−6
4

−7

2
0

−4
r/a=0.2

−7

2

−5

0
n

5

0

−8

−5

0
n

5

−8

Fig. 1: Mode spectrums of the drift wave dominant state (left) and the D’Angelo mode dominantstate (right). Here m and n are the azimuthal and parallel mode numbers, respectively. The
spectrums are calculated from the potential fluctuations at r/a = 0.2.

References[1] P. H. Diamond, et. al., Nucl. Fusion, 49, 045002 (2009).
[2] H. Hasegawa, et. al., Nature, 430, 755 (2004).
[3] N. D’Angelo, Phys. Fluids, 8, 1748 (1965).
[4] Y. Kosuga, et. al., Plasma Fusion Res. 10, 3401024 (2015).
[5] S. Inagaki, et. al., Scientific Reports, 6, 22189 (2016).
[6] N. Dupertuis, et. al., Plasma Fusion Res. 12 1201008 (2017).
[7] N. Kasuya, et. al., Plasma Phys., 15, 052302 (2008).

2

Cross scale energy transfer in parallel flow shear driven
turbulence
Y. Kosuga[1,2][1] Research Institute for Applied Mechanics, Kyushu University, Fukuoka, Japan[2] Research Center for Plasma Turbulence, Kyushu University, Fukuoka, Japan

Flows along the magnetic field (parallel flows, in short) are ubiquitous in space andlaboratory plasmas. In particular, in fusion plasmas, parallel flows are observed in scrape
off layer plasmas[1] and ITB plasmas[2]. In these regions, parallel flows with the Mach
number close to unity are observed. These flows can be source of fluctuation, which impact
transport processes in other channels. Indeed, coupled dynamics of parallel momentum
transport and particle transport is observed on recent studies[3].
A key to understand the underlying dynamics is the interplay between drift wave
turbulence and parallel flow shear driven turbulence. While previous studies reveal abasic feature of parallel flow shear driven mode, there nonlinear physics remains elusive.
Once nonlinear coupling is invoked, cross-scale energy transfer comes into play. In thisprocess, energy in the primary mode can be distributed to different scales, including
larger scales. As a result, large scale structure, such as zonal flows[4,5] can be generated
in parallel flow shear driven turbulence.
In addition to zonal flows, several other secondary structures can arise as a consequenceof nonlinear interaction. In this work, we discuss generation of large scale convective cells
in parallel flow shear driven turbulence. Nonlinear energy transfer among primary parallel
flow shear driven mode and secondary convective cells is formulated. Performing closure
calculation, we obtain the nonlinear growth rate of convective cells. As a limiting cases,
we discuss the generation of i.) zonal flows, poloidally elongated cells and ii.) streamers,
radially elongated cells. It is shown that both structures can be realized in parallel flow
shear driven turbulence. Since two different structures are possible, a relevant question iswhich type of structure is formed. In other words, what is the selection rule of a structure
in parallel flow shear driven turbulence? From deterministic point of view, we discuss thatzonal flows are preferred. This is due to the fact that for cells with finite poloidal wave

number, sideband in parallel velocity fluctuation tends to prohibit the accumulation ofenergy in large scale. In a extreme case, we argue that while energy can be pumped into
zonal flows, radially elongated cells with finite poloidal wave numbers do not gain energy
and rather pump smaller scales. Relation to the observed high frequency component ofparallel flow velocity fluctuation is discussed.
Stimulating discussion with Drs. S.-I. Itoh, K. Itoh, S. Inagaki, T. Kobayashi, P.H.
Diamond, L. Wang, and Z. Guo is acknowledged. This work was supported by Grantsin-Aid for Scientific Research of JSPF of Japan (JP15H02155, JP15K17799), Kyushu
University Interdisciplinary Programs in Education and Projects in Research Development(TT26705) and QR program (28315).

References[1] X. Garbet, et al., Phys. Plasmas 6 3955 (1999)[2] X. Garbet, at al., Phys. Plasmas 9 3893 (2002)[3] S. Inagaki, et al., Sci. Rep. 6 22189 (2013)[4] P.H. Diamond, et al., Plasma Phys. Control. Fusion 47 R35 (2005)[5] Y. Kosuga, et al., submitted to Phys. Plasmas (2016)

直線装置における分布形成に関するドリフト波不安定性の非線形シミュレーション
Nonlinear simulation on profile formation with drift-wave instability
in linear devices
Y. Matsui 1, N. Kasuya 1,2, M. Sasaki 1,2, S. Abe 1, S. Inagaki 1,2 and M. Yagi 3
1
Interdisciplinary Graduate School of Engineering Science, Kyushu University, Kasuga,
Fukuoka 816-8580, Japan
2
Research Institute for Applied Mechanics, Kyushu University, Kasuga, Fukuoka 816-8580,
Japan
3
National Institute for Quantum and Radiological Science and Technology, Obuchi,
Rokkasho, Aomori 039-3212, Japan
E-mail: matsui.yousuke@riam.kyushu-u.ac.jp
Turbulent structural formation is one of the keys for understanding transport phenomena inmagnetized fusion plasmas. We have been carrying out the research of turbulent structuralformation with resistive drift-wave instability in linear devices [1]. Numerical Linear Device(NLD) [2] is a numerical simulation code of resistive drift-wave turbulence in linear devices,
and is used for quantitative predictions of experimental observations. In experiments in linear
device PANTA, bifurcation between turbulent structures, as streamer [3] and solitarydrift-wave [4], has been observed. In addition, a theory predicts the bifurcation of the positionof the steep wave front of the solitary drift-wave [5]. In this research, stability and structuralformation condition of resistive drift-wave turbulence are investigated by scaning of theplasma parameters using the NLD code to understand the parameter dependency and effect ofthe mean profile formation with turbulent structures.
A numerical simulation code has been developed for resistive drift wave in linear devices [2].
The plasma geometry is assumed to be a simple cylinder and the magnetic field has only theaxial component with uniform intensity. According to the experiments in PANTA device, we
focus on high-density and low-temperature plasmas in an argon discharge. Neutral particles
also exist inside the plasma and its effect must be taken into account. The code is based on thethree-field (density, electrostatic potential, and parallel velocity of electrons) reduced fluid
equations, which consist of the continuity equation, the charge conservation equation, and themomentum conservation equation. The model is an extension of the Hasegawa-Wakatanimodel by taking neutral particle effects into account. The periodic boundary conditions are
assumed in the azimuthal and axial directions, which allow the spectral expansion. We
describe the drift wave by nonlinear simulations using parameters; electron temperature Te = 2[eV], plasma radius a = 7 [cm], length of the device  = 4 [m]. Simulations are carried out
with a time independent particle source. The density profile in plasma is formed
self-consistently by balance of the particle supply, turbulence fluxes, and collisional transport.
In the edge region, plasma particles are lost by the collisional transport. Focusing on basic
nonlinear processes for turbulence caused by the density gradient, we do not consider the
moment source induced by a particle source here. The modes that has azimuthal mode numberm = 0 are linearly stable, but the nonlinear interaction of drift waves can drive these modes.
Using these parameters, the time evolution of each mode is calculated by solving the set of
equations.

Time evolutions of fluctuations are calculated with the fixed particle source term. Timeevolution of the energy of each Fourier component of the density fluctuation is shown in Fig.
(a). This is the case when the streamer is formed. The dominant modes, which are nonlinearly
coupled to each other, have (m, n) = (2, 1) and (3, 1), where m and n are azimuthal and axial
mode number, respectively. The density profile is flattened in the nonlinear evolution to makethe modes less unstable, and only these two modes are unstable in the saturated state. The
growth rate is evaluated by using the instantaneous background density and potential profiles.
In this case, the effect of the potential formation is little. The growth rate dependencies onmagnetic field B and collision frequency in of (m, n) = (2, 1) mode and (3, 1) mode arecalculated as in Fig. (b) and (c). The ion-neutral collision frequency in makes the drift-waves
stable. In the saturated state, it is important to distinguish which kind of state, solitary drift
wave or streamer, is formed, and co-existence of (2,1) and (3,1) modes or sole-existence of
one dominant mode is related to the selection. The tendency in Fig. (b) and (c) shows somedifference, and smaller B and larger in is possible to give the condition for sole-existence of
one dominant mode. In addition, it is necessary to consider the effect of the background
distribution formation, which is self-consistently taken place by nonlinear transport processes.
The ion-neutral collision frequency is also related to zonal flow damping, which gives aselection rule of turbulent structure formation [2]. By taking account both of the quasi-linear
and nonlinear properties of the instabilities, the structural bifurcation observed in theexperiment should be understood. In the solitary drift wave state, the steep gradient position is
selected in accordance with the radial structure of the fundamental mode [5]. The mode
structure is found to depend on the radial profile of the background density profile, so theselection will be investigated by nonlinear simulations with variation of the particle source
profile.

FIG. (a) Time evolutions of the energy of each Fourier component of the density fluctuation.
Dependency of the growth rate on magnetic field strength B and ion-neutral collision
frequency in of (m, n) = (b) (2, 1) mode and (c) (3, 1) mode.
Acknowledgements
This work is supported by JSPS KAKENHI Grant Number JP24760703 and JP16K06938, by
thecollaboration
programof
NIFS(NIFS15KNST082,
NIFS16KNXN323,
NIFS13KOCT001) and of RIAM of Kyushu University.
References[1] S. Oldenbürger, et al. Plasma Phys. Control. Fusion 54, 055002 (2012)[2] N. Kasuya, et al. Phys. Plasma 15, 052302 (2008)[3] T. Yamada, et al., Nature Phys. 4, 721 (2008)[4] H. Arakawa, et al., Plasma Phys. Control. Fusion 51, 085001 (2009)[5] M. Sasaki, et al. Phys. Plasmas 22, 032315 (2015)

Turbulence simulation of structural formation considering neutral particle profiles inlinear devices
Satoshi Abe1, Naohiro Kasuya 1,2, Makoto Sasaki 1,2, Yosuke Matsui 1,
Shigeru Inagaki 1,2, Masatoshi Yagi 3
1.Department of Advanced Energy Engineering Science, Interdisciplinary Graduate School of Engineering
Sciences, Kyushu University
2. Research Institute for Applied Mechanics, Kyushu University
3.National Institute for Quantum and Radiological Science and Technology
6-1 Kasuga-koen, Kasuga, Fukuoka 816-8580, Japan*E-mail: abe@riam.kyushu-u.ac.jp

ABSTRACT
Fundamentals of turbulent transport have been studied using linear devices. Neutral particles play an importantrole in turbulent structural formation in the devices, so their effects are investigated by numerical simulations oflinear device PANTA. Density profiles of neutral particles are calculated by a particle simulation using theMonte Carlo method. Then, the obtained profiles are introduced into a simulation code of resistive drift wave
turbulence. Comparison of cases with the ion - neutral collision frequencies constant or not-constant in radius
shows that non-uniformity of the neutral particle profiles can affect selection of fluctuation spectrum in the
turbulent state in the linear devices.
KEY WORDS
turbulence, simulation, linear device, collisional plasma, neutral particle, drift-wave instability, structuralformation.
1 INTRODUCTION
Turbulent transport is one of the keys for power generation by nuclear fusion, and basic experiments using
linear devices are being conducted for the understanding [1]. It is known that control of neutral particles isimportant for the turbulent state in the linear devices, and structural bifurcation has been shown in the simulationof resistive drift wave turbulence with the neutral particle as the control parameter [2]. There are neutral particle
distribution in the radial and axial direction, which must be taken into consideration. The purpose of this researchis to evaluate effects of the neutral particle profile on the plasma turbulent structure. Numerical simulations ofresistance drift wave turbulence are carried out with neutral particle profiles calculated by Monte Carlo method
in a cylindrical plasma for that purpose.
2 SETUP FOR NUMERICAL SIMULATION
Basic researches on plasma turbulence have been carried out in linear device PANTA at Kyushu University [3].
Figure 1 shows the schematic of the device configuration. Argon plasmas are generated by helicon wave with
device length ~ 4 [m] and plasma diameter ~ 10 [cm], plasma density at the center n0 ~ 1 × 1019 [m3] and theelectron temperature Te ~ 3 [eV]. Detailed analyses of plasma instability has been carried out. This plasma is thetarget for analysis in this research.
To calculate the neutral particle profile the particle code by the Monte Carlo method is used [4]. The Monte
Carlo method is a numerical technique obtaining an approximate solution by iterative calculation using random
numbers. The density and temperature profile of neutral particles are calculated, based on experimental
parameters in PANTA. Particles are supplied from the production tube and exhausted by the pump. Reemission
on the end plate and recombination in the plasma are considered. There are baffle plates installed in the linear
device. Collisional processes between neutral particles, neutral particles - ions (elastic; exchange of momentum
and charge), neutral particles - electrons (ionization) are taken into consideration. The weight of each neutral
particle is calculated, associated with flight of the particle. The elastic collisions between neutral particles
changes only their flight directions, and do not change energies and weights. The neutral - ion collisions affectthe particle direction, but do not change the weights. The weights are changed only with ionization processes bycollision with electrons. The test flight particles are traced until their weight become less than a defined
threshold. If the pressure of the neutral particles in the device is sufficiently high and the effective ionization rate
is low, the collisions between neutral particles plays an important role in determining the neutral particle
distribution, so this process is calculated by iteration to obtain the nonlinear convergence.
To calculate the drift wave turbulence fluid code NLD is used [2]. The set of three-field simplified fluid
equations including the effect of neutral particles is used to calculate the fluctuations of density N, electrostatic
potential φ, and electron flow velocity in the magnetic field direction V. Spectral expansion with Fourier modesFig.1 Schematic of linear device PANTA
are taken in the azimuthal and axial directions. By giving a particle source profile, nonlinear calculation iscarried out to determine the self-consistent density distribution, and the relationship between excited modes andsaturation states sustained by nonlinear couplings is studied.

3 RESULTS
3.1 Neutral Particle Profile
Firstly, two-dimensional neutral particle profiles are calculated using the Monte Carlo code. The profile with Te= 3 [eV] and axial flow velocity normalized by the thermal velocity M = 0.3 is shown in Fig.2. The radial neutral
particle profiles are obtained by averaging the two dimensional profile in the axial direction. The dependency on
Te and M is evaluated as in Fig.3. The neutral density in the radial direction is almost constant with the condition
Te = 3 [eV], but the higher electron temperature (Te =5 [eV]) gives decrease of the density near the center andstronger density gradient, because of increase of the ionization rate. In addition, the larger axial flow gives
decrease of the absolute value of the density. In this way, we obtained the difference of the neutral density profile
depending on the plasma parameters.
3.2 Drift Wave Turbulence
We have calculated the neutral density profile, and is introduced in the calculation of drift wave instabilities.
The profiles in the radial direction are almost constant in the case with Te = 3 [eV], and small at the center,peaked near the plasma edge and almost constant outside the plasma with Te = 5 [eV]. These profiles are fitted as
the red lines A and B in Fig.3. The difference of the saturation states of drift wave turbulence using the ion neutral collision frequency in profiles considering the neutral profiles are calculated. Figure 4 shows the time
evolutions of mode energies, where m and n are the azimuthal and axial mode numbers, respectively. In the casewith constant νin = 0.06, the modes with (m, n) = (2,1) and (3,1) become dominant, and, (2,1) + (3,1) mode
coupling forms a streamer [5] in the saturated state. With the profile B, the (2,1) mode is dominant, and theamplitude of the (3,1) mode becomes small, though this mode is linearly unstable. This selection is found to berelated to the formation of the mean potential profile. In this way, it is found that non-uniformity of the neutral
particle profiles can affect selection of dominant modes in the turbulent state in linear devices.
4 SUMMARY
The effects of neutral particle on structural formation of resistive drift wave turbulence in the linear device were
investigated. The density profiles of neutral particles were calculated by a particle simulation, and thedependencies on the electron temperature and axial flow velocity were examined. As the electron temperature
increases, the density profile becomes not flat in radius. As the result of introducing the neutral particle radial
profile into the turbulence code, the formation of the mean potential acts to stabilize one of the modes, and gives
selection of mode spectrum for saturation. In the future, quantitative analyses to find the selection rule on
turbulent structure formation should be made.

Fig.2 2-D neutral particle profile with Te = 5 [eV] and M = 0.3.

Fig3 Radial profiles of neutral particle density obtained by the Monte Carlo code. Fitted curves introduced into
the turbulence code are also shown.

Fig.4 Time evolutions of fluctuation energies with (a) flat A and (b) non-flat B profiles in Fig. 3.

ACKNOWLEDGMENTS
This work is partly supported by JSPS KAKENHI Grant Numbers JP24760703 and JP16K06938, by thecollaboration program of NIFS (NIFS15KNST082, NIFS16KNXN323, NIFS13KOCT001) and of RIAM of
Kyushu University.
REFERENCES
1 G. R. Tynan, et al., Plasma Phys. Control. Fusion 51 (2009) 113001.
2 N. Kasuya, et al., Phys. Plasmas 15 (2008) 052302.
3 S. Inagaki, et al., Sci. Reps. 6 (2016) 22189.
4 M. Ignatenko, et al., Jpn. J. Appl. Phys. 46 (2007) 1080.
5 T. Yamada, et al., Nature Phys. 4 (2008) 721.

Kinetic Simulation for Edge Plasma and Detached Divertor
T. Takizuka
Graduate School of Engineering, Osaka University
Kinetic effects play an important role in edge plasmas with steep gradient and including theopen-field region. Complete kinetic behaviors can be simulated only by the kinetic modeling.
One of the most powerful kinetic models is the particle simulation. The primary role of particle
simulation is to clarify the fundamental physics and to establish a numerical model of the
kinetic effect applicable to divertor codes with fluid modeling. Several results of PARASOL
simulation for the kinetic effects in edge plasmas are presented.

1. Introduction
The detached divertor plasma is considered a solution for the heat control to reduce thedivertor-plate erosion in fusion reactors with huge output power [1,2]. Synthetic reviews of thedivertor plasma detachment will be published in a special issue of Plasma Physics and Controlled
Fusion reviewing divertor plasma detachment in magnetic fusion devices [3].
There have been many experimental achievements on the detached divertor plasmas. To understandthe physical mechanisms of the detachment observed in experiments, a lot of effort has been put into
numerical studies using the comprehensive divertor simulation codes based on the plasma fluid
modeling, SONIC, SOLPS, UEDGE, EDGE2D, EMC3-EIRENE etc. The agreement in the results
between experiments and simulations has been rather fair but not satisfactory so far. Since these codes
are widely applied in the predictive simulations for future fusion reactors, it is indispensable toimprove the modeling in these codes.
Kinetic effect is one of key issues for improving the modeling for edge plasmas with steep gradientand including the open-field region. The mean free path is roughly given by lmfp0 [m] ~ 1016 × T2/n
[eV2/m-3], which is comparable to the connection length of SOL plasma in a fusion reactor withpossible parameters, LSOL ~ 100 m, T ~ 300 eV, n ~ 4×1019 m-3 and lmfp0 ~ 20 m. Even in the detacheddivertor plasma, the kinetic effect by energetic particles still plays an important role. Assuming nT
nearly constant along B, lmfp at an energy of E in the detached divertor region is estimated as lmfp ~(Tdiv/T) (E/T)2 lmfp0. Considering a detached-divertor plasma with Tdiv ~ 2 eV near the divertor plate
within the parallel connection length Ldetach ~ 2 m, lmfp becomes longer than Ldetach for energetic
particles with E/T > 4 (E > 1.2 keV). These energetic particles directly arrive at the divertor plate.
Plasma-wall interactions are influenced by the kinetic effects.
Several models of the kinetic effects are included into the plasma fluid equation, such as plasma-wallboundary conditions, heat transport coefficients parallel to B etc. Complete kinetic behaviors can be
simulated only by the kinetic modeling. One of the most powerful kinetic models is the particle
simulation [4-10]. The particle modeling for edge plasma simulations can easily treat the plasma-wallboundary conditions, and is very flexible in introducing A&M processes, compared with thecontinuum modeling of the Boltzmann equation. A review of the kinetic modeling for edge plasmaand detached divertor will be found in the above special issue of PPCF [11]. In the present
proceedings paper, several topics are picked up from this review article.
2. Particle simulation for whole SOL-divertor plasma
The particle model simulates the whole SOL and divertor plasmas including the sheath in front of
the divertor plate. Although the system size L is much larger than the Debye length λ! = ε! T! /e! n!
in real plasmas, usual particle simulations, such as PARASOL [9-11], with the grid size of order of λDare available to study such plasmas by setting artificially smaller values of L/λD = 102 ~ 104. This
artificial setting is allowable because the characteristics of SOL/divertor plasmas under thequasi-neutral condition, except in the sheath region, are determined mainly by collisionality L/lmfp andnormalized ion Larmor radius ρLi/a but insensitive to the L/λD value as shown in Fig. 1. For divertor
plasma simulations, the small L/λD value can be set by artificially elongating λD keeping the size L.
The permittivity ε0 in the Poisson equation, – ∇2 φ = e (Σi Zini – ne)/ε0, is artificially enlarged for thissake. As an artificially elongated λD is used, we adopt, at the same time, a “collision cut-off technique”
near the wall to keep the collisionless sheath condition valid.

Fig. 1 Potential profiles fordifferent L/λD values. The sameprofiles are established except
for sheath region with the sizeof ~3λD.

We estimate the computation time tcomp necessary for a particle simulation. The total number ofparticles Ntot in a 1D system is L/Δ times of the mean number of particles in a cell Ncell, and it becomes
Ntot = (L/Δ)2 × Ncell in a 2D system. The characteristic time τ for a SOL-divertor plasma equilibrium is
L/Cs (Cs: sound speed), and the required time steps is Kt ~ L/CsΔt ~ 𝑚! /𝑚! (𝐿/Δ), because vteΔt ~ Δ~ λD in the particle simulation. Then tcomp is estimated as tcomp = Ntot × Kt × FLOΔt / FLOPS (FLOΔt:number of floating-point operations for a single time-step advance of a particle, FLOPS: effective
operating speed of a computer system). For example a small 1D case (L/Δ = 103, Ncell = 102, Kt = 105,
and FLOΔt = 103), it takes about 3 hours using a 1 GFLOPS computer. For a 2D case with parameters
similar to the above, tcomp is still within several hours by the use of a 1 TFLOPS parallel computer.
When we try a large and complex 2D system (L/Δ = 104, Ncell = 103, Kt = 107, and FLOΔt = 104), tcomp
becomes up to 100 days even if using a 1 PFLOPS massively parallel computer.
3. PARASOL simulation results of edge plasmas
3.1 Heat conduction parallel to the magnetic field
In the plasma fluid model, the heat conduction parallel to B is given by a collisional diffusive one,
the Spitzer-Harm expression, q// = qSH ≡ – κ// ∇//T. When lmfp > L//, q// is no more the diffusive transport
but becomes the free-stream-like transport; q// = α qFS ≡ α nT vt. To bridge the gap between short lmfp
case and long lmfp case, a harmonic average has been applied; 1/q// = 1/qSH + 1/αqFS. Here α is theflux-limiting coefficient, but this value is not conclusively obtained.
The α value is directly measured with the PARASOL simulation [12]. Hot plasma is supplied from the
hot source region in the central SOL, and is lost to the divertor plates. Electron energy is cooled byradiation in the radiation cooling region near the plate (Fig. 2 (a)). For the various cases of collisionality L///lmfp (10-3 ~ 103) and radiation fraction frad, qe// for electron is measured in the intermediate
region (Fig. 2 (b)). The PARASOL simulation results agreed with the classical theory, qe// = qSH, in thecollisional regime, L///lmfp >> 1. On the other hand in the collisionless regime, L///lmfp << 1, qe// becomesmuch smaller than qSH. When comparing qe// with qFS, the ratio qe///qFS decreases in inverse proportion
to L///lmfp in the collisional regime, while it saturates at a certain value in the collisionless regime. Thissaturated value was measured as αe. It was found that the variety of αe value is dependent on frad; αe<< 1 for small frad (Tediv ~ TeSOL) and αe ~ 1 for large frad (Tediv << TeSOL) as shown in Fig. 2 (c).

Fig. 2 Poloidal profiles of(a) Te and (b) qe// for various
cases of L///lmfp and frad. Hot
electrons are supplied in the
hot source (H.S.) region andcooled in the radiation cooling(R.C.) region. Flux limiting
factor αe measured in theintermediate region is shownin (c) as a function of frad.

3.2 SOL flow
The SOL flow pattern is varied by the null-point location in tokamak divertor plasmas. Many
simulation studies on the SOL flow have been carried out using fluid codes, but the agreement withexperiments were not satisfactory so far. PARASOL 2D simulation was carried out to clarify thephysics of SOL flow pattern variation [13].
Figure 3 shows the 2D structures of V// in a tokamak with an aspect ratio A = 5.5 and the downward!"#𝑉∇! . The red region corresponds to the co-flow to the plasma current and the blue region corresponds
to the counter-flow. The flow pattern is changed by the location of a null point (NP); the upper-null!"#(UN) or the lower-null (LN) configurations. For the UN case with 𝑉∇!
away from NP (Fig. 3 (a)), V//
is directed to the divertor plate both in the inner (high-field side, HFS) and outer (low-field side, LFS)
SOL regions and the stagnation point (V// = 0) is located symmetrically at the bottom. On the other
hand for the LN case (Fig. 3 (b)), V// in the outer SOL region has a backward flow pattern (red region).
The stagnation point moves below the mid-plane of the outer SOL.
Based on these PARASOL simulation results, a new model for the edge plasma flow by introducing!"#
the “ion-orbit-induced flow” was presented [14]. Due to 𝑉∇!
in a tokamak, an ion with the co-v// isshifted outward (LFS), while an ion with the counter-v// is shifted inward (HFS). A banana orbit of atrapped ion drifted out from the core generates the co-flow (red region) in the outer SOL region. In a
UN configuration, this “ion-orbit-induced flow” is additive to the original SOL flow towards the upper
divertor plate in the outer SOL. On the contrary, in a LN configuration, the “ion-orbit-induced flow” issubtractive from the original SOL flow towards the lower divertor plate in the outer SOL. Transit ionsas well play a role in the “ion-orbit-induced flow”. This “ion-orbit-induced flow” causes the remarkable difference in flow patterns for the UN- and LN-divertor configurations.

Fig. 3 Structures of V// in atokamak with A = 5.5. For UN
case (a), the “ion-orbit-induced
flow” is additive to the original
SOL flow towards the upper
divertor plate in the outer SOL.
For LN case (b), the “ion-orbitinduced flow” is subtractive fromthe original flow towards thelower plate in the outer SOL.

5. Summary and discussion
The particle simulation solving the full ion orbit is the closest to first principles, and has a highability to understand the physics of edge plasmas and detached divertor. Although the kinetic effectbecomes less important in the very collisional detached divertor region, the kinetic effect by energetic
particles still plays an important role. In addition, the down-stream detached divertor plasma is
influenced by the upstream SOL plasma, in which the kinetic effect is one of the essential factors.
Due to tremendous computational resources it seems too challenging for the particle simulation towork directly on the divertor design of future reactors. The primary role of the particle simulation is toclarify the fundamental physics and to establish a numerical model of the kinetic effect applicable todivertor codes with fluid modeling. Benchmark test between particle and fluid simulations before andafter installing the kinetic-effect model is important. The second role is to study the whole tokamak ofsmaller size or in short time scale.
Time-separated integration of the codes can be a future project for the reliable design oftokamak-reactor divertor; 3D particle code for the turbulent diffusion in the shortest time scale (say ~1ms), 2D particle code for the parallel transport in a shorter time scale (10~100 ms), 2D comprehensive
divertor code to respond the sec order edge phenomena, and 1.5D transport code to follow fully theoperation period.
Acknowledgements The author thanks many collaborators on the simulation and modeling for edge
plasmas in fusion devices: Dr. K. Shimizu (QST), Dr. K. Hoshino (QST), Dr. N. Hayashi(QST), Dr. H.
Kawashima (QST), Dr. M. Yagi (QST), Mr. M. Hosokawa (ITER), Dr. A. Loarte (ITER), Dr. A. Froese, Mr.
S. Azuma, Dr. K. Ibano (Osaka U.), Dr. S. Togo (U. Tsukuba), Prof. A. Fukuyama (Kyorto U.), Prof. Y.
Ueda (Osaka U.), and Prof. Y. Ogawa (U. Tokyo). Majority of simulation studies using PARASOL havebeen carried out at Japan Atomic Energy Agency, and some studies at Kyushu University, Kyoto
University, and ITER organization.

References[1] ITER Physics Basis, Chap. 4: Power and particle control, Nucl. Fusion 39 (1999) 2391.
[2] Progress in the ITER Physics Basis, Chap. 4: Power and particle control, Nucl. Fusion 47 (2007) S203.
[3] Reviewing divertor plasma detachment in magnetic fusion devices, edited by R. Maingi, T. Rognlien,
to be published in Plasma Phys. Control. Fusion.
[4] R.W. Hockney, J.W. Eastwood, Computer Simulation Using Particles (New York, McGraw-Hill, 1981).
[5] C.K. Birdsall, IEEE Trans. Plasma Sci. 19 (1991) 65.
[6] R. Cohen, X.Q. Xu, Contrib. Plasma Phys. 48 (2008) 212.
[7] J.P. Verboncoeur, Plasma Phys. Control. Fusion 47 (2005) A231.
[8] D. Tskhakaya, K. Matyash, R. Schneider, F. Taccogna, Contrib. Plasma Phys. 47 (2007) 563.
[9] T. Takizuka, AIP Conf. Proc. 1237 (2010) 138.
[10] T. Takizuka, Plasma Sci. Technol. 13 (2011) 316.
[11] T. Takizuka, “Kinetic effects in edge plasma: kinetic modeling for edge plasma and detached divertor”,
to be published in Plasma Phys. Control. Fusion 59 (2017).
[12] A. Froese, T. Takizuka, M. Yagi, Plasma Fusion Res. 5 (2010) 026.
[13] T. Takizuka, K. Shimizu, N. Hayashi, M. Hosokawa, M. Yagi, Nucl. Fusion 49 (2009) 075038.
[14] T. Takizuka, K. Hoshino, K. Shimizu, M. Yagi, Contrib. Plasma Phys. 50 (2010) 267.

Research on behavior of divertor particle flux with a
one-dimensional numerical model and an analytical one
Satoshi Togo, Tomonori Takizukaa) , Makoto Nakamurab) , Kazuo Hoshinoc) , Kenzo Ibanoa) , Yue Lid) ,Yuichi Ogawad) , Mizuki Sakamoto
Plasma Research Center, University of Tsukuba, 1-1-1 Tennodai, Tsukuba 305-8577, Japan
Graduate school of Engineering, Osaka University, 2-1 Yamadaoka, Suita 565-0871, Japanb)
National Institutes for Quantum and Radiological Science and Technology, 2-166 Omotedate, Obuchi-aza, O-aza,Rokkasho 039-3212, Japanc)
National Institutes for Quantum and Radiological Science and Technology, 801-1 Mukoyama, Naka 311-0193, Japand)
Graduate School of Frontier Sciences, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa 277-8568, Japana)

Behavior of particle flux at the divertor plate Γt as a function of that from the core plasma Γsep has beenstudied by using a one-dimensional fluid model and a reduced analytical one focusing on the eﬀects of ion
temperature anisotropy, radial width of the scrape-oﬀ-layer (SOL) and radial plasma diﬀusion. Bifurcation of Γt
has been observed numerically and analytically. The eﬀect of ion temperature anisotropy on Γt was small in thecurrent calculation condition. In wider SOL case, the change of Γt became slower against Γsep . The region ofmultiple solutions became narrower in larger plasma radial diﬀusivity case.

1 Introduction
In order to design divertors which will escape serious
damages due to the plasma particle flux, it is necessary
to develop a numerical simulation code packages which
can reproduce the detached divertor plasmas observed inexperiments of existing torus devices. Numerical resultson the behavior of the divertor particle flux from theconventional code packages, however, do not always agreewith the experimental ones. Therefore, various eﬀects,
such as kinetic eﬀects, supersonic flow, improved atomic
and molecular processes and plasma wall interactions,have been studied in order to improve the numerical resultsof the conventional code packages [1–3]. Compared tothese code packages, one-dimensional (1D) models are
more convenient to gain physical insights because of its
simplicity and short calculation time. In this paper, westudied the behavior of particle flux at the divertor plateΓt as a function of that from the core plasma Γsep byusing a 1D fluid model [4–7] focusing on the eﬀects of ion
temperature anisotropy, radial width of the scrape-oﬀ-layer(SOL) and radial plasma diﬀusion in divertor (DIV) region.
We also used a reduced analytical model in order tounderstand the structure of the numerical solutions.

2 Models
2.1 Numerical model
We have been developing a 1D fluid model of the
SOL-DIV plasma and neutrals which distinguish between
the parallel, T i,∥ , and perpendicular, T i,⊥ , components ofauthor’s e-mail: togo@prc.tsukuba.ac.jp

ion temperature [4–7]. Because of the space limitation, we
omit the explanations on the details of the fluid equations
for the plasma and neutrals in this paper. The fluid
equations and notations are the same as those used in Ref.
[7]. In the isotropic T i model, the following equationsof parallel momentum transport and isotropic ion energy
transport were solved in stead of the anisotropic ones;
(
)
∂
∂V
∂
2
(mi nV) +mi nV + nT i + nT e − ηi
∂t
∂s
∂s
(1)
= Mm ,()∂ 1
3
2
mi nV + nT i
∂t 2
2
(
)
∂ 1
5
∂T i
∂V
3
(2)
+
mi nV + nT i V − κi
− ηi V∂s 22∂s∂s3me n (T e − T i )∂= Qi +− V (nT e ) .
miτe∂sHere, ηi is the limited Braginskii ion parallel viscosity;
ηi =

ηBr
i

ηBr | ∂V |1 + iβnT∂si, ηBr
i = 0.96nT i τi ,(3)

where the limiting factor was set to be β = 0.7 according toRef. [5]. The ion heat source term Qi is given as follows;()
3(νiz + νcx )Qi = Qcore
+
ϵ
n
+
T
n
FC
n,recy
i
n,diﬀ
i
2
)
(
3
1
(4)
− mi V 2 + T i (nn νcx + nνrc )
2
2
(
)
1
3
− mi V 2 + g⊥i T i ν⊥ ,22where the ion heat source from core plasma Qcore
andi2
radial diﬀusion frequency ν⊥ = D⊥ /dDIV
have non-zerovalues in SOL and DIV regions, respectively. The
coeﬃcient g⊥i , which represents the eﬀect of radial heatdiﬀusion, is set to be 3.

2.2 Reduced analytical model
In order to understand the structure of the numericalsolutions and obtain steady state solutions including
unstable ones, we have been developing an analyticalmodel. This model is motivated by the fact that neutrals
exist in a very thin region in front of the plate compared tothe system length in the numerical results. In this model,
reduced neutral equations are solved in homogeneousbackground plasma of density n, flow velocity V andsingle temperature T described by the sheath boundary
conditions. For a given value of T and tentative total target
heat flux Pt [W], total ion flux is given by Γt = Pt /γT [/s]where γ is the sheath heat transmission factor. Then, the√flow velocity and density are estimated as V = 2T/miand n = Γt /2V A∥ where A∥ = (dSOL /LSOL ) Asep sin θ is the
cross section of the flux tube normal to the magnetic field.
In the homogeneous plasma described above, thefollowing reduced neutral equations are solved;
inn
nout
ninn
n,recy dnn,recy
n,recy
,
=
,
λr
dx
λr

(5)

d2 nn,diﬀ
nn,diﬀ= −Dn 2 + nn,recy νcx .
2
dx
λd

(6)

dnout
n,recydx−Dn=−Here, the volume recombination is neglected.
Thecharacteristic lengths, which are both constants for a given
background plasma, are given as follows;√
Vn,recyDnλr =, λd =.
(7)
νcx + νiz + νL,recy
νiz + νL,diﬀ
The solutions to Eqs. (5) and (6) becomenn,recy = nn,recy,t
[
(
(
)
)]
x−L
x+L
× exp
+ exp −,λr
λrnn,diﬀ = f nn,recy − f gnn,recy,t
[
(
)
(
)]
x−L
x+L
× exp
+ exp −,λd
λd
V sin θn,Vn,recy + (g − 1)Cn f vn

−1
 1
Dn /λr + Cn vn
1  νcx
f =  2 − 2 
, g=
,
Dn
Dn /λd + Cn vn
λd λr

nn,recy,t =

(8)

(9)

Pitch angle θ
Parallel system length 2L/ sin θ
Parallel SOL length LSOL / sin θ
Parallel DIV length LDIV / sin θ
SOL width dSOL
DIV width dDIV
Separatrix surface area Asep
Radial particle diﬀusivity D⊥
Carbon impurity fraction rimp
Heat flux from core Psep
Ion flux from core Γsep
Ion temperature T i

6◦
44 m
35.2 m
4.4 m
1 (2) cm
5 cm
40 m2
0.5 (2) m2 /s
0.01
4 MW
0.5-4.5×1022 /saniso. (iso.)Nn,diﬀ = 2 f nn,recy,t (λr − gλd ) A∥ / sin θ.

(13)

Particle fluxes from ionization, charge exchange,
radial diﬀusion and core become as follows, respectively;
Γiz,r = Nn,recy νiz , Γiz,d = Nn,diﬀ νiz ,(14)

Γcx,r = Nn,recy νcx , Γcx,d = Nn,diﬀ νcx ,(15)Γ⊥ = −2nν⊥ LDIV A∥ / sin θ,(16)

Γsep = Γt − Γiz,r − Γiz,d − Γ⊥ ,(17)

where LDIV is poloidal length of DIV region. Heat
fluxes from ionization, charge exchange, radial diﬀusion,
impurity radiation and core become as follows;()
3
Piz = (ϵFC − ϵiz ) Γiz,r + T − ϵiz Γiz,d ,
(18)
2
)]
[
(
3
1
2
Pcx = ϵFC − mi V + T Γcx,r
2
2
(19)
1
2
− mi V Γiz,d ,
2
(
)
1
2
⊥
P⊥ = − mi V + 3g T Γ⊥ ,(20)
2
Prad = −2Lz (T ) rimp n2 λd A∥ / sin θ,(21)

Psep = Pt − Piz − Pcx − P⊥ − Prad .

(22)

Then, the value of Pt is iteratively controlled such thatobtained Psep has the prescribed value.
(10)

3 Results(11)

where L is the poloidal length from the stagnation point to
the divertor plate. By integrating these neutral densities
over x, the total number of the neutral particles are
determined as follows;
Nn,recy = 2nn,recy,t λr A∥ / sin θ,

Table 1 Basic calculation conditions. The parentheses represent
the compared ones.

(12)

We investigated the eﬀects of ion temperature anisotropy,
SOL width dSOL and radial particle diﬀusivity D⊥ on the
behavior of Γt as a function of Γsep . Table 1 shows
the basic and compared calculation conditions we used.
The results are shown in Fig. 1. For all calculation
conditions, we obtained lower (2 × 1017 ∼ 3 × 1019
/m3 ) and higher (8 × 1019 ∼ 5 × 1020 /m3 ) densityRpΓt (1022 /s)

25
20
15
10
5
0
1.0
0.8
0.6
0.4
0.2
0.0

because Eq. (21) remarkably underestimated the impurity
radiation power. Improvement of the impurity radiation
model in reduced analytical model is one of our future
works. In the numerical upper solutions, we obtained therollover of Γt which represents the onset of the detachedplasma.
In the current calculation conditions, we did notobserve a remarkable diﬀerence of the results fromanisotropic and isotropic T i models. This is due to theeﬀects of ion parallel viscosity limitation for the lower
solutions and the high collision frequency for the upper
ones. With wider SOL width, the changes of Γt and otherparameters became slower against Γsep because A∥ became
simultaneously large and required larger Γsep in order tomake the same plasmas. In larger D⊥ case, the solutions
shifted like the wider SOL width case only for the upper
ones because the radial diﬀusion was not large due to low
density in the lower ones. Note that the region of multiple
solutions became narrower in this case, which indicatesmuch larger D⊥ may extinguish that region.

(a)

(b)

nt (/m3)

1021
1020

(c)

1019
1018
1017
106

Vt (m/s)

(d)
105
104

Tt (eV)

103
103
102

(e)

4

101
100
10-1

0

1

2

3

4

5

Γsep (1022 /s)

Fig. 1 (a) Particle flux Γt , (b) momentum flux amplification R p ,(c) density nt , (d) flow velocity Vt and (e) (electron)
temperature T t at the sheath entrance as functions ofparticle flux from core plasma Γsep . Circular, square,
triangular and rhombic symbols represent the numericalsolutions with basic calculation condition, isotropic T i ,wider SOL width (dSOL = 2 cm) and larger radialdiﬀusivity D⊥ = 2 m2 /s, respectively, where closed andopen ones correspond to lower and upper stable solutions.
Solid, broken and dotted lines stand for analyticalsolutions with basic calculation condition, wider SOL
width and larger radial diﬀusivity, respectively.

branches of stable numerical solutions. Here, we callthe former the lower solutions (closed symbols) and thelatter the upper ones (open symbols) according to theway used in Ref. [8]. These branches share some rangeof Γsep (e.g. 1.1 ∼ 1.9 × 1022 /s for basic calculation
conditions) and the actual branch of a solution depends
on the branch the initial condition belongs to in thisrange. By using the reduced analytical model, we obtained
S-shaped solutions. This kind of bifurcated behavior of thesolutions agrees with those pointed out by earlier studies
including the intermediate branch in the analytical solution
which is not found by the numerical model because it isphysically unstable [8, 9]. Although the analytical solution
well reproduced the numerical lower solutions, it did not
reproduce the upper ones even qualitatively. This wasConclusion

We investigated the behavior of Γt as a function of Γsepby using a 1D fluid model and a reduced analytical one
focusing on the eﬀects of anisotropic T i , radial widthof the SOL and radial plasma diﬀusion. Bifurcation
and rollover were observed in numerical solutions whichwere reproduced by the analytical model as for the lower
ones. The eﬀect of T i anisotropy was small in thecurrent calculation condition. In wider SOL case, thechange of Γt became slower against Γsep . The region ofmultiple solutions became narrower in larger plasma radialdiﬀusivity case.
Improvement of the impurity radiation model in
reduced analytical model is one of our future works. In
addition, we will study the eﬀect of anisotropic T i withmore realistic conditions such as inhomogeneous magnetic
field and finite plasma current.

References[1] A.V. Chankin et al., J. Nucl. Mater. 390-391 (2009)319.
[2] K. Hoshino et al., J. Nucl. Mater. 415 (2011) S549.
[3] M. Groth et al., J. Nucl. Mater. 415 (2011) S530.
[4] S. Togo, et al., J. Nucl. Mater. 463 (2015) 502.
[5] S. Togo, et al., J. Comput. Phys. 310 (2016) 109.
[6] S. Togo, et al., Contrib. Plasma Phys. 56 (2016) 729.
[7] S. Togo, et al., Proc. 13th Burning Plasma Simulation
Initiative (BPSI) Meeting (Kasuga, Japan 2015) 3-3.
[8] M. Sugihara, et al., J. Nucl. Mater. 128-129 (1984)114.
[9] H.X. Vu and A.K. Prinja, Phys. Fluids 5 (1993) 1809.

Vacancy cluster growth in tungsten induced by hydrogen
Kazuhito Ohsawa1, Hideo Watanabe, 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 have excellent properties as plasma facing
materials (PFMs), e.g., low hydrogen (H) solubility and high melting point etc. However,a large amount of H isotope (D and T) retention has been reported in tungsten materials.
In particular, tritium inventory is an important subject associated with safety of fusion
reactors. H isotopes are not retained in W crystal lattice because of its low H solubility.
However, H isotopes can be accommodated in vacancy and vacancy-type lattice defects.
Stable configurations of multiple H atoms trapped in a W mono-vacancy have beenreported in the previous works [1].
In the present paper, we pay attention to stable configurations of di-vacancy in W
lattice and H atoms trapped in the di-vacancy. Beside, stability of the di-vacancy andvacancy-hydrogen (V-H) complex are investigated. Formation energy of di-vacancy isestimated to be negative. So, it is expected that di-vacancy nucleation does not takeplace [2]. On the other hand, V-H complexes combine because the binding energy isestimated to be positive. Therefore, the presence of H or H isotopes will enhance thegrowth of vacancy-type lattice defects.
2. Methodology
Nucleation of di-vacancy and binding of V-H complexes are investigated in terms offirst-principle calculations on the basis of density functional theory. The calculations are
performed by using Vienna ab initio simulation package (VASP) with PBE typepotential. The cut-off energy of plane wave is 350eV. We used large simulation cellcomposed of 432 lattice points (6x6x6 BCC lattice) in the present calculations in order toreduce the effects of periodic boundary condition imposed on the simulation cell. Lattice
relaxations are iterated until resultant force of each atom is lowered than 0.003 eV/Å.
Binding energy is estimated, as follows.

E B  2 EV   EV2 

EB  EVHm   EVHn   EV2 Hmn  ,

(1)
(2)

where V and H are represented vacancy and hydrogen, respectively. Integer m and nindicates the number of hydrogen atoms trapped in each vacancy and function E is thecohesive energy of the system including such lattice defects. Equations (1) and (2)indicate binding energies of two empty vacancies and V-H complexes, respectively.
Positive binding energy indicates attractive interaction.
3. Results
We calculate binding energies between two vacancies at first and second nearest
neighbors in a variety of BCC transition metals, that is, the two vacancies line up in the<111> or <100> direction. Fig. 1 shows the binding energies corresponding to stabledi-vacancy structure. Di-vacancy of <111> arrangement is more stable for in Cr, Mo, andW, while di-vacancy of <100> arrangement is favorable in Fe. Fig. 1 also indicatesbinding energies of V-H complexes in the BCC metals. Stable configuration of
di-vacancy hydrogen complex ( V2 H12 ) in Fe is <111> arrangement, that is, thedi-vacancy structure changes as the number of H atoms increases in the di-vacancy.
Fig. 2 shows binding energies of V-H complexes in tungsten in detail. We investigate
binding energies between V-H complex ( V1H n ) for n=1, 2, 3, and 6, where n is thenumber of H atoms trapped in a W vacancy. Di-vacancy hydrogen complex ( V2 H 2 n ) is
more stable than isolated V-H complex ( V1H n ). Only, the binding energy between two
empty vacancies is estimated to be negative. Therefore, it is expected that the presenceof H atoms enhances binding of isolated vacancies in W specimen.

Fig. 1: Binding energies of vacancies

Fig. 2: Binding energies of a variety ofand V-H complexes in BCC metals.

V-H complexes in tungsten.

4. Discussion
In usual metals, vacancy cluster growth is observed by absorbing isolated vacancies
in an appropriate high temperature. If formation energies of such vacancy-type latticedefects are assumed to be proportional to the surface area of the defects, these
phenomena can be explained adequately. Vacancy clusters grow so as to reduce theformation energy. This assumption is supposed to be valid for many kinds of BCC
transition metals, as shown in Fig. 1. However, di-vacancy formation in W will not takeplace exceptionally, according to the present works. On the other hand, V-H complexesin W can combine because the binding energies are positive. Fig. 3 is the schematic viewof the effect of H on the vacancy binding in W lattice, that is, H atoms trapped in W
vacancy enhances di-vacancy formation. Positron annihilation measurement is asensitive experimental method to investigate vacancy-type lattice defects in metals.
Usually, di-vacancy formation in metals is proved by drastic extension of the positron
life-time at vacancy diffusion temperature. We hope that the behaviors of vacancies in
W specimen will be revealed by such experiment.
In the previous paper, we pointed out an abnormal property associate with H atom
structure in W vacancy. For example, symmetric octahedral structure is stable for six H
atoms in usual BCC transition metals. However, an abnormal structure of H atoms waspredicted in W vacancy, as shown in Fig. 4 [1]. Similarly, abnormal properties are
expected associated with di-vacancy formation [2] and enhancement of vacancy cluster
growth by H presence in W specimen. We do not conceive an appropriate reason why Wmaterials have such abnormal properties but it will be an important subject of PFMs inthe future.

Fig. 3: Schematic view of binding of V and

Fig. 4: (a) Octahedral structure of

V-H complex in W specimen. It is predicted

six hydrogen atoms in a vacancy.

that Di-vacancy formation does not occur.

(b) Abnormal structure of H atoms.

References[1] K. ohsawa et al., Phys. Rev. B 85 (2010) 184117.
[2] L. Ventelon et al., J. Nucl. Mat. 458 (2015) 187.

Simulation study on internal transport barrier formation using gyrofluid model
M. Yagi, Y. W. Cho1, H. Seto and T. S. Hahm1
National Institutes for Quantum and Radiological Science and Technology,
Rokkasho Fusion Institute, 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. In QST, the project for
multi-scale simulation with transport/MHD is going on. On MS_Turb Project (IFERC-CSC),
the quasi-linear modification of neoclassical flows is taken into account for GF-ITG code andwe revisit the internal transport barrier (ITB) formation investigated by Tokunaga et al[1]. We
report the simulation results on ITB formation by the tuned heat source.
2. Model equations
GF-ITG code consists of three field equations (ion continuity equation, parallel momentum
equation and ion temperature evolution equation)[1]
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[2,3]µ1 =0.66ε 1/2ν i(1+1.03ν *i1/2 + 0.31ν *i )(1+ 0.66ε 3/2ν *i )(6)κ nc =%
1 " 1.17 − 0.35ν *i1/2− 2.1ν *i2 ε 3 '2 3$
1/2
1+ ν *iε # 1+ 0.7ν *i&(7)For the normalization, Bohm time t B = a 2 χ B , χ B = cTe eB and the minor radius a areused [1].
3. Quasi-linear effect on transport
In the previous study [4], we have neglected thequasi-linear (QL) term in the convective nonlinearity
of density assuming the effect is small, namely
Φ, 𝑛 = 0, where 𝑛 = Φ − Φ . However, for 𝑚 ≠ 0
and
-./,/
-0

𝑛≠0

,

the

term

[Φ, 𝑛]+,, =𝑖𝑘3 Φ+,, appears in Eqs.(1) and (3). Figure 1

shows the time evolution of ion temperature profile.
The dependence of quasi-linear effect on transport at
𝑡 = 5 is shown. In the simulation, the Gaussian type ofthe heat source 𝑆 = 𝑃𝑒𝑥𝑝(− 𝑟 − 𝑆0 > Δ> ) with 𝑃 =

Fig.1

Time

evolutionofion

temperature profile.

2.0, 𝑆0 = 0.4, Δ = 0.05 is applied at 𝑡 = 2 (shown by the dashed curve) then, the time
evolution of ion temperature profile is examined. Note that the neoclassical ballooning mode
(NBM) is turned off for this case. It is found that the QL effect is small contribution around
ITB region for the case without NBM.
4. Dependence of heat source location on ITB formation
Next, we examine the dependence of heat source location
on ITB formation. In this simulation, the QL effect of theconvective nonlinearity of density is NOT taken into
account. The source location is controlled by the parameter𝑆0 . Figure 2 shows the heat source profile. 5 cases with
𝑆0 = 0,03,0.4,0.5,0.6 are investigated.

Figure 3 shows

the ion temperature profile at 𝑡 = 12 for each case. In thiscase, the heat source is applied at 𝑡 = 5 (shown by thedashed curve). It is found that the weak ITB forms for
𝑆0 = 0.4. Note that the q minimum position is locatedFig.2 Heat source profile.

φ!(r,m,z=0)mr

Fig.3 Time evolution of ion temperature

Fig.4 Turbulence intensity on (r,m)

profile for each location of heat source.

space.

at 𝑟 = 0.6 (shown by the blue curve). Figure 4 shows the turbulence intensity (in thelogarithmic scale) at 𝑡 = 20 on (r,m) space. It is seen that inside of ITB, the turbulence issuppressed, on the other hand, the low m modes develop outside of ITB which produce thelarge convection.
5. Discussion
We perform the flux driven simulation taking account of the QL modification of neoclassicalflows and investigate the time evolution of ion temperature profile. At first, the effect of QL
term in the convective nonlinearity of the density is examined. It is shown that this effect isnegligible around ITB region if NBM is not taken into account. Then, we investigate thedependence of heat source location on ITB formation. A new finding is that the weak ITB is
formed when the Gaussian type heat source is applied at 𝑟 = 0.4.
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] S. Tokunaga et al., Nucl. Fusion 49 075023(2009).
[2] 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.
[3] F. L. Hinton and R. D. Hazeltine, Rev. Mod. Phys. 48 239(1976); κ nc is given byeq.(6.136).
[4] M. Yagi, H. Seto, Y. W. Cho and T. S. Hahm, the procs of 13th Burning Plasma
Simulation Initiative (BPSI) Meeting, Dec. 10-11, 2015, Kyushu Univ., Kasuga,
Japan. http://www.riam.kyushu-u.ac.jp/sosei/bpsi/bpsi13.html

14th Burning Plasma Simulation Initiative (BPSI) Meeting

H. Seto and M. Yagi

Linear analysis of non-ideal ballooning mode instability
with real electron inertia
H. Seto and M. Yagi
National Institutes for Quantum and Radiological Science and Technology,
2-166 Obuchi-Omotedate, Rokkasho-mura, Aomori, 039-3212, Japan
Shaping eﬀects on non-ideal ballooning mode instabilities are numerically investigatedby the use of two diﬀerent geometries, a shifted circular and a D-shape equilibrium,with the same pressure shear and magnetic shear profile based on a linearized GradShafranov equilibrium. The result is consistent with a local analysis indicating thatthe ellipticity κ > 1 stabilizes ballooning-type instabilities (BIs) by weakening thepoloidal magnetic curvature and the positive triangularity δ > 0 may also stabilize
BIs by deepening the magnetic well on the bad curvature plane.

1

Introduction

The H-mode discharge operation [1] is indispensable for a DEMO reactor from the view point of theconfinement eﬃciency but the intermittent large heat flux driven by edge localized modes should be avoided
for the engineering constraint of the heat load on plasma facing components [2]. It is therefore one of theimportant issues to understand shaping eﬀects on the edge plasma stability, since the ellipticity κ > 1 andthe positive triangularity δ > 0 improve the plasma confinement property. According to a local analysisof shaping eﬀects on the current diﬀusive ballooning mode instability [3], the ellipticity stabilizes theballooning-type instabilities (BIs) explicitly by weakening the poloidal magnetic curvature and the positive
triangularity may also stabilize BIs implicitly by deepening the magnetic well on the bad curvature plane.
In this paper, impacts of shaping eﬀects on the non-ideal ballooning modes (NIBMs) are numericallyinvestigated by BOUT++ code [4] with two diﬀerent geometries that have the same pressure shear andthe magnetic shear profile. One has a shifted circular geometry and the other has a D-shaped geometry
and the both are based on a linearized Grad-Shafranov equilibrium (LGSE) [5]. For describing NIBMs,
we employ a 3-field reduced MHD model consisting the resistive ballooning mode (RBM) [6,7] and thecollisionless ballooning mode (CBM) [8,9], where the other non-ideal eﬀects are neglected [10].
The rest of this paper is organized as follows. In section 2, brief descriptions of the LGSE and the
3-field RBM/CBM model are given. Section 3 gives simulation results and a summary.

2Analytic equilibrium model and 3-field RBM/CBM model

The Grad-Shafranov equation Eq.(2.1) is a nonlinear partial diﬀerential equation (PDE) for the poloidal
flux function ψ requiring two flux functions such as the pressure gradient dp/dψ, the flux function IdI/dψ,
(
)
1
dp
dI
R2 ∇ ·∇ψ= µ0 R2
−I
,
(2.1)
2
R
dψ
dψ
where R is the major radius and µ0 is the permeability in vacuum. For comparison with the analytic
analysis, an analytic solution of Eq.(2.1) with the fixed p(r) and q(r) profile is desirable, where r is theminor radius. To linearize Eq.(2.1), we define Fourier series for all quantities and expand them with respect
to the inverse aspect ratio ε = r/R. With these approximations, the horizontal position R and the vertical
position Z can be expressed with respect to non-orthogonal torus (NOT) coordinates (r, ϑ, ζ),
)
(
∞
∑
r
∆(r)
S
(r)
k
cos ϑ − ε2− ε2cos [(k − 1)ϑ] + O(ε3 ) ,(2.2)
R(r, ϑ) = Rax ε0 1 − ε1
Rax
Rax
Rax
k=2
(
)
∞
∑
Sk (r)1 r
2
3
sin ϑ − εsin [(k − 1)ϑ] + O(ε ) ,Z(r, ϑ) = Rax ε
(2.3)
Rax
Rax
k=2

1

14th Burning Plasma Simulation Initiative (BPSI) Meeting

H. Seto and M. Yagi

where the subscript “ax” stands a value at the magnetic axis, ∆(r) is the Shafranov shift defined to be∆(0) = 0, Sk (r) is the shaping coeﬃcients, ϑ is the poloidal angle and ζ is the geometrical toroidal angle.
The red colored tags εn show that labeled terms are assumed to be n-th order with respect to ε.
By the use of these orderings, the nonlinear PDE (2.1) can be finally expanded into the following Fourier
series of linearized ordinary diﬀerential equations (ODEs) [5],cos(0ϑ) :cos(1ϑ) :cos(kϑ) :r(2 − s) β 0p(r)I20 (r) + 2 2 += 0, β(r) = 2, B(r, ϑ) = ∇ζ × ∇ψ + I∇ζq Rax
2Bax /2µ0()()
d r3 0
r
r2
2
0
∆
+
R
rβ
−
= 0,
ax
dr q 2Rax
q2
3 − 2s 0
1 − k2Sk00 (r) +Sk (r) +Sk (r) = 0 for k ≥ 2,
r
r2

(2.4)
(2.5)
(2.6)

where 0 stands derivatives with respect to r and I2 (r) is the 2nd-order correction for the flux function I, B
is the magnetic field. Eq.(2.4) and Eq.(2.5) are 1st-order ODEs, where the former describes the radial force
balance and the latter yields the Shafranov shift profile ∆(r). Eq.(2.6) is the 2nd-order ODE determining
the shaping coeﬃcient profile Sk (r). If the safety factor is defined by q(r) = q0 + q1 (r/r1 )l , Eq.(2.6) has
analytic solutions satisfying a boundary condition Sk (a) at r = a and a non-singular condition at r = 0,
rk−1 2 F1 (αkl , βkl ; γkl ; −(q1 /q0 )(r/r1 )l )Sk (r) = Sk (a) k−1,
la
2 F1 (αkl , βkl ; γkl ; −(q1 /q0 )(a/r1 ) )
√
√
k
2
k
2
k2
k2
αkl = −1 + −
−
+
1,
β
=
−1
+
+
− + 1,
kl
2
2
l
l
l
l
l
l

(2.7)
k
γkl = 1 + 2 ,
l

(2.8)

where a is the radial position of the last closed flux surface and 2 F1 (a, b; c; z) is the hypergeometric function.
It is notable that relations between Sk (a) and well-known shaping factors are roughly written by
S2 (a) ' a

1−κ
,
1+κ

S3 (a) ' −

aδ
,
4

(2.9)

where the squareness S4 (r) and the other higher mode coeﬃcients have been neglected.
BOUT++ code uses field-aligned coordinates (x, y, z) [11] based on the local orthogonal flux (LOF)
coordinates (ψ, θ, ζ) rather than the NOT coordinates [4], where θ is the poloidal angle in the LOF coordinates. A coordinate transform from the LOF to the NOT coordinates is, therefore, required to generate
shaped equilibriums with the LOF coordinates by the use of Eqs.(2.2)-(2.9). For the radial coordinateslabel, the relation between ψ(r) and r(ψ) is obtained by solving Eq.(2.10) in the NOT coordinates,
(
)
I
dψ
I
1 ∂Z ∂R ∂R ∂Z=
dϑ
−.
(2.10)dr
2πq
R ∂r ∂ϑ∂r ∂ϑ
The NOT poloidal angle with respect to the LOF coordinates ϑ(ψ, θ) is then evaluated from∫ router∂ϑ∂ϑ∇r · ∇ϑϑ(ψ, θ) = θ −dr,(ψ, θ) =, router = r(ψouter ),∂r∂r∇r · ∇rr(2.11)

where the boundary condition ϑ(ψouter , θ) = θ has been applied implicitly.
To evaluate the impact of shaping eﬀects on NIBMs, we employ the linearized RBM/CBM model [10]∇2⊥ Ak1∂ ∇2⊥ φ1
2
+ b0 × κ0 · ∇p1 ,(2.12)=−B∇k0∂t B02B0∂p1
b0 × ∇⊥ φ1 · ∇⊥ p0=−,(2.13)∂tB0)∂ (Ak1 − d2e ∇2⊥ Ak1 = −∇k φ1 + η∇2⊥ Ak1 ,(2.14)∂twhere t is the time, φ is the electrostatic potential, Ak is the parallel vector potential, κ is the magnetic
curvature, de is the electron skin depth, η is the resistivity. In Eqs.(2.12)-(2.14), the subscript ‘0’ stands theequilibrium part, the subscripts ‘1’ stands the perturbed part of physical quantities normalized with poloidal
Alfvén unit respectively. As described before, the other non-ideal eﬀects are neglected for simplicity.
2

14th Burning Plasma Simulation Initiative (BPSI) Meeting

0.3

0.3

q
α

r/a

(c)

0.2

0.2

0.1

0.1

Z/Rax

Z/Rax

r/a

0.0

(d)

Fig.(a): pressure and safety factor profiles
Fig.(b): pressure and magnetic shear profiles
Fig.(c): shifted circular equilibriumFig.(d): D-shaped equilibrium (κ ∼ 1.3, δ ∼ 0.3)
Fig.(e): resistivity dependence of growth rate ofn=80 CBM/RBM for two geometries

0.0

−0.1

−0.1

−0.2

−0.2

−0.3

−0.3
0.8

3

s

0.5
1.2
0.4 (b)
1.1
0.3
1.0
0.2
C 0.9
0.1
D
0.0
0.8
0.70 0.75 0.80 0.85 0.90

p/pax

1.0
3.4
0.8 (a)
3.2
0.6
3.0
0.4
C 2.8
0.2
D 2.6
0.0
0.70 0.75 0.80 0.85 0.90

0.9

1.0

R/Rax

1.1

H. Seto and M. Yagi

1.2

0.8

0.9

1.0

R/Rax

1.1

1.2

Simulation results and summary

Figs.(a)-(d) show the shifted circular and the D-shaped equilibrium with the same pressure and safety factor
profile based on the LGSE, where we have used following parameters for the both equilibriums, Rax =2.0 [m], a = 0.5 [m], Bax = 2.0[T], 0.7 ≤ r/a ≤ 0.9, βax = 5.0×10−4 , p(r)/pax = 0.5(1−tanh[50(r/a−0.8)]),q(r) = 1.5 + 1.5(r/0.8a)2 . The peak of the pressure shear factor is α = −q 2 R0 β 0 ∼ 0.45 at r/a = 0.8 andthe corresponding magnetic shear factor s = (r/q)q 0 ∼ 1, which is ideally stable since the ideal stability
boundary is αcrit ∼ 0.6s ∼ 0.6. We employ equally spaced grids for 1/80 sector of torus, where the grid
widths in directions perpendicular to the magnetic field line are about 10 times finer than characteristic
lengths relevant to the electron skin depth.
The shaping eﬀects on the n = 80 CBM/RBM growth rate are summarized in Fig.(e). Since theresistivity regime relevant to the edge plasma lies roughly from η ∼ 10−9 for Te ∼ 103 [eV] to η ∼ 10−7 forTe ∼ 102 [eV], the ellipticity and the triangularity can stabilize NIBMs and their stabilizing eﬀects become
stronger in the collisionless regime where CBM is dominant and its growth rate is almost independent of
the resistivity.

Acknowledgement
This work is partly supported by MEXT KAKENHI Grant No. 16K18342. Computations were carried out
on the Helios supercomputer at the Computational Simulation Centre of the International Fusion Energy
Research Centre (IFERC-CSC).

Reference[1] F. Wagner et al., Phys. Rev. Lett. 49, 1408 (1982), [2] N. Asakura et al., Nucl. Fusion 53, 123013(2013), [3] M. Yagi et al., J. Phys. Soc. Jpn. 66, 379 (1997), [4] B.D. Dudson et al., J. Plasma Phys. 81,
365810104 (2015), [5] J.P. Graves, Plasma Phys. Control Fusion 55, 074009 (2013), [6] B.A. Carreras etal., Phys. Rev. Lett. 50, 503 (1983), [7] L. Garcia et al., Phys. Plasmas 6, 107 (1999), [8] F. Porcelli,Phys. Rev. Lett. 66, 25 (1991), [9] J.W. Conner, Plasma Phys. Control. Fusion 35, 757 (1993), [10] H.
Seto et al., Plasma Fusion Res. 11, 1203122 (2016), [11] M.A. Beer et al., Phys. Plasmas 2, 2687 (1995).

3

Particle simulation of plasma heating and production
by electron cyclotron waves in tokamak plasmas
Tomohide Ikeda, Atsushi Fukuyama
Graduate School of Engineering, Kyoto University

1

PurposePlasma production by electron cyclotron waves has been conducted in tokamak, though the physics
mechanism has not been understood well. Various phenomena are coupled with each other - electron
accelation at the elecron cyclotron resonance, ionization by collisions with neutral particles, current
drive by toroidal electric field and electron cyclotron waves, suppresion of particle losses by verticalmagnetic field, and production of closed magnetic configuration. In order to analyze these pheomena
self-consistently, we are developing particle simulation model.

2

Simulation model

Based on the two-dimentional electrostatic particle code pic2des p[1], we have extended the velocity
space to three dimention and make it possible to describe electromagnetic waves by Maxwell’s equation.
The procedure of simulation is as follows: set position and velocity of particles, calculate the current
density on grid, solve normalized Maxwell’s equation to obtain the electric field E and the magnetic
field B at the particle position, and solve relativitic equation of motion. Furthremore, to avoid thereflection of electromagnetic waves at boundary, no-reflection boundary condition and vacuum region
near the boundary are assumed. In the following simulation, time and space is normalized by typicalplasma frequency and Debye length respectively, and the ratio of the light speed and the electron
thermal speed is assumed to be 10.

3

Simulation result

First, we show the propagation of electromagnetic waves excited by waveguide in vacuum. The twodimentional structure of electric field Ey is shown in Fig 1. Horizontal axis and vertical axis are in
the x and y direction respectively. No reflection of waves at boundary was confirmed.

y

Ey
500

Ey
500

Ey
500

Ey
500

Ey
500

400

400

400

400

400

300

300

300

300

300

200

200

200

200

200

100

100

100

100

100

0

0

0

0

0

100

200

x

300

t̄ = 16

400

500

0

100

200

300

t̄ = 32

400

500

0

100

200

300

t̄ = 48

400

500

0

100

200

300

t̄ = 60

400

500

0

0

100

200

300

400

500

t̄ = 72

Fig 1: The propagation of electromagnetic waves in vaccum without reflection at the boundary

Next, linear plasma density is assumed and the reflection at the cutoﬀ where the plasma frequency
equal to the wave frequency. In a plasma the average number of electrons and ions per cell are 32,
and it is assumed that the density increase linearly in the x direction, with n̄(x = 0) = 0.5 andn̄(x = 500) = 1.5, and the wave frequency ω̄ is 1.0. Reflection of the EM wave at the cutoﬀ n̄ = 1 isshown in Fig 2.

1

Ey
500

Ey
500

Ey
500

Ey
500

Ey
500

400

400

400

400

400

300

300

300

300

300

200

200

200

200

200

100

100

100

100

100

0

0

0

0

0

100

200

300

400

500

0

100

200

300

400

500

0

100

200

400

500

0

100

200

300

400

0

500

0

100

200

300

400

500

t̄ = 100t̄ = 80t̄ = 60t̄ = 40t̄ = 20300

Fig 2: Reflection of electromagnetic waves at the plasma cutoﬀ

Finally, static magnetic field which incleases linearly in the x direction is introduced, and X mode
waves are excited from the weak magnetic field side. The magnitude of the static magnetic fieldB¯y (0) = 4.0, B¯y (500) = 8.0 and the wave frequency ω̄ = 6.0. Cylotron resonance is located atx ∼ 250. Time evolution of the spatial profile of wave’s electromagnetic field is shown in Fig 3, andthe electron temperature is shown in Fig 4.
Ez
500

Ez
500

Ez
500

Ez
500

Ez
500

400

400

400

400

400

300

300

300

300

300

200

200

200

200

200

100

100

100

100

100

0

0

0

0

0

100

200

300

400

500

0

100

t̄ = 14

200

300

400

500

0

100

t̄ = 28

200

300

400

500

0

100

t̄ = 42

200

300

400

0

500

0

100

t̄ = 56

200

300

400

500

t̄ = 70

Fig 3: Time evolution of spatial structure of wave electromagnetic field : absorption by cyclotron resonance isobserved
1.6

Tperpe(x)1.6

Tperpe(x)1.6

Tperpe(x)1.6

Tperpe(x)

1.6

1.2

1.2

1.2

1.2

1.2

0.8

0.8

0.8

0.8

0.8

0.4

0.4

0.4

0.4

0.4

0.0
0

100

200

300

400

t̄ = 14

500

0.0
0

100

200

300

400

t̄ = 28

500

0.0
0

100

200

300

400

t̄ = 42

500

0.0
0

100

200

300

400

t̄ = 56500

Tperpe(x)

0.0
0

100

200

300

400

500

t̄ = 70

Fig 4: Electron temperature averaged over y : heating by cyclotron resonance is observed

In order to reduce the discrete noise due to electromagnetic fluctuation, it is necessary to increase thenumber of particles per cell. Furthremore, to carry out the analytis with higher resolution, larger scaleof calculation is needed. Therefore we have implumented high speed cauculation by parallelization,
MPI and OpenMP. Using computer cluster with 8 nodes and 64 cores, computation time with numberof proceeses npr in MPI, and number of threads nth in OpenMP are shown in Table 1. Computation
time without any parallelization is 360s, while parallel calculation 14 times higher than the casewithout parallelization is realized. Parallel number with 2 process and 4 threads per node provides
best perforance. Computation times as a function of the node number is shown in Table 2. It is
confirmed that with the increase of the number of nodes, the computation time is monotonicallyreduced.

2

HH
HH npr
HH
nth
H

1

2

4

8

1
2
4
8

52 s
33 s
26 s
26 s

33 s
26 s
25 s
34 s

28 s
27 s
35 s
46 s

41 s
36 s
49 s
77 s

the number of node
computation times

1
84 s

2
50 s

4
33 s

8
25 s

Table 2: Computation time as a function of the numberof nodes with 2 processes and 4 threads per node

Table 1: Computation time as a function of thenumber of processes npr and the number of threads nth
on one node4

Remaining issues

The remaining issues include to develop physical model, like plasma production by ionization, toconfirm the conversion to the Berstein wave which requires large-scale computation, and to enhancecomputational performance, full-scale parallelization with and introdution of region splitting.
[1] H.Naitou, S,Satake, J. Plasma Fusion Res. Vol.89, No.4 (2013)245-260.

3

Start-up simulation of burning plasma
by integrated modeling code TASK
Kazuki Nagai1 , Atsushi Fukuyama1 , Hideo Nuga2
1

Graduate School of Engineering, Kyoto University
2
National Institute for Fusion Science

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 the time
evolution of momentum distribution functions of fuel ions has been developed by using the
integrated code TASK for tokamak core plasmas.

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 the start-up phase. In order to improvethe momentum resolution and the computational performance, TASK/FP[1] was parallelizedwith respect to the minor radius, the momentum and the particle species.
The purpose of this study is to simulate ITER burning plasma start-up including time
evolution of momentum distribution functions of fuel ions by the use of the integrated code
TASK for tokamak core plasma.

2

Numerical model

The TASK/FP describes the time evolution of momentum distribution functions fs (p, θ, ρ, t)
for particle specie s by solving the bounce averaged drift kinetic equation of the Fokker-Plancktype∂fs= E(fs ) + C(fs ) + Q(fs ) + D(fs ) + Ss ,(1)∂twhere p is the magnitude of momentum, θ the pitch angle, ρ the normalized minor radius,E(fs ) acceleration term by DC electric field, C(fs ) non-Maxwellian Coulomb collision term,Q(fs ) quasi-linear diﬀusion term due to interaction with waves, D(fs ) radial diﬀusion term,and Ss source and sink term (NBI and fusion reaction).
In calculating the Coulomb collision term, the background particle is assumed to be Maxwellian
velocity distribution in conventional Fokker-Planck analysis. In the present analysis, the background particle specie is a non-Maxwellian velocity distribution, the solution of equation (1);therefore this model can describe the collisional eﬀect correctly.
The radial diﬀusion term aﬀects the radial transport of particles and heat. Two transportmodels are used. One is the model which assumes that the radial diﬀusion coeﬃcients areparabolic in space and fixed value in time, and the other is the CDBM turbulent transport1

model[2] which is considered to be more realistic and depends on plasma density and temperature gradient. Since fast ions which have a large cyclotron radius feel spatially averaged
turbulence fluctuation, the radial diﬀusion constant is assumed to be inversely proportional
to the magnitude of momentum. In addition we have reproduced H-mode transport behavior
by reducing the diﬀusion constant in edge region ρ ≥ 0.85.

3Numerical results

We consider a tokamak plasma which consists of electron, deuterium ion, tritium ion,and α particle (helium ion). We assume initial temperature of 2.5 keV, and heating by the
deuteron NBI 50 MW. The fuel ion (deuterium and tritium) is continuously injected into theplasma simulating intermittent pellet injection. Changing the collision model (backgroundMaxwellian distribution or non-Maxwellian distribution) and the radial diﬀusion model, thesimulation was performed in order to look for the operation scenario in which a steady state
Q = 10 is maintained where Q is the power multiplication factor.
Figure 1 shows the 30 seconds time evolution of the power multiple factor for two collisionmodels. Since the background Maxwellian distribution model assumes that the backgroundparticles have a Maxwellian momentum distribution using the temperature evaluated by theaverage kinetic energy which is higher than the bulk temperature, the tritium ions are mainly
heated by collision with NBI ion. Since the existence of fast particles is considered accurately
in the non-Maxwellian distribution model, however the electrons are heated at the beginning.
The large diﬀerence of the power multiplication factor is attributed to the initial collisional
heating.
Figure 2 shows the momentum distribution function of the deuterium ion for every 5 sec.
While the NBI peak can be seen in t = 6 sec and t = 9 sec immediately after heating start, it
is diﬃcult to observe at t = 15 sec because α particle heating dominates. Figure 3 shows thetime evolution of volume averaged density and temperature for the deuterium ion. Figures 4
and 5 show the radial profiles of the density and temperature of the deuterium ion for every
1 second.
The volume averaged density saturates, while the temperature increases exponentially. andthe steady state cannot be achieved with constant heating power due to thermal runaway.
Therefore it is conclude that it is necessary to optimize and control heating power and fuelsupply using pellets in order to realize steady state.

4Remaining issues

The remaining issues are simultaneous heating simulation of core plasma including IC
heating and EC heating as well as NBI heating, maintaining steady state by controlling
heating power and fuel injection[3], and comparison with the analysis using 1D diﬀusivetransport code TASK/TR.

2

20
16

8

Non-linear collision model

40

Heating time t[s]

12

0.3

8

Temperature 0.2

4

0.1

00

30

25

20

15

10

5

0

0.8

0.0

15sec

-3

10

-40

-30

-20

-10

50

p

12

20

30

20

t=9secDensity nD[×10 m ]

40

0
-50

8

Fig.3.Time evolution of density andtemperature

50

p

4

Heating time t[sec]Fig.1.Time evolution of energy multiple

t=6sec0.4

Density

0

p

10

20

30

40

50

40

30

0.6

0sec

0.4

0.2

0.0
0.0

0.2

0.4

0.6

0.8

1.0

Normalized minor radius ρ

20

10

Fig.4.Radial density profile
-40

-30

-20

-10

50

t=15sec
p

0

p

10

20

30

40

50

40

30

20

10

0
-50

-40

-30

-20

-10

0

10

20

30

40

p
factor
Fig.2.Momentum distribution function

50

40
Temperature TD[keV]

0
-50

15sec

30
20
10 10sec
0 0sec
0.0

0.2

0.4

0.6

0.8

1.0

Normalized minor radiusρFig.5.Radial temperature profile

References[1] H. Nuga, A. Fukuyama, Progress in nuclear science and technology. Volume 2, pp78-84(2011)[2] A. Fukuyama, K. Itoh, S-I. Itoh, M.Yagi, M. Azumi, Plasma Phys. Control. Fusion Volume
37 611-631 (1995)[3] Y. Miyoshi, Y. Ogawa, Plasma and Fusion Res. Volume 9, 1405015 (2014)

320

Temperature TD[keV]

Background Maxwellian collision modelDensity nD[×10 m ]Energy multiplication factor Q

12

Construction of reduced transport model by gyro-kinetic simulation
with kinetic electrons in helical plasmas
S. Toda, M. Nakata, M. Nunami, A. Ishizawa1 , T. -H. Watanabe2 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 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. This is because the turbulent transport induces a large amount of the particle and heatloss in toroidal plasmas. The optimization by improving the turbulent transport due to the plasma
micro-instabilities is expected. Recently, a large number of the gyro-kinetic simulations which isapplied to the turbulent transport have been done in toroidal plasmas. The gyro-kinetic simulation
results in tokamak and helical plasmas have been studied with the experimental observations. It
is well-known that helical systems have the advantage in maintaining a steady state for the longplasma discharge. Gyro-kinetic simulations of helical plasmas need a large number of mesh points
along the field line to catch the helical ripple structure.
The GKV-X code [1] solving the gyro-kinetic equation has been used to examine the ion temperature gradient (ITG) mode and zonal flows in the Large Helical Device (LHD) for studying the
turbulent transport in helical plasmas [2]. The gyro-kinetic simulation with the adiabatic electronis performed for the high ion temperature LHD discharge (shot number 88343 [3]). This reduced
model for the ion heat diﬀusivity is proposed [4] using the GKV-X code for the transport simulation. The reduced model is the function of the linear growth rate for the ITG mode and thezonal flow decay time [5, 6]. It is costly to carry out linear calculations of the growth rate by the
gyro-kinetic simulation at each time step of the dynamical transport code such as TASK3D [7],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 tokamak plasmas, thegyro-kinetic simulations at each time step are globally done in the dynamical transport simulation. How to apply the reduced model of the turbulent heat diﬀusivity for the ITG mode derived
from the gyro-kinetic simulation to the transport code has been shown in helical plasmas with alow computational cost [8]. The calculation by substituting the the linear growth rate integrated
over the poloidal wavenumber space and the zonal flow decay time to the reduced model [4] reproduces the results of the reduced model itself within allowable errors. The ion heat flux by the ITG
mode instability agrees with the experimental results [2, 4]. The transport simulation results using
the adiabatic electron for the Ti profile do not contradict with the experimental observation [8].
However, the gyro-kinetic simulation with the kinetic electrons shows the larger ion heat flux thanthe experimental results in the LHD [9]. The value of the electron energy flux by the gyro-kinetic
simulation with the kinetic electrons is close to that in the experimental results of the LHD [9].
The eﬀect of the kinetic electrons induces the enhancement of the linear growth rate of ITG modes.
On the other hand, the eﬀect of the finite β reduces the linear growth rate of the ITG modes.
To compare the gyro-kinetic simulation results with the experimental results with respect to theturbulent ion heat diﬀusivity, these eﬀect of the kinetic electrons and the finite β on the plasma
instability should be included.
In order to examine the eﬀect of the kinetic electrons, the reduced model of the ion heat diﬀusivity
for the ITG mode is needed to be constructed by solving the gyro-kinetic equation in terms of theelectron in addition to the ion. The finite β eﬀect is included in the gyro-kinetic simulation. The
plasma parameter region of the poloidal wavenumber is studied 0.1 ≤ k˜y ≤ 0.5, where the ITG

mode becomes unstable, where k˜y = ky ρi . The same method with the gyro-kinetic analysis using
the adiabatic electrons [4] is adapted to construct the reduced model of the ion heat diﬀusivity as
the function of the linear growth rate and the characteristic quantity of zonal flows. The nonlinear
gyro-kinetic simulation with the kinetic electrons is performed to derive the value of the ion heatdiﬀusivity.
We examined the turbulence driven by the microinstabilities in LHD plasmas by using the gyrokinetic local flux tube code GKV-X [1]. The gyro-kinetic equation for non adiabatic part is solvedin the flux tube coordinate (x, y, z), where x = a(ρ − ρ0 ), y = aρ0 (q(ρ) − ζ)/q0 and z = θ withthe flux coordinate (ρ, θ, ζ). Here, q0 is the safety factor at the flux surface of ρ0 and a is theminor radius. The density, the electron and the ion temperatures at t = 2.233s in the high-Ti
discharge #88343 [3] are used. The GKV-X simulation with the kinetic electrons is performed forthe three dimensional equilibrium field configuration with R = 3.75m, using the plasma profiles inthe VMEC calculation. The finite β eﬀect is included, where β = 0.3% at ρ = 0.65, for example.
To reduce computer resources and perform the nonlinear GKV-X simulation in the wider plasma
parameter region, we take smaller Fourier modes and the smaller number of the grid points inthis simulation than the case in [9]. The total Fourier mode numbers in the x and y direction are
16 and 16. The total Fourier mode numbers in the k̃x and k̃y direction are 8 and 6 in the region−0.5 ≤ k̃x ≤ 0.5 and 0.0 ≤ k̃y ≤ 0.5, where k̃x (= kx ρi ) is the normalized radial wavenumber. The
total grid number in the parallel direction, the parallel velocity direction and the perpendicularvelocity direction are 256, 64 and 16, respectively. We carry out nonlinear GKV-X nonlinear
simulations at ten radial points between ρ = 0.46 and ρ = 0.80. The time evolution of the total
ion energy flux Qi at ρ(= r/a) = 0.65 is examined as the results of the gyro-kinetic simulation
with the kinetic electrons. The value√of time t is normalized by R/vti , where R is the major radiusand vti is the ion thermal velocity, Ti /mi . The saturation of the ion energy flux is obtained in
the nonlinear phase. The averaged value of the ion energy flux in the time interval 50 < t < 100 is
0.12 MW/m2 . In ref. [9], the averaged value of the ion energy flux in the time interval 50 < t < 80
is about 0.13 MW/m2 , when the Fourier mode numbers in the k̃x and k̃y direction are 256 and 43.
Even if the number of the Fourier modes is only 8 and 6 in the k̃x and k̃y directions, respectively,
the close value of the ion energy flux can be obtained in this study. The time evolution of theto the total ion heat flux Qi is also studied. The
ratio of the electromagnetic contribution Qemiem
averaged value of |Qi |/Qi in the time interval 50 < t < 100 is 0.33%, because of the low beta
3% plasma. The saturation of the electron heat flux is seen and the averaged value of the ion
energy flux in the time interval 50 < t < 100 is 0.070 MW/m2 . The averaged value of the electron
heat flux in the interval 50 < t < 80 is about 0.058 MW/m2 . However, the time averaged value
in the time window 50 < t < 100 of the ratio of the electromagnetic contribution to the total
electron heat flux becomes about 30%, even if the low beta 3% plasma. If the larger numberof the modes is taken, the electromagnetic part becomes much smaller [9]. Therefore, we focus
on the construction of the reduced model for the ion heat diﬀusivity
in ⟨this article.
⟩ )The time(
evolutions of the squared turbulent potential fluctuation, T(

squared zonal flow potential, Z = Σkx⟨

ϕ̃kx ,ky =02⟩)= Σkx ,ky ̸=0ϕ̃kx ,ky2/2

and the/2 are obtained, where ϕ̃ is the normalized

electrostatic potential fluctuation which is defined as ϕ̃ = ϕ/(Ti ρi /(eR)). Here, the bracket ⟨⟩represents the averaged values along the magnetic field line. The nonlinear saturation is shown in
the time evolution of T and Z.
Based on the temperatures and density radial profiles and field configuration from the LHD experimental results of high-Ti phase at t = 2.233s, the electron and ion temperature gradients, R/LTeand R/LTi , the normalized density gradient R/Ln and the safety factor q radially change. The
radial profiles of the ion heat diﬀusivity χ̄i are compared by the gyro-kinetic nonlinear simulation
with the kinetic electrons and the adiabatic electrons. Here, the bar¯shows the averaged value inthe time interval 50 < t < 100 of the nonlinear saturation phase. The values of χ̄i by the simulation
with the kinetic electrons are from two to three times larger than those by the simulation withthe adiabatic electrons. The normalized electron and ion temperature gradients are artificiallyaltered from 0.8 to 1.2 times the experimental values at ρ = 0.46, ρ = 0.50, ρ = 0.54, ρ = 0.58,
ρ = 0.62, ρ = 0.65, ρ = 0.68, ρ = 0.72, ρ = 0.76 and ρ = 0.80. The transport coeﬃcients χ̄i /χGBi
are compared with a model function of T̄ and Z̄ in figure 1. A fitting function [4] for the ion heatdiﬀusivity is defined byχ̄i
C1 T̄ α
=
F(
T̄
,
Z̄)
≡
(1)
,
1
χGB
C2 + Z̄ 2 /T̄
iwith α = 0.24, C1 = 0.11 and C2 = 1.8 × 10−3 , where χGBis the gyro-Bohm diﬀusivity. Thei
GB
relative error for fitting χ̄i /χi by F is 0.31, in spite of the wide parameter range[ in the nonlinear]
simulation, where the relative error is defined as the root mean square of (χ̄i /χGBi )/F − 1 .
Therefore, the value of the ion heat diﬀusivity is well reproduced by the model function eq. (1).
The gyrokinetic simulation with the kinetic electrons is performed toconstruct the transport model for the turbulent ion heat diﬀusivity.
The values of the ion heat diﬀusivity are derived from the nonlinear
gyro-kinetic simulation for high-Ti mode in the LHD, where the ITG
mode is destabilized. The nonlinear simulation results for the ionheat diﬀusivity with the kinetic electrons are compared with thoseby the simulation with the adiabatic electrons. The model function
for the ion heat diﬀusivity is shown in terms of the squared turbulent electrostatic potential fluctuation and the squared zonal flow
potential. The reduced model for the ion heat diﬀusivity will be
constructed by the electromagnetic gyro-kinetic simulation with kinetic electrons. The reduced transport models for the electron heatdiﬀusivity and the particle diﬀusivity in addition to ion heat diﬀusivity will be studied when the total Fourier mode number increases.
These is left for future study. This work was partly supported by
JSPS KAKENHI Grant Number 23561002, the NIFS Collaboration
Research program (Plasma Simulator), NIFS16KNST093 and theCollaborative Research Program of Research Institute for Applied
Mechanics, Kyushu University.

Figure 1: Comparison ofχ̄i /χGB from the nonlinear gyro-kinetic simulation
with the model function
F(T̄ , Z̄).

[1] T. -H. Watanabe and H. Sugama Nucl. Fusion 46 24 (2006)[2] M. Nunami, T. -H. Watanabe, H. Sugama and K. Tanaka, Phys. Plasmas 19 042504 (2012)[3] K. Tanaka, et al., Plasma Fusion Res. 5 S2053 (2010)[4] M. Nunami, T. -H. Watanabe and H. Sugama, Phys. Plasmas 20 092307 (2013)[5] H. Sugama and T. -H. Watanabe, Phys. Plasmas 13 012501 (2006)[6] S. Ferrando-Margalet, H. Sugama and T. -H. Watanabe Phys. Plasmas 14 122505 (2007)[7] M. Yokoyama, et al., Plasma Fusion Res. 7 2403011 (2012)[8] S. Toda et al., Journal of Physics: Conference Series 561 012020 (2014)[9] A. Ishizawa et al., Nucl. Fusion 55 043024 (2015)