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

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

11 月 29 日(水)9:30 – 9:40 はじめに
村上 (京大)9:40 – 9:45 事務連絡糟谷 (九大)(座長：糟谷)9:45 – 10:10 講演 1-1
登田 (核融合研)
Reduced models for electron and ion heat diffusivities by gyro-kinetic
analysis with kinetic electrons in helical plasmas
10:10 – 10:35 講演 1-2
山口 (核融合研)
Integrated transport simulation of LHD plasmas with transport models
based on GK simulations and experiment database
10:35 – 10:55 休憩
10:55 – 11:20 講演 1-3
村上 (京大)Simulation study of fusion reaction rate and comparison withexperimental results in LHD
11:20 – 11:40 講演 1-4
前田 (京大)Integrated transport simulation of LHD plasma by TASK3D
11:40 – 11:55 講演 1-5
山本 (京大)Effects of electron cyclotron heating on the toroidal flow in helical plasmas
11:55 – 13:00 昼休み(座長：本多)13:00 – 13:25 講演 2-1
滝塚 (阪大)PARASOL simulation of ELM heat flux in tokamak SOL-divertor plasma
13:25 – 13:50 講演 2-2
東郷 (筑波大)、代理発表 滝塚
Effect of parallel viscosity approximation in inhomogeneous magnetic
fields on SOL-divertor plasma profiles1

13:50 – 14:15 講演 2-3
伊庭野 (阪大)
Simulation of vapor shielding at a solid wall by a weighted PIC method
14:15 – 14:40 講演 2-4
大澤 (九大)Simulation for stability of di-vacancy in tungsten
14:40 – 15:00 休憩(座長：福山)15:00 – 15:25 講演 2-5
本多 (量研機構)
Application of genetic algorithms to modelings of plasma physics
15:25 – 15:50 講演 2-6
Hahm (Seoul Univ.)
Effects of 3D compressional magnetic fields on zonal flows in tokamak
plasmas15:50 – 16:15 講演 2-7
小菅 (九大)Pattern selection in drift wave turbulence with parallel flows
16:15 – 16:40 講演 2-8
佐々木 (九大)Suppression and enhancement of turbulence by geodesic acoustic modes
16:40 – 16:55 講演 2-9
大野 (九大)Analysis of parameter dependency of ion-temperature-gradient instability
using a gyro-fluid model in linear devices
16:55

散会

19:00 – 22:00 懇親会 (博多前炉ばた 一承にて)2

11 月 30 日(木)9:30 – 9:35 事務連絡(座長：村上)9:35 – 10:00 講演 3-1
福山 (京大)Progress in development of module interface of integrated code
10:00 – 10:25 講演 3-2
糟谷 (九大)On effects of impurity and neutral in magnetically confined plasmas
10:25 – 10:45 休憩
10:45 – 11:10 講演 3-3
沼波 (核融合研)
Simulation studies on transport of multi-species plasmas including
impurity ion
11:10 – 11:35 講演 3-4
矢木 (量研機構)
Simulation study on internal transport barrier formation using gyro fluid
model (II)11:35 – 11:50 まとめと議論 糟谷 (九大)Summary and discussion
11:50 – 13:00 昼休み
13:00 – 17:00 核融合エネルギーフォーラムサブクラスター会合
17:00

散会

3

Reduced models for electron and ion heat diﬀusivities by gyro-kinetic
analysis 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. A large number of the gyro-kinetic simulations which are applied to the turbulent transport have been performed in toroidal plasmas. The gyro-kinetic analysis results in tokamak andhelical plasmas have been studied with the experimental observations. Gyro-kinetic simulationsof helical plasmas require a large number of mesh points along the field line in order to capture
the helical ripple structure. Therefore, the gyro-kinetic analysis in helical plasmas consumes thelarger computer resources than for tokamaks. The reduced model, which reproduces the nonlinear
gyro-kinetic analysis results, is needed for the transport simulation in helical plasmas.
The GKV code [1] has been used to examine the ion temperature gradient (ITG) mode andzonal flows in the Large Helical Device (LHD) for studying the turbulent transport. The gyrokinetic simulation with the adiabatic electron is performed for the high ion temperature LHD
discharge (shot number 88343 [2]). The ion energy flux by the ITG mode instability agrees with theexperimental results for the high-Ti plasmas at t = 2.2s in the high-Ti LHD discharge #88343 [3].
The reduced model for the ion heat diﬀusivity is proposed [3] by the simulation with the adiabatic
electron. This reduced model is the function of the linear growth rate for the ITG mode and thezonal flow decay time [4]. How to apply the reduced model of the turbulent ion heat diﬀusivity
from the gyro-kinetic simulation with the adiabatic electron to the transport code has been shownin helical plasmas [5]. The simulation with the kinetic electron shows the larger ion energy flux
than the experimental results in the high-Ti plasmas [6]. On the other hand, the electron andion energy fluxes of the simulation results with the kinetic electron are close to those of theexperimental results in the low-Ti #88343 plasmas [7] at t = 1.9s in the LHD discharge #88343.
To show the reduced transport model for the turbulent electron and ion heat diﬀusivities, theeﬀect of the kinetic electron on the plasma instability has been included. If the small number for
the wave-length is taken, the reduced model for the ion heat diﬀusivity was proposed. At thattime, the electromagnetic part in the total electron heat flux becomes large around 30%, in spiteof the low beta plasma studied [8]. In this study, the reduced models of the electron and ion
heat diﬀusivities for the ITG mode are proposed by solving the gyro-kinetic equation in termsof the electron in addition to the ion to examine the eﬀect of the kinetic electron. The similar
method with the gyro-kinetic analysis using the adiabatic electron [3] is adopted. The nonlinear
gyro-kinetic simulation is performed to evaluate the electron and ion heat diﬀusivities, when thenumber of the wave-length taken in the simulation is larger than the study to estimate only theion heat diﬀusivity [8]. The eﬀect of the kinetic electron on the linear response of zonal flows isstudied. The linear gyro-kinetic simulation is also performed for the reduced model in order toreproduce the values of the heat diﬀusivities by the nonlinear analysis.
The turbulence driven by the microinstabilities is studied in LHD plasmas, using the gyro-kinetic
local flux tube code GKV [1]. Based on the temperatures radial profiles, density radial profiles,
and field configuration from the LHD experimental results of the high-Ti plasmas at t = 2.2s andof the low-Ti plasmas at t = 1.8s, 1.9s [2], the electron and ion temperature gradients, R/LTe andR/LTi , the density gradient R/Ln , and the safety factor q radially change. The GKV simulation

with the kinetic electrons is performed for the three dimensional equilibrium field configurations
with R = 3.75m for the high-Ti plasmas and with R = 3.6m for the low-Ti plasmas. The β value
is 0.3% at ρ = 0.65 for the high-Ti plasmas. To reduce the electromagnetic part in electron energy
flux and the particle flux, we take the larger number for the Fourier modes and the larger numberof the grid points for the parallel direction in this simulation than in the study [8]. The total
Fourier mode numbers in the k̃x and k̃y directions are 41 and 12, where k̃x and k̃y are the radialand poloidal wavenumbers. The total grid number in the parallel direction, the parallel velocity
direction and the perpendicular velocity direction are 512, 64 and 16. In the high-Ti plasmas, thenonlinear GKV simulations are performed at ten radial points between ρ = 0.46 and ρ = 0.80.
The nonlinear simulations are also carried out at the ten radial points between ρ = 0.65 andρ = 0.80 in the low-Ti plasmas. The time evolution of the electron and ion energy fluxes, andthe particle flux at ρ = 0.68 for the high-Ti plasmas is examined. The saturation of the electron
and ion energy fluxes and the particle flux is obtained in the nonlinear phase. The temporalem
traces for the ratio of the electromagnetic contributions Qemand Γem to the electron and ion
e , Qienergy flux Qe , Qi and the particle flux Γ are examined. The time averaged values of Qem
e /Qe ,emem
Qi /Qi and Γ /Γ are 6.3%, −0.76% and 8.1% at the low-β plasma (≃ 0.3%) case. Because
the larger number of the modes is taken, the ratio of the electromagnetic part to the electronenergy flux significantly becomes smaller than the simulation result [8]. The averaged values ofthe electron and ion energy flux at ρ = 0.65 in the time interval 50 < t < 100 for the high-Ti
plasmas are 0.036 MW/m2 and 0.11 MW/m2 . The time averaged value of the particle flux isaround −1.0 × 1019 /(m2 s). In the simulation results [6], the averaged values of the electron andion energy flux in the time interval 50 < t < 80 are about 0.046 MW/m2 and 0.13 MW/m2 , when
the Fourier mode numbers in the k̃x and k̃y directions are 169 and 43. The time averaged value ofthe particle flux is around −1.0 × 1019 /(m2 s). Even if the number of the Fourier modes is smaller
than the simulation results [6], the close values of the electron and ion energy flux and the particle
flux
potential
fluctuation,⟩ )⟨( are obtained⟨ in this⟩study.
) The time evolutions of the squared turbulent(T= Σk̃x ,k̃y ̸=0ϕ̃k̃x ,k̃y2

/2 and the squared zonal flow potential, Z = Σk̃x

ϕ̃k̃x ,k̃y =0

2/2 arestudied. The nonlinear saturation is seen in the time evolution of T and Z.
The nominal values of the LHD experimental results at t = 2.2s for high-Ti plasmas and at
t = 1.8s, 1.9s for the low-Ti plasmas are used. The radial profiles of the electron and ion energy
flux, and the particle flux are examined as the nonlinear simulation results. In this study, thereduced model of the particle diﬀusivity is not shown, because the stationary condition of the
turbulent particle flux with the neoclassical flux with no particle source is violated at ρ < 0.5 in
the high-Ti plasmas and in the low-Ti plasmas [7]. The transport coeﬃcients χ¯e /χGBand χ̄i /χGBbyiithe nonlinear simulation are compared with model functions of T̄ and Z̄, where the bar¯represents
the averaged value in the time interval of the nonlinear saturation phase. A(fitting function
) for
GB
αe
ξ
the electron heat diﬀusivity is defined by χ¯e /χi = Fe (T̄ , Z̄) ≡ C1e T̄ / C2e + Z̄ /T̄ withis the gyro-Bohm
αe = 0.19, C1e = 6.8 × 10−2 , C2e = 2.1 × 10−2 and ξ = 0.10, where χGB
i
by
F
is
0.079,
wheretherelative
error is defined
diﬀusivity. The relative error[ for fitting χ¯e /χGBii ]GB
as the root mean square of (χ¯e /χi )/Fe − 1 . A fitting function for the ion heat diﬀusivity is(1)

also defined by χ̄i /χGB
= Fi (T̄ , Z̄) ≡ C1i T̄ αi / C2i + Z̄ 2 /T̄ with αi = 0.41, C1i = 0.13 andiby Fi is 0.15. The exponent of Z in Fi is
C2i = 4.9 × 10−2 . The relative error for fitting χ̄i /χGBifound to be larger than that in Fe , therefore the eﬀect of zonal flows in χi is stronger than that inχe . The values of the electron and ion heat diﬀusivities are well reproduced by the model functions.

To represent the contributions of the turbulent fluctuations and zonal flow fluctuations to
the turbulence level by the linear gyro-kinetic
analysis, the simulation is performed with the
kinetic electron. The
profiles )the mixing( radial∫
length estimate L ≡ (γ̃k̃y /k̃y2 )dk̃y typicallyintegrated over 0.05 ≤ k̃y ≤ 0.1 with the adiabatic and the kinetic electrons are studied for Figure 1: The comparison of the nonlinear simulation
the high-Ti plasmas, where γ̃k̃y is the linear results, χ̄e and χ̄i with the model predictions.
growth rate. The turbulence fluctuation T̄ isapproximated by T̄ = CT La with the coeﬃcients CT = 6.6 × 10 and a = 1.6. The level of the turbulence is determined by the interaction
between the
The linear zonal flow response function is defined by⟨ turbulence⟩ and⟨ the zonal flows.
⟩
Rk̃x (t) ≡ ϕ̃k̃x ,k̃y =0 (t) / ϕ̃k̃x ,k̃y =0 (t = 0) . The linear zonal flow response is examined using thefield configuration for the high-Ti and low-Ti phases in the LHD #88343 discharge. Note that
we set the value of k̃x around 0.25, because there is a peak of the wavenumber spectra aroundk̃x = 0.25
as the results of the nonlinear simulation. The zonal flow decay time is defined by∫ τfτZF ≡ 0 dtRk̃x (t). Here, we set τf = 30R/vti , because the zonal flow decay time does not change
significantly for τf > 30R/vti . The correlation time of the turbulence in the nonlinear simulation
result is shorter than 30R/vti , thus Rk̃x (t) for τf > 30R/vti does not influence the excited zonal
flow level. The zonal flow fluctuation Z̄ is approximated by the linear simulation results with therelation, Z̄ b /T̄ c = Cz τ̃ZF with Cz = 0.91, b = 0.16 and c = 0.27.
The reduced models which represent the electron and ion heat diﬀusivities in terms of the lin=/χGB
ear plasma
electron are) shown as χRMie( parameters by) the simulation with the kinetic
(
B
B
B1e
B3e
B3i
RM
GB
B1i
2e
2i
A1e L / A2e + τ̃ZF /L, where the coeﬃcients, and χi /χi = A1i L / A2i + τ̃ZF /L
are given by A1e = C1e CT e
Cz−ξ/b = 1.3 × 10, A2e = C2e CT
Cz−ξ/b = 2.0, A1i =1−c/(2b) −1/(2b)α +1−c/(2b) −1/(2b)= 1.8 × 10. The exponents= 2.6 × 102 and A2i = C2i CT
Cz
C1i CT iCz
are given by B1e = αe a = 0.30, B2e = ξ/b = 0.62, and B3e = a(1 − cξ/b) = 0.63, B1i = αi a = 0.66,
B2i = 1/(2b) = 3.1 and B3i = a(1 − c/(2b)) = 0.26. By the linear gyro-kinetic simulation, the values of χRMand χRM
are obtained. The electron diﬀusivity χRMand χ̄e are compared in Fig. 1(a)eieRMand the ion diﬀusivity χi and χ̄i are also compared in Fig. 1(b). The circles and the boxes show
the results in high-Ti and low-Ti phases. The reduced models reproduces the nonlinear simulation
results χ¯e for the relative error 0.21 and χ̄i for the relative error 0.20.
This work was partly supported by the NIFS Collaboration Research programs (Plasma Simulator), NIFS16KNST093, NIFS16KNTT035 and the Collaborative Research Program of ResearchInstitute for Applied Mechanics, Kyushu University, 29FP-4.
α +1−cξ/b

1−cξ/b[1] T. -H. Watanabe and H. Sugama Nucl. Fusion 46, 24 (2006)[2] K. Tanaka, et al., Plasma and Fusion Research 5, S2053 (2010)[3] M. Nunami, T. -H. Watanabe and H. Sugama, Phys. Plasmas 20, 092307 (2013)[4] H. Sugama and T. -H. Watanabe, Phys. Plasmas 13, 012501 (2006)[5] S. Toda et al., Journal of Physics: Conference Series 561, 012020 (2014)[6] A. Ishizawa et al., Nucl. Fusion 55, 043024 (2015)[7] A. Ishizawa et al., Nucl. Fusion 57, 066010 (2017)[8] S. Toda et al., Plasma and Fusion Research 12, 1303035 (2017)

Study of Energetic Particle Confinement by Triton Burn-up
Experiment in LHD
S. Murakami1 , M. Isobe2,3 , K. Ogawa2,3 , T. Nishitani2 , S. Kamio2 , Y. Fujiwara2 , N. Pu3 ,H. Kawase3 , M. Osakabe2,3 , Y. Saito1 , S. Maeta1 , H. Yamaguchi2 , LHD Exp group2
1 Department of Nuclear Engineering, Kyoto University,
2 National Institute for Fusion Science, National Institutes of Natural Sciences
3 School of Physical Sciences, SOKENDAI (The Graduate University for Advanced Studies)e-mail: murakami@nucleng.kyoto-u.ac.jp
The energetic particle confinement is investigated by triton burn-up experiment in LHD. A significant increase of triton burn-up ratio can be seen in the strongly inward shifted configuration, in
which the energetic particle confinement becomes nearly optimum in the LHD configurations. The
triton burn-up simulation is also performed by GNET code, and we obtain a similar increase of thetriton burn-up ratio depending on the magnetic configurations in LHD.
Introduction
The deuterium plasma experiment has been started from 2017 campaign in the Large Helical
Device (LHD)[1]. The study of the energetic particle is one of the important issues in the deuterium
plasma experiment of LHD. The nuclear fusion between deuterium ions produces the 2.5MeV neutron
or 1MeV triton almost equal probability. Also, the produced 1MeV tritons react with the deuterium
and produce 14MeV neutrons. The D-T fusion cross section has a peak at the energy of about
150keV and, then, the 1MeV tritons would react mostly after slow down to that energy. Therefore,
the confinement of 1MeV tritons can be seen by measuring the produced neutrons[2].
Because of the three-dimensional magnetic configuration, the behaviors of trapped particles arecomplex, and the confinement of the energetic particle is an important issue in the heliotrons. The
LHD has flexibility in the magnetic configuration controlling the poloidal magnetic field, etc. It is
pointed out that the inward shift of the magnetic axis position in major radius, Rax , improves theconfinement of the energetic particles[3]. We have to confirm this fact experimentally to design the
future fusion reactor in heliotron configurations. The confinement of energetic particle has been
investigated in a hydrogen operation phase of LHD. The charge-exchange neutral particle analyzer(NPA) and scintillator-based fast-ion loss detector (FILD) have played a primary role in the measurement of energetic particles[4]. However, we can obtain an only limited physics information fromthese diagnostics. NPA can detect energetic particles having a limited pitch-angle along the line ofsight and strongly depends on the density profile of neutral particles. We need information of theenergetic particle in the important core region of plasma.
Triton burn-up experiment
Time-resolved 14 MeV neutron flux is measured by scintillatingfiber detectors. More than 20MW of NBI heating power is injected
into deuterium plasma, and we obtain the electron and ion temperatures more than 3keV with the plasma density about 3×1019 m−3 .
Figure 1 shows the typical time evolutions of total and 14MeV neutron rates. We can see that a build-up rate of 14 MeV neutron flux isslower than that of total neutron rate. It comes from cross-section
curve for D-T reaction and the peak of the triton fusion reaction
starts after the energy slow down about 150keV.
Triton burn-up ratio is evaluated in stellarator/heliotron devices
for the first time. It is defined as a ratio of secondary DT neutron
yield to total neutron yield. The scintillating-fiber detector is cal- Figure 1: Typical time evolutions of total and 14MeV neutron
ibrated by using results measured with calibrated NAS. Therefore, rates.
triton burn-up ratio can be evaluated in every discharge in LHD.

1

Confinement property of helically trapped energetic ions largely depends on magnetic field configurations. The triton burn-up ratios are plotted as a function of Rax (3.55m, 3.6m, 3.75m, and3.9m) in Figure 2. A significant increase of triton burn-up ratio can be seen in the strongly inward
shifted configuration. It is pointed out that the inward shift of the magnetic axis position in major
radius, Rax , improves the confinement of the energetic particles[3]. Figure 3 shows the typical orbitof trapped particle in the Rax =3.75m (left) and 3.55m (right) configurations. We can see a nearly
optimum of confinement of the trapped particle in the Rax =3.55m configuration, where the large
burn-up ratio is obtained.

[1] M. Osakabe, et al., Fusion Sci. Technl. 71, 199 (2017).
[2] M. Isobe, et al., Rev. Sci. Instrum. 81 (2010) 10D310.
[3] S. Murakami et al., Nuclear Fusion 42, Nuclear Fusion
42 (2002) L19.

Figure 2: Triton Burn-up ratio as afunction of Rax (3.55m, 3.6m, 3.75m,
and 3.9m).

Figure 3: Typical orbit of trapped particle in the Rax =3.75m (left) and 3.55m(right) configurations.

v⊥/v1MeV

Simulation analysis
We perform the triton burn-up simulation of the deuterium
experiment of LHD and evaluate the D-T fusion reaction rates
to compare with the experimental results of the 14 MeV neutron diagnostic system. The radial profile of the 1 MeV triton
production rate due to the D-D fusion reaction between Dbeam and D-thermal ions is evaluated by GNET. As a result,
we obtain a center-peaked triton production rate.
Next, we solve the 5-D drift kinetic equation for the 1MeV
tritons by GNET and evaluate the slow down distribution
in the real and velocity space. Figure 4 shows the slowing
down velocity space distribution of the 1MeV tritons in LHD(Rax = 3.6m). We can see a relatively large loss region near
the boundary of the passing and the helically trapped particles. About 25% of tritons escape with almost the initialenergy of 1 MeV by the prompt orbit loss due to drift motion
immediately after their birth (t < 10−5 s). After the prompt
orbit loss, the collisionless diffusive loss (10−5 < t < 10−2 )and, then, the collisional diffusive loss (t > 10−2 s) become
dominant.
It is found that the confinement of the 1MeV tritons is improved by the strongly inward shifted configuration of LHD(Rax = 3.5m) and that the triton burn-up ratio is increased
to about 0.1% (Fig. 2). This tendency agree with the experimental results and it well agrees also quantitatively if we add
a constant factor 0.14% as an additional effect.

28
26
24
22
20
18
16
14

1
0.8
0.6
0.4
0.2
0
-1

-0.5

0
v||/v1MeV

0.5

1[4] M. Isobe et al., Fusion Science and Technology 58, 426Figure 4: Slowing down velocity space
435, 2010.
[5] K. Ogawa et al., Nuclear Fusion 53, 053012, 2013.
[6] S. Murakami et al., Nuclear Fusion 46, S425 (2006).
[7] H. Yamaguchi and S. Murakami, Nucl. Fusion 56 (2016)
026003.

2

distribution of the 1MeV tritons in
LHD (Rax = 3.6m).

Integrated transport simulation of LHD plasma by TASK3D
S. Maeta, S. Murakami, H. Yamaguchi, A. Fukuyama, K. Nagaokaa , H. Takahashia ,H. Nakanoa , M. Osakabea , K. Tanakaa , M. Yokoyamaa , and LHD Exp Groupa
Department of Nuclear Engineering, Kyoto University, Kyoto 615-8540, Japana National Institute for Fusion Science,322-6 Oroshi-cho, Toki 509-5292, Japan

introduction
Various physical phenomena with wide ranges of time and space scales are connectedin the fusion plasma. Thus, in order to predict the performance of the fusion plasma, anintegrated simulation combining various physics models is required. Many integrated
simulation code has been developed for tokamak and helical plasmas. TASK3D is an
integrated simulation code for helical plasma and has been applied to the analyses ofthe plasma transport of LHD plasmas[1, 2].
The D-plasma experiment has been started from March 2017 in LHD, where thedeuterium NBI heating beams with the power more than 30MW is injected to the deuterium plasma. Main objects of this project are to improve plasma performance, tomake clear the isotope eﬀect on the energy confinement and/or turbulent transport inhigh T i plasma and to understand energetic ion confinement by D-D nuclear fusion neutron measurements.
In this paper, we verify the turbulent transport model of deuterium plasma in LHD
using TASK3D. We evaluate the NBI heat deposition of deuterium plasma by GNET
code[4, 5]. An integrated simulation of deuterium plasma has been performed assuming a transport model including transport improvement obtained in the H/He plasma
simulation. In addition, we investigate the ECH plasma of LHD applying TASK3D andcompare with the experimental results. We perform the integrated transport simulation
by assuming several turbulent transport models.

Heat transport simulation of deuterium plasma and ECH plasma
In previous studies, hydrogen experiments are well reproduced by assuming gyro-Bohm(gB)
for electron and gB+gradT for ion turbulence transport coeﬃcients in hydrogen experiments[2]and H/He experiment results are well reproduced assuming Ae f f depending modelχTi B = χTi B × exp(−kF(Ae f f − c))(1)where c = 1.4, k = 2.0 and F(x) is step function. We have selected hydrogen anddeuterium plasma with similar density profiles (Fig.1). We compare the temperature
profiles of experiment and simulation results to clarify the isotope eﬀect on the turbulent

1(a) ne(b) Ti(c) TeFigure 1: Experimental result of deuterium or hydrogen plasma.
transport. We investigate the heat transport of hydrogen and deuterium plasma using
TASK3D and GNET.
First, we performed hydrogen plasma simulation assuming gB+gradT for ion andgB for electron turbulence transport model. We determine the coeﬃcient of gB+gradT
model so that the simulation result agree well with experimental results of hydrogenplasmas. Simulation result of hydrogen plasma (Fig.2) shows similar radial profile of T ivalues.
Next, we performe the deuterium plasma simulation. We apply the gB+gradT modeland Ae f f depending models to the deuterium plasma assuming the coeﬃcient of thesemodel same as that of hydrogen plasma. Fig.3, 4 and 5 show simulation results ofdeuterum plasma. Simulation results with Ae f f depending model agree well with theexperimental results of T i profiles in all three density cases.

Summary
We have developed the integrated simulation code TASK3D to predict the plasma performance in the 3D magnetic configuration plasma. We have investigated deuterium
experimental plasma by the integrated simulation code TASK3D. We determine the coeﬃcient of gB+gradT model so that the simulation result agree well with experimental
results of hydrogen plasmas. Simulation results assuming the present turbulent transportmodel based on the H/He plasma good agree with the deuterium experimental results.

References[1] S. Murakami, et al., Plasma Phys. Control. Fusion 57 (2015) 054009.
[2] A. Sakai, et al., Plasma Fusion Res. 9 (2014) 3403124.
[3] K. Nagaoka, et al., Plasma Fusion Res. 11 (2016) 2402106
2[4] H. Yamaguchi and S. Murakami, Plasma Fusion Res. 9 (2014) 3403127.
[5] S. Murakami, et al., Nuclear Fusion 46 (2006) S425.

(a) low density(b) middle density(c) high density

Figure 2: Simulation result of hydrogen plasma.

(a) T i(b) T e(c) χiFigure 3: Simulation result of deuterium plasma of low density.

3

(a) T i(b) T e(c) χiFigure 4: Simulation result of deuterium plasma of midlle density.

(a) T i(b) T e(c) χiFigure 5: Simulation result of deuterium plasma of high density.

4

Eﬀects of electron cyclotron heating on the toroidal flowin helical plasmas
Y.Yamamoto, S.Murakami, H.Yamaguchi1 , Chang Ching-Chieh, H.Takahashi1
1Graduate School of Engineering, Kyoto University
National Institute of Fusion Science, Toki 509-9252, Japan

Abstract
We study the behaviors of supra-thermal electrons by the electron cyclotron heating, which
can generate the j × B torque and the collisional torque in the toroidal direction in the LHD
plasma. We apply the GNET code, which can solve a linearized drift kinetic equation in 5-D
phase space. As a result, we find that the obtained torque by ECH in the toroidal direction iscomparable to that by NBI (neutral beam injection) and that its direction is opposite to NBI
torque in the inner region.

1

Introduction

Recently, spontaneous toroidal flows in the ECH plasma have been observed in many tokamaks andhelical devices e.g. JT-60U, HSX, and LHD. Many studies have been done experimentally[1] andtheoretically[2]. Particularly, in LHD, when ECH was applied into the NBI heated plasma, the radial profile of the toroidal flow velocity changes drastically and the direction is reversed in the coreregion as shown in figure 1. This change of the toroidal flow by ECH has not yet been understoodwell.
In the previous studies, we have found that ECH generates aradial flux of energetic electrons in non-axisymmetric plasmas[3].
Therefore, the radial bulk ion current must flow to cancel thisenergetic electron current, and the bulk ion current could generate
the toroidal torque by the j ×B force. In this study, we investigate
the behaviors of energetic electrons by ECH, which can generate
the radial current making the j × B torque in the LHD plasma.
Also, we evaluate the collisional torques, by the collision between
energetic electrons and bulk plasma. We apply the GNET code,
which can solve a linearized drift kinetic equation for energetic
Figure 1: The toroidal flow
electrons by ECH in 5-D phase space[3]. Then, we calculate thechange in the LHD experitoroidal flow by solving the radial diﬀusion equation including thement.
toroidal component of the j × B torque, collisional torques, andthe NBI torque, which is evaluated by FIT-3D code[4].
50NBI
NBI+ECH30
20

φV [km/s]

40

10

0

-10

-20
-1

-0.5

0

0.5

1

r/a

2

Simulation model

In order to study ECH, we apply the GNET code, which can solve a linearized drift kinetic
equation for energetic electrons in 5-D phase space. We can write the gyrophase averaged electron
distribution function asf (r, v∥ , v⊥ ) = fMax (r, v 2 ) + δf (r, v∥ , v⊥ )
(1)

1

3.5

1.0e-07

3.5

5.0e-08

2.5

3.5

5.0e-072.5

20.0e+001.5
12
0.0e+00
1.5
1-5.0e-08

0.5-5.0e-07

0.5

0-1.0e-07
-3

-2

-1

0
v|| / vthe

125.0e-0820.0e+00
1.5
1-5.0e-08

0.5

0

3

1.0e-073v⊥ / vthev⊥ / vthe

1.0e-063v⊥ / vthe

32.5-1.0e-06
-3

-2

-1

0
v|| / vthe

1

2

0

3

-1.0e-07
-3

-2

-1

0
v|| / vthe

1

2

3

(a) r/a ∼ 0.00:The inner re- (b) r/a ∼ 0.15:The region of (c) r/a ∼ 0.25:The outer region of the deposition.
the deposition.
gion of the deposition.

Figure 2: The velocity distribution of δf
0.15

0.15
collision
j x B
total

collision
j x B
total

2

T [N/m ]0.1

2

T [N/m ]

0.1

0.05

0

-0.05

-0.1

0.05

0

-0.05

0

0.1

0.2

0.3

0.4

r/a

(a) Axisymmetry

0.5

-0.1

0

0.1

0.2

0.3

0.4

0.5

r/a

(b) LHD

Figure 4:
Figure 3: The j ×B and collisional torques by ECH(1MW). torques.

NBI and ECH

where fMax represents Maxwellian, and r, v∥ and v⊥ are the particle position vector, the paralleland perpendicular velocity, respectively. The drift kinetic equation for δf is given by
∂δf
∂δf
∂δf
+ (v d + v ∥ ) ·+ v̇ ·− C(δf ) − L(δf ) = S ql (fMax )
∂t
∂r
∂v

(2)

where v d is the drift velocity. C(δf ), L(δf ) and S ql (δf ) are linear collision operator, orbit loss,
and the quasi-linear diﬀusion operator for the absorption of ECH, respectively.

3

Results

Figure 2 shows the velocity destribution of δf at each minor radius. We can see that moreenergetic electrons can be found outer / inner regions apart from the heating point(r/a ∼ 0.15).
Thus, the radial electron flux is produced and it makes ion return current.
The j × B and collisional torques cancel each other in axisymmetric plasma[5]. We first
calculate the torques in the axisymmetric magnetic configuration. As shown in figure 3(a), we
obtain this cancellation of the toroidal torques. Next, we study the behaviors of energetic electronsand evaluate the torques assuming the LHD plasma. Figure 3(b) shows that the j × B torque islarger than the collisional torque suﬃciently in LHD. Figure 4 shows the torques with experimental
heating parameters. As a result, we find that the obtained torque by ECH is almost same order
as NBI torque and that its direction is opposite to NBI torque in the inner region.
We evaluate the toroidal flow velocity Vϕ solving the radial diﬀusion equation
(
)
∂Vϕ
∂Vϕ
1
1 ∂
rD
+
=
(TECH + TNBI + TNTV )
(3)
∂t
r ∂r
∂r
mi ni R2

where T and D are torques by ECH, NBI or Neoclassical Toroidal Viscosity (NTV) and radialdiﬀusivity of toroidal flow. TNTV is given as
TNTV = −mi ni µti

δB 2⟨RVϕ ⟩
B02(4)

where µti is the viscosity coeﬃcient and we set its value as the order of 105 . In the calculation, we
assume that ECH is superposed into the NBI plasma. Figure 5 shows the change of toroidal flows
by ECH. The j × B torque of ECH makes the toroidal flow reduce / enhance at the inner/outerof heating point. It agrees with the experiment observations qualitatively.
D=0.1, µ=4x105
NBI
NBI+ECH(4MW)
NBI+ECH(6MW)
NBI+ECH(8MW)40

50
Vφ [km/s]

Vφ [km/s]

50

30

40
30

20

20

10

10

0

0

0

0.1 0.2 0.3 0.4 0.5 0.6
r/a

(a)

D=0.05, µ=4.5x105
NBI
NBI+ECH(4MW)
NBI+ECH(6MW)
NBI+ECH(8MW)60

D=0.01, µ=5x105
NBI
NBI+ECH(4MW)
NBI+ECH(6MW)
NBI+ECH(8MW)60
50
Vφ [km/s]

60

40
30
20
10

0

0.1 0.2 0.3 0.4 0.5 0.6
r/a

0

0

0.1 0.2 0.3 0.4 0.5 0.6
r/a

(b)

(c)

Figure 5: Change of toroidal flows by ECH

4Conclusion

We have evaluated the behaviors of energetic electrons by ECH applying GNET code in order tomake clear the mechanism of the toroidal flow change in LHD plasma heated by ECH and NBI.
The obtained j ×B torque is the same order as NBI torque and its direction agree with experimentobservation.

References[1] M.Yoshida, et al., Nucl. Fusion 49 115028 (2009)[2] P. H. Diamond, et al., Nucl. Fusion 53 104019 (2013)[3] S. Murakami, et al., Nucl. Fusion 40 693 (2000)[4] S. Murakami, et al., Fusion Technol. 27 256 (1995)[5] M. N. Rosenbluth and F. L. Hinton, Nucl. Fusion 36 55 (1996)

3Eﬀect of parallel viscosity approximation in inhomogeneous
magnetic fields on SOL-divertor plasma profiles
Satoshi Togo, Tomonori Takizukaa) , Mizuki Sakamoto, Naomichi Ezumi, Yuichi Ogawab) ,Kenzo Ibanoa) , Kunpei Nojiri, Yue Lib) , Yousuke Nakashima
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)
Graduate School of Frontier Sciences, University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa 277-8568, Japana)

The eﬀect of parallel viscous approximation on plasma profiles is studied on a one-dimensional scrape-oﬀ
layer (SOL)-divertor system by comparing the profiles from two plasma fluid models; the Braginskii model andthe anisotropic-ion-pressure (AIP) model. In a collisional plasma, profiles from two models agree fairly wellwhile a clear deviation is obtained in a collisionless plasma. In the latter case, the viscous flux-limiting model
was used beyond its criterion for physical validity which can be important in order to reproduce the behavior of
detached plasmas with the Braginskii model.

1 Introduction
For designing divertors of future fusion devices, it
is expected to develop a simulation code for plasmas
in scrape-oﬀ layer (SOL) and divertor regions which
considers all of important physical processes and iscomputationally realistic. Various SOL-divertor plasma
code packages have been developed such as SOLPS(B2-EIRENE) [1], SONIC [2, 3] and UEDGE [4]. These
code packages can basically reproduce experimental
results fairly well. However, in some cases, especiallydetached plasmas, the agreement between experimental
and numerical results is not satisfactory [5, 6].
Conventional code packages are based on a plasma
fluid model derived by Braginskii [7] and the kinetic
eﬀect, which becomes essential in collisionless conditions,is considered by introducing some flux-limiting models
to the heat and viscous fluxes [8]. For the parallelviscous flux, which is an approximation of the anisotropic
part of the ion pressure, it is possible to exclude it (ormake it automatically considered) by directly introducing
the anisotropic ion pressure (AIP) to the plasma fluid
model. In this paper, the eﬀect of parallel viscous
approximation on plasma profiles is studied by comparing
the profiles from the Braginskii model and the AIP one ona one-dimensional SOL-divertor system.

2 Models
2.1 Anisotropic-ion-pressure (AIP) model
Basic plasma fluid equations of the AIP model are given as
follows [9, 10];∂n
∂ ( nV )
+B
= S,
(1)
∂t
∂s B
author’s e-mail: togo@prc.tsukuba.ac.jp

[
)]
∂ 1(
∂
2
(mi nV) + B
mi nV + pi,∥ + pe∂t∂s Bpi,⊥ + pe ∂B= Mm ,
+
B
∂s
(
)
∂ 1
1
mi nV 2 + pi,∥
∂t 2
2
[ (
)]
3
∂ 1 1
mi nV 3 + pi,∥ V + qi,∥+B
∂s B 22pi,⊥ V + qi,⊥ ∂B+B∂spi,⊥ − pi,∥ me pe − pi,∥∂pe
= Qi,∥ +
,
+
−V
τrlx
mi
τe
∂s
(
)
∂pi,⊥
pi,⊥ V + qi,⊥ ∂B∂ pi,⊥ V + qi,⊥
+B
−
∂t
∂s
B
B
∂s
pi,∥ − pi,⊥ 2me pe − pi,⊥+,
= Qi,⊥ +
τrlx
mi
τe
(
)
[ (
)]
∂ 3
∂ 1 5
pe + Bpe V + qe∂t 2∂s B 2
3me pi − pe∂pe.
= Qe +
+V
mi
τe
∂s

(2)

(3)

(4)

(5)

Here, the pressure is defined by pσ = nT σ in which σ ∈{(i, ∥) , (i, ⊥) , i, e}. The eﬀective-isotropic ion temperature()
is defined by T i ≡ T i,∥ + 2T i,⊥ /3. Other notations are
the same as Ref. [11]. The ion-pressure relaxation time isgiven by τrlx = 2.5τi [12] in which the ion-ion Coulomb
collision time τi is estimated by using T i . The electron
conductive heat flux qe is estimated by a harmonic averageSHof the Spitzer-Härm heat conduction qSHe = −κe√ (∂T e /∂s)and a free-streaming heat flux qFS= nT e T e /me as
e
[
(
)]−1
SH
FS
qe = 1/qe + 1/ αe qe
in which the heat-flux limiting
factor of electron is set to be αe = 0.5. The ion conductive
heat fluxes, qi,∥ and qi,⊥ , are estimated by the Spitzer-Härm()(1/3) κiSH ∂T i,∥ /∂s and qSHheat conductions as qSHi,⊥ =i,∥ = −()− (2/3) κiSH ∂T i,⊥ /∂s , so that the total ion conductive heatSH
SH(∂T i /∂s). A virtual divertor
flux is qi = qSHi,∥ + qi,⊥ = −κimodel is used for the sheath boundary conditions [13, 14].

2.2 B2 code

in which δpi is approximated by the following viscousflux;δpi ≈ −ηi B−1/2∂ ( 1/2 )B V ≡ πi .
∂s(7)Here, the following flux-limiting model is used inestimating the eﬀective parallel-ion viscosity ηi ;()
ηcl ∂V
1 + Ωη ηi = ηcl ≡ 0.96pi τi , Ωη =.
βpi ∂s(8)

In this study, the viscous-flux limiting factor is set to be
β = 4/7 [8]. Note that ηi always becomes a finite positive
value because Ωη is positive or zero.
For the Bohm criterion at the sheath edge, V ≥ cs is
imposed in which V is extrapolated if it tends to be higher√
than the ion sound speed cs ≡ (T i + T e ) /mi .
By summing Eqs. (3) and (4), the equivalent equationof ion energy in the Braginskii model is obtained asfollows;()∂ 1
3
mi nV 2 + pi
∂t 2
2
[ (
)]
5
∂ 1 1
(9)
mi nV 3 + pi V + πi V + qi
+B∂s B 22
3me pe − pi∂pe= Qi +
.
−V
mi
τe
∂s
Here, Qi = Qi,∥ + Qi,⊥ and the ion conductive heat flux qi
is estimated by the Spitzer-Härm heat conduction qSH=i−κiSH (∂T i /∂s). Other equations are the same as those of
the AIP model.

2.3 Calculation conditions
In this study, one-dimensional (parallel-to-B direction
the coordinate of which is denoted by s) plasmas in a
SOL-divertor-like system are simulated by the B2 code

1.5
(b)
B/B0

The B2 code [15] is one of the most famous codes for
the SOL-divertor plasma based on the Braginskii model.
Although the B2 code is a two-dimensional code, it isapplied to a one-dimensional system in this study byneglecting all of radial transports. By rewriting Eq. (2)()
in terms of the isotropic, pi ≡ pi,∥ + 2pi,⊥ /3, and()
anisotropic, δpi ≡ 2 pi,∥ − pi,⊥ /3, parts of the ionpressure instead of pi,∥ and pi,⊥ , the equivalent equationof parallel momentum of plasma in the Braginskii model
is obtained as follows;()∂ mi nV 2∂(mi nV) + B
∂t
∂s
B
(6)
(
)
∂
∂
3/2
−3/2
+B
B δpi = − (pi + pe ) + Mm ,
∂s
∂s

1.0
0.5
0.0
0.0

0.2

0.4

0.6

0.8

1.0

s/L

Fig. 1 (a) Schematic picture of a cross section of a tokamak and(b) normalized magnetic field B/B0 as a function of thenormalized parallel coordinate s/L.

and the AIP model. Figure 1 (a) shows a schematic
picture of a poloidal cross section of a tokamak plasma.
Plasmas in the outer (low-field side) SOL-divertor-likesystem are simulated roughly taking into account thedependence of the magnetic field strength on the major
radius, B ∝ R−1 . Figure 1 (b) shows the parallel-to-B
profile of the normalized magnetic field strength B/B0
actually used in this study. The parallel system length isset to be L = 2.8 m although it is much shorter than a real
tokamak SOL-divertor system. The ratio of the maximum
to the minimum magnetic field strength is Bmax /Bmin =2. As boundary conditions, a mirror symmetry condition(∂n/∂s = ∂T/∂s = V = 0) is imposed at the center s/L = 0
for simplicity and a sheath-edge condition is set at thetarget s/L = 1. In order to examine the viscosity behavior
excluding the source eﬀect, localized source terms forparticle, S , and energy,[ Qi and ]Qe , are artificially given
in proportion to exp −20 (s/L)2 . The ion heat source Qi
is kept to have the same value as the electron one Qe andis divided into Qi,∥ = Qi /3 and Qi,⊥ = 2Qi /3 assuming
isotropic heat sources for the AIP model. The momentum
source Mm is set to be zero. The absolute values ofsource terms are changed in order to control the plasma
collisionality which is estimated by L/λmfp (λmfp represents
the mean free path of ion-ion Coulomb collisions) justin front of the sheath by the AIP model in this study.
Deuterium is chosen as the ion species.

3

Results

Direct comparisons of plasma parameter profiles between
the B2 code and the AIP model are shown in Figs. 2 and

0.0
4n (1018/m3)

0.5

(b)

2

(c)

B2
AIP

M

0
6
3

(d)

5
0
10
5
0
0.0

0.2

0.5
0.0
1.5

0.4

0.6

0.8

0.5
0.0
4

(c)

(d)

20
10

(e)

20
10
0.2

s/L
Fig. 2 Parallel-to-B profiles of (a) normalized magnetic field
B/B0 , (b) plasma density n, (c) Mach number M, (d) ion
temperature T i and (e) electron temperature T e of the B2
code (solid lines) and the AIP model (broken ones) for acollisional plasma (L/λmfp ≈ 20).

In the collisional case, Fig. 2, a fairly good agreement
between two models is obtained in all plasma parameters
except for the vicinity of the target. A supersonic flow due
to a magnetic nozzle [16] in the diverging-magnetic-field(DMF) divertor region s/L > 0.85 is obtained in bothmodels. The deviation in the vicinity of the target mightcome from the diﬀerent treatments of the Bohm criterion
at the sheath entrance.
In the collisionless case, Fig. 3, on the other hand,a clear deviation is seen in the upstream DMF region in
profiles of n and M in addition to the vicinity of the target.
Because the temperature profiles are almost the same inthe upstream DMF region, it is presumed that the parallelviscosity approximation of the anisotropic part of the ionpressure is responsible for this deviation.

0.4

0.6

0.8

1.0

s/L
Fig. 3 Parallel-to-B profiles of (a) normalized magnetic field
B/B0 , (b) plasma density n, (c) Mach number M, (d) ion
temperature T i and (e) electron temperature T e of the B2
code (solid lines) and the AIP model (broken ones) for acollisionless plasma (L/λmfp ≈ 0.2).

4
3 for collisional (L/λmfp ≈ 20) and collisionless (L/λmfp ≈0.2) plasmas, respectively. Here, in estimating the Mach
number M profiles, the parallel√temperature T i,∥ is used for()
the ion sound speed as cs ≡T i,∥ + T e /mi in the AIP√
model while cs ≡ (T i + T e ) /mi is used in the B2 code.

B2
AIP

2

0
0.0

1.0

(b)

1.0

0
30

(e)

(a)

1.0

0
30
Ti (eV)Ti (eV)

B/B0

1.0

0
10

Te (eV)

1.5(a)

Te (eV)Mn (1018/m3)

B/B0

1.5

Discussion

By combining Eqs. (1)-(4), we derive the following
generalized expression for the viscous flux;()
gen
1 + Ωη δpi∂ ( 1/2 )B V= −ηgen B−1/2
( ∂s
)
3V
∂V ∂pi ∂pe+ ηgen
mi nV
+
+
4pi
∂s
∂s
∂s
(
)
1
∂
− ηgen2qi,∥ − qi,⊥(10)
2pi ∂s2qi,⊥ − qi,∥ ∂B− ηgenBpi∂s57+ τrlx mi V 2 S − τrlx V Mm
9
9
(
)
2
+ τrlx 2Qi,∥ − Qi,⊥ ,910pi τi ≈ ηcl ,9
ηgen ∂V ηgen V ∂B 2me τrlx
gen
Ωη =
+
+
,
(4/7) pi ∂s
8Bpi ∂s 3mi τeηgen =

(11)
(12)

(Pa)

0.5
0.0
-0.5

gen

-1.0
1
Ωη

gen

(a)

(b)

0
-1
0.0

0.2

0.4

0.6

0.8

1.0

s/L
(
)
Fig. 4 Parallel-to-B profiles of (a) 1 + Ωηgen δpi (solid line), theRHS of Eq. (10) (broken line), πgen(chain line) and (b)iΩηgen for the collisional case.

(Pa)

20

-20

gen

-40
80
Ωη

(a)

0

however, Ωη becomes smaller than −1 in 0.2 ≤ s/L ≤0.4 as shown in Fig. 5 (b). It means that the eﬀective()
gengen −1
parallel-ion viscosity ηeﬀ = ηgen 1 + Ωηto reproduce
the behavior of δpi in this region becomes negative whichis physically unstable. Therefore, by using the absolutevalue of the flow velocity gradient, |∂V/∂s|, in estimating
the eﬀective parallel-ion viscosity as is done in the B2
code, Eq. (8), the problem becomes solvable. In thiscase, however, the validity of the viscosity approximation
is lost leading to deviated profiles of plasma parametersfrom those from the AIP model as shown in Fig. 3.
The deviation in the density profiles in the upstream
region as shown in Fig. 3 could lead to a diﬀerencein the amount of the total radial particle diﬀusive flux
in a real situation and, thus, the particle flux onto thedivertor plate can be diﬀerent aﬀecting the condition ofgen
the detached plasma. Therefore, Ωη > −1 can be an
important criterion in order to reproduce the behavior of
detached plasmas with the Braginskii model.

Acknowledgement(b)

Authors would like to thank Dr. D. Reiser of the
Forschungszentrum Jülich GmbH for fruitful discussion.

40
0
-40
0.0

gen

Ωη
0.2

0.4

0.6

0.8

= -1
1.0

s/L
(
)
Fig. 5 Parallel-to-B profiles of (a) 1 + Ωηgen δpi (solid line), theRHS of Eq. (10) (broken line), πgen(chain line) and (b)iΩηgen for the collisionless case.

in which the steady state (∂/∂t = 0) is assumed. Note thatgen
the second and third terms in Ωη are much smaller than
the first term and unity, respectively, and are negligible.
gen
Therefore, the essential diﬀerence of Ωη from Ωη is thatthe absolute value of the flow velocity gradient, |∂V/∂s|, is
not used. Note also that the first term in the RHS of Eq.
gen(10), which is denoted by πi hereafter, has almost thesame form as πi .
()
gen
In Fig. 4 (a), profiles are shown of 1 + Ωη δpi , thegen
RHS of Eq. (10) and πi estimated in the AIP modelfor the collisional case.( It shows) that the RHS of Eq.
gen(10) agrees well with 1 + Ωη δpi (i.e. Eq. (10) is
gensatisfied) and πi qualitatively approximates the behavior
(
)
gen
gen
of 1 + Ωη δpi . In this case, Ωη is positive or at leastgenclose to zero as shown in Fig. 4 (b) and, thus, Ωη ≈ Ωηholds. Therefore, πi in the B2 code approximates δpi wellin this case leading to a good agreement with the AIP
model.
Figure 5 shows the same kinds of profiles as Fig. 4 forgen
the collisionless case. It again shows that πi qualitatively()
genapproximates the behavior of 1 + Ωη δpi . In this case,

References[1] R. Schneider et al., Contrib. Plasma Phys. 46 (2006)
3.
[2] H. Kawashima et al., Plasma Fusion Res. 1 (2006)
031.
[3] K. Shimizu et al., Nucl. Fusion 49 (2009) 065028.
[4] T. D. Rognlien et al., J. Nucl. Mater. 196-198 (1992)
347.
[5] A. V. Chankin et al., J. Nucl. Mater. 390-391 (2009)319.
[6] K. Hoshino et al., J. Plasma Fusion Res. Ser. 9 (2010)
592.
[7] S. I. Braginskii, Reviews of Plasma Physics, vol. 1,
Consultants Bureau, New York, 1965, p. 205.
[8] W. Fundamenski, Plasma Phys. Control. Fusion 47(2005) R163.
[9] S. Togo et al., submitted to Contrib. Plasma Phys.
[10] S. Togo et al., 16th International Workshop on
Plasma Edge Theory in Fusion Devices (2017) O-05.
[11] S. Togo et al., Contrib. Plasma Phys. 56 (2016) 729.
[12] E. Zawaideh et al., Phys. Fluids 29 (1986) 463.
[13] S. Togo et al., J. Nucl. Mater. 463 (2015) 502.
[14] S. Togo et al., J. Comput. Phys. 310 (2016) 109.
[15] B. J. Braams, NET report no. 68 (EUR-FU/XII-80/87/68) (1987).
[16] R. H. Cohen and D. Ryutov, Contrib. Plasma Phys.
44 (2004) 111.

Simulation study of vapor shielding at a solid wall by aweighted PIC method
K. Ibano1∗ , A. Tanaka1, Y. Kikuchi2, S. Togo3, H. T. Lee1, Y. Ueda1, and T. Takizuka1
1

Graduate school of Engineering, Osaka University, Suita 565-0871 Japan
Graduate school of Engineering, University of Hyogo, Himeji 671-2280 Japan
3
Plasma Research Center, University of Tsukuba, Tsukuba, Japan
2

Introduction
When the plasma facing surface suffers intense plasma loads, vaporized wall material interacts with theincoming plasma flux and “shields” the remaining wall. This phenomena is known as vapor shielding and
now interests are being increased as an inherent mitigation [1–5]. In particular, transient heat loads in ITER
are expected to be considerably larger than ones in current devices [6], e.g., edge localized mode (ELM) canreach 0.1~1 GW/m2 for 0.1~1 ms. If the vapor shielding effectively mitigate the plasma energy, less amountsof heat loads are absorbed to the wall. In case of the vapor shielding at a solid wall, several physicalphenomena occur; melting of the surface, molten layer instability, droplet ejections, vapor emission and their
ionization. Thus, there are multiphase matters, i.e., plasma, gas, liquid, and solid phases and their interactions.
Effects of these reactions are examined previously by analytic modeling [7] and simulation works [8,9].
MHD fluid modeling were applied for these simulation works. In the fluid model, the boundary conditions
are usually set on a plasma-sheath entrance. However, the vapor shielding happens mainly in vicinity of asurface, and the reactions in the sheath region should be solved. Thus, we employ a Particle-In-Cell (PIC)
simulation. In the PIC simulation, kinetic behaviors of super particles are solved. The super particles arecomputation particle which represents numbers of real particles. Then, the PIC can simulate the particle
behavior in the sheath region. The simulation model is further described in the next section. In order tovalidate our simulation code, a comparison is taken with experimental observation from a plasma gun (PG);
observation of vapor shielding at a helium plasma irradiation on an Al coated W specimen. Heat flux
dissipation in the PG was observed for the Al coated W specimen due to the expanding Al clouds.

Model Description
In this study, a PIC code, named PIXY[10,11], is applied for the simulation. The PIXY code solves
trajectories of all particles (ions, electrons, and impurities) for 1 dimension in space and 3 dimensions invelocity (1d3v). Magnetic pre-sheath and Debye sheath can be simulated with this coordination. In the 1dsimulation space, central region is set as a Langevin heat bath region where electrons and ions are supplied.
The super particles move toward the boundary of the system. As they reached to the boundary, their energy
is counted as the heat flux to the wall. Then, a heat transfer calculation is taken for the wall component using
the heat flux as a surface boundary by solving a 1d heat transfer equation. The heat transfer calculation
determines the surface temperature, then evaporation rate of surface materials is determined. These
evaporated particles are fed into the plasma side. The interaction between the ejected vapor particles andincoming plasma are also solved. Several Atomic and Molecular reactions are treated by Monte-Carlotechniques; ionization, recombination, and radiations via line, Bremsstrahlung, and recombination.
Probabilities of these reactions in a time-step are calculated based on their reaction rates in the OPEN-ADAS
database [12] and compared with a random-number. When the reaction happens, corresponding energybalance is maintained by loss or gain of a randomly-selected electron in a same cell. As well as the A&M
reactions, ion-neutral collision is also simulated by a Monte-Carlo method. The ion-neutral elastic crosssections are calculated from the semi-classical approximation explained in [13] where the dipole
polarizability is taken from ref [14]. Once reaction is determined via a comparison between the reaction rate
and a random number for a neutral particle, an ion in a same cell is randomly selected as a binary collisionpartner. In order to avoid over counting of this collision effects, trajectories of neutral particles are calculated
in a three-dimensional space. Then, if neutral particles leave from a plasma tube, they are counted as pumped
out. In the Plasma Gun simulation, a 10 cm diameter plasma tube is used as the boundary.
---------------------∗ Corresponding author. E-mail: kibano@eei.eng.osaka-u.ac.jp, Phone: +81 6 6879 7234, Fax: +81 6 6879 7867

One issue is vapor ejection rate which changes significantly with surface temperature. If same superparticle weights (number of particles represented by a super particle) are used throughout the calculation,
numbers of super particles become significantly large. Thus, we use the weighted particle model. In thismodel, particle weights are different for all super particles. Then, the significant difference in emission rates
of impurity particles can be treated.

Results
In the plasma gun, vapor shielding by Al vapor cloud was observed in a helium plasma gun irradiation onto
an Al coated W specimen. The typical plasma parameters of the plasma gun are; electron temperature Te ~30 eV, density ne ~ 1021 m-3, and energy density 1 MJ/m2 for 0.2 ms. The detail of the PG experiment can befound in ref. [15,16]. The specimen was a thin W specimen with a thickness of 50 m coated by an Al thinlayer with a thickness of 3m. Specimen temperature was monitored by a pyrometer to the back surface.
Steady-state temperature of the W and Al/W specimens after the shots were 2250 K and 1600 K, respectively.
In the plasma gun simulation, the plasma with the 30 eV and 3 ×1021 m-3 is simulated. As Al ejection
starts, excitation, recombination and Bremsstrahlung radiation are observed. Then, as the PG pulse fully
develops and reaches to the specimen, massive evaporation of Al due to the ablation starts. Increasingradiation due to larger Al density dissipates the incoming plasma energy, then a temperature drop near
specimen is observed. Finally, the Al clouds expands to the entire calculation region and the PG plasmaenergy is effectively dissipated.
Simulated temporal change of the heat flux q” during the pulse is shown in Fig. 1. Developing heat flux
as the PG plasma reaches to the steady states is seen until 20 s after the pulse starts. Strong drop ofheat flux is seen at 100 s. In contrast, PG heat flux to the W specimen keeps constant during the PG
pulse. The heat flux dissipation is well seen from this result. A remarkable heat flux dissipation of ionsdue to ion-neutral collision, which is comparable to the radiation cooling without radiation trapping, isobserved. Thus, although the radiation transport and trapping models need to be included for the precise
heat flux estimation, it can be concluded that the ion-neutral collision is also important during the vapor
shielding as well as the radiation cooling on electrons.

Figure 1 Temporal change of heat flux during the plasma gun pulse for the plane W and the Al coated W
specimens. Individual heat flux for the Al coated W specimen with ion-neutral collision model is also
shown.[17]Summary
A weighted PIC code, PIXY, is applied for simulations of heat flux dissipation by vapor from a solid wall
(vapor shielding) observed at the plasma gun experiments using an Al coated W specimen. In the simulation,
1d3v PIC model including various A&M reactions combined with a heat transfer calculation at the wall isused. The weighted particle method is applied to treat the vapor emission flux largely varied by the wall
temperature. The plasma gun pulse was well simulated by the code and energy dissipation was reproduced.
---------------------∗ Corresponding author. E-mail: kibano@eei.eng.osaka-u.ac.jp, Phone: +81 6 6879 7234, Fax: +81 6 6879 7867

It was found that the ion-neutral collision as well as the radiation cooling is the important reactions for theenergy dissipation at the vapor shielding.
Acknowledgements. This work is performed on "Plasma Simulator" (FUJITSU FX100) of NIFS with the supportand under the auspices of the NIFS Collaboration Research program (NIFS16KNXN320). This work is partially
supported by the “ZE Research Program, IAE ZE29B-01”.

References
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]

J.K. Tripathi, T.J. Novakowski, G. Joseph, et al., J. Nucl. Mater. 464 (2015) 97–106.
T. Sizyuk, A. Hassanein, Nucl. Fusion. 54 (2014) 23004.
Y. Kikuchi, I. Sakuma, Y. Kitagawa, et al., J. Nucl. Mater. 463 (2015) 206–209.
G.G. van Eden, T.W. Morgan, D.U.B. Aussems, et al., Phys. Rev. Lett. 116 (2016) 135002.
G.G. van Eden, V. Kvon, M.C.M. van de Sanden, et al., Nat. Commun. 8 (2017) 192.
G. Federici, A. Loarte, G. Strohmayer, Plasma Phys. Control. Fusion. 45 (2003) 1523–1547.
D.I. Skovorodin, A.A. Pshenov, A.S. Arakcheev, et al., Phys. Plasmas. 23 (2016) 22501.
H. Wurz, I. Landman, B. Bazylev, et al., J. Nucl. Mater. 233 (1996) 798–802.
A. Hassanein, Fusion Eng. Des. 60 (2002) 527–546.
K. Ibano, D. Nishijima, J.H. Yu, et al., Nucl. Mater. Energy. 12 (2016) 278–282.
K. Ibano, S. Togo, T.L. Lang, et al., Contrib. Plasma Phys. 56 (2016) 705–710.
H.P. Summers, M.G. O’Mullane, Atomic data and modelling for fusion: The ADAS Project, in: AIP
Conf. Proc., 2011.
C.H. Liu, J.G. Wang, R.K. Janev, J. Phys. B At. Mol. Opt. Phys. 43 (2010) 144006.
P. Schwerdtfeger, Table Exp. Calc. Static Dipole Polariz. [online] (2015).
http://ctcp.massey.ac.nz/dipole-polarizabilities.
I. Sakuma, Y. Kikuchi, Y. Kitagawa, et al., J. Nucl. Mater. 463 (2015) 233–236.
I. Sakuma, Y. Kikuchi, W. Isono, et al., Plasma Fusion Res. 10 (2015) 1205089.
K. Ibano, Y. Kikuchi, A. Tanaka, et al., Contrib. to Plasma Phys. to be published---------------------∗ Corresponding author. E-mail: kibano@eei.eng.osaka-u.ac.jp, Phone: +81 6 6879 7234, Fax: +81 6 6879 7867

Simulation for stability of di-vacancy in tungsten
Kazuhito Ohsawa1, Ken Toyama2, Yuji Hatano3, Masatake Yamaguchi4, Hideo
Watanabe1
1Research Institute for Applied Mechanics, Kyushu University
2 Institute for Materials Research, Tohoku University
3Hydrogen Isotope Research Center, University of Toyama
4Center for Computational Science and e-system, Japan Atomic Energy Agency

1. Introduction
Divertor armor tiles installed in fusion reactors are exposed to intense plasma
particle irradiation. So, it is planning that the divertor armor tiles are covered bytungsten (W) or its alloys because they have excellent properties as plasma facing
materials (PFMs), e.g., low hydrogen (H) solubility, high melting point and high thermal
conductivity etc. However, a large amount of H isotope retention has been reported in
the W materials. In particular, tritium (T) inventory is an important subject associatedwith safety of fusion reactors. Usually, H isotopes are not retained in W crystal latticebecause of its low H solubility and large diffusivity. However, H isotopes areaccommodated in irradiated zone of the W materials. Vacancy-type lattice defects andimpurities are supposed to be plausible trap sites for the H isotopes [1].
In the present report, we pay attention to stable configuration of di-vacancy in W
lattice and other BCC transition metals for comparison. Besides, we investigate theeffect of H presence on the di-vacancy stability in those transition metals. Di-vacancies
in metals are usually stable. But formation energy of di-vacancy, i.e., binding energy of
two vacancies, in W is estimated to be negative [2]. So, it is expected that di-vacancy
nucleation does not occur in W. On the other hand, the presence of H or H isotopes will
enhance the growth of vacancy-type lattice defects.
Positron annihilation lifetime measurement is effective method to observe
vacancy-type lattice defects in metals. So, we calculate the positron lifetime to examinedi-vacancy formation in W.
2. Computational method
Binding energy of two vacancies for di-vacancy formation and binding energy ofvacancy-hydrogen complexes are defined in Eq. (1) and (2), respectively.

E B  2 EV  EV2 
EB  EVH m   EVH n   EV2 H mn  ,

(1)
(2)

where V and H represent vacancy and hydrogen atom, respectively. Therefore, V 2
indicates di-vacancy and VHm represents vacancy-hydrogen complex trapping mhydrogen atoms in a vacancy. Similarly, V2Hm+n indicates that di-vacancy trapping m+nhydrogen atoms. The function E is the cohesive energy of the system including suchlattice defects. Positive binding energy corresponds to attractive interaction between
the lattice defects. The cohesive energies are calculated in terms of first-principle
calculations on the basis of density functional theory. The calculations are performed byusing Vienna ab initio simulation package (VASP) with PBE type potential. The cut-offenergy of plane wave is 350eV. We used large simulation cell composed of 432 lattice
points (6x6x6 BCC lattice) in the present calculations in order to reduce the artificial
effects of periodic boundary condition imposed in the simulation cell. Lattice relaxations
are iterated until resultant force of each atom is lowered than 0.003 eV/Å.
Positron lifetime  is estimated in the overlap integral









 1  ca 0 2 n (r )n (r ) g (n n )dr(3)

V

where n  and n  are positron and electron charge densities, respectively. The
constants a 0 and c are classical radius of electron and light velocity, respectively. In thepresent calculation, the enhancement factor g (n n ) in Eq. (3) is estimated frompositron lifetime measured from bulk W experimentally. The electron charge density
around W atom is given from reference [3].
3. Results
Two types of configurations are possible for di-vacancy in BCC lattice, that is, thetwo vacancies line up in the <111> or <100> direction. We calculate binding energies of
two vacancies in a variety of BCC transition metals in Eq. (1). <111> arrangement of
two vacancies is stable in V (vanadium), Cr, Mo, and W, while <100> arrangement isfavorable in Fe, Nb and Ta. In addition to it, binding energies of vacancy-hydrogencomplex VH6 are also calculated in Eq. (2). Fig. 1 shows the binding energies of two
vacancies for di-vacancy formation and binding energies of VH6 for V2H12 formation in
the BCC metals. Fig. 2 shows binding energies of VHn in W in detail, where n=0, 1, 2, 3,
and 6 in the present work. These calculations indicate that di-vacancy is usually stable
than isolated two vacancies in BCC transition metals. However, binding energy of two
vacancies in W is negative. So, di-vacancy in W is only unstable. Besides, the presenceof H atoms contributes to stabilization of the di-vacancy.

Lifetime of positron trapped in di-vacancy is calculated in Eq. (3). Positron lifetime inbulk W lattice is experimentally measured as 118 ps [4]. According to it, theenhancement factor g (n n ) in Eq. (3) is estimated to be 1.91. So, lifetime of positron
trapped in di-vacancy introduced in W lattice is estimated to be 210 ps, which is good
agreement of experimental results.

Fig. 1: Binding energies of vacancies

Fig. 2: Binding energies of a variety ofand vacancy-hydrogen complexes invacancy-hydrogen complexes in W.

BCC metals.

4. Discussion
Vacancy cluster growth and the following void swelling is one of the troublesome
phenomena for metals in irradiation circumstance, which causes degradation of
structural materials. If formation energies of vacancy-type lattice defects are assumed
to be proportional to the surface area of the defects, vacancy cluster growth will occur so
as to reduce the formation energies. Actually, di-vacancy is more stable than isolated
two vacancies in many kinds of metals, as shown in Fig. 1. However, di-vacancy in W is
unstable exceptionally, according to the present calculations. On the other hand, thepresence of H and H isotopes contribute to binding of vacancies in W and other BCC
transition metals. Therefore, we infer that di-vacancy formation does not occur inannealing process in W specimen even if sufficiently high density vacancies are
introduced in the W specimen. On the other hand, the presence of H and H isotopes will
enhance the growth of vacancy-type defects, as shown in Fig. 3. Hydrogen plasma andirradiation circumstance coexist in fusion reactors. Therefore, the influence of H on suchvacancy-type lattice defects is one of the important subjects in the field of PFM.

Positron annihilation measurement in Fig. 4 is a sensitive experimental method toinvestigate the structure of vacancy-type lattice defects in metals. The size of
vacancy-type lattice defect is reflected to the positron lifetime. Usually, di-vacancy
formation in metals is proved by drastic extension of the positron lifetime at vacancy
diffusion temperature in annealing process. We hope that the behaviors of vacancies in
W specimen will be revealed by the positron annihilation measurement.

Fig. 3: Schematic view of bonding of
vacancy and vacancy-hydrogen complex in
W specimen. It is expected that di-vacancy
formation does not occur in W.

Fig.4: Schematic view of the positronannihilation lifetime measurement.

Acknowledgement
This work was supported by JSPS KAKENHI Grant Number 17K06993, and thecollaboration programs of Research Institute for Applied Mechanics, Kyushu University
and Hydrogen Isotope Research Center, Organization for Promotion Research,
University of Toyama.
References[1] K. ohsawa et al., Phys. Rev. B 85 (2010) 184117.
[2] L. Ventelon et al., J. Nucl. Mat. 458 (2015) 187.
[3] F. Herman et al., Atomic Structure Calculations, Prentice-Hall Inc. (1963).
[4] T. Toyama et al. J. Nucl. Mater. 499 (2018) 464.

Application of Genetic Algorithms to Modeling of Plasma
Physics
M. Honda
National Institutes for Quantum and Radiological Science and Technology, Naka,
Ibaraki 311-0193 Japan

Introduction
The genetic algorithm (GA) based optimizer is a powerful tool for global
optimization, such as finding a global maximum, not local maxima, which are typicallymore than a gradient-based method can manage. Sometimes it is advantageous comparedto other schemes with respect to the convenience, partly because a merit function is notnecessarily differentiable. Taking advantage of this feature, it is applied to physicalmodelings of magnetically confined toroidal plasmas. Here, we show some examples
such as searching the positions of the magnetic axis and fitting an analytic function to themeasurement points. Through the modelings, the wide applicability of the genetic
algorithm based optimizer is elucidated.

Genetic Algorithms
GAs are a class of metaheuristics inspired by the evolutionary biology. GAs makeuse of the ideas of the evolutional biology such as the natural selection, inheritance andvariability and turn them into the algorithm suited for computational calculations. GAs
are commonly used to robustly generate solutions to optimization and search problems
and, as such, they are known as GA-based optimizers. A GA-based optimizer has theadvantage in that it is applicable to global and multi-dimensional problems and does not
require a fitness function to be differentiable. By contrast, it is computationally slower
than other typical algorithms that can be applied to global problems and local
optimization methods. Capability of a global optimization is traded off for a slow
computation. In this paper, the public domain GA package PIKAIA [1, 2] is used as aGA-based optimizer.

Searching the magnetic axis
Searching a magnetic axis of an axisymmetricequilibrium in tokamaks computed by a GradShafranov solver is a sort of global optimization
problem. Equilibria we use are adjusted such that 𝜓𝜓

on the outermost nested closed flux surface becomeszero and the contour is concave upward, and thus themagnetic axis is the point where 𝜓𝜓 becomes

negative maximum. Hence, the magnetic axis may belocated at 𝒙𝒙 where the fitness function, 𝑓𝑓 𝒙𝒙 =−𝜓𝜓 for 𝒙𝒙 = (𝑅𝑅, 𝑍𝑍) , is maximized. Namely, anindividual has the information of the position in thedomain and the corresponding 𝜓𝜓 can be evaluatedFig. 1: Searching the magnetic axis.
The black crosses are the individuals
at the 20th generation and the red ×
mark is the fittest one.

as the fitness. An evolutionary run is carried out over 200 generations and the population
size is 100. The snapshots of the evolutionary runs at the 20th generation are shown in
Fig. 1 for the normal magnetic shear plasma. The majority of the individuals is confined
inside the plasma at this stage. Until around 30 generation, the fittest individual reaches
the magnetic axis.

Profile fitting
We now attempt to t a given function through the discrete measurement data in termsof the density and temperature proles. The double transport barrier formation has been
often achieved in JT-60U, which means that the edge and internal transport barriers are
simultaneously formed. In such cases, it is too difficult to use a polynomial function forprole fitting due to Runge's phenomenon, as discussed in Introduction. To fit such data,
the model would take the form:𝑔𝑔 𝜌𝜌; 𝒙𝒙 = 𝑥𝑥0 tanh[− 𝑥𝑥6 (𝜌𝜌 − 𝑥𝑥7 )] + 𝑥𝑥: tanh[− 𝑥𝑥; (𝜌𝜌 − 𝑥𝑥< )] + 𝑥𝑥=

where the radial coordinate 𝜌𝜌 is typically defined by the square root of the normalized

toroidal magnetic flux and 𝜌𝜌 ∈ [0,1]. The function 𝑔𝑔 is infinitely differentiable withrespect to 𝜌𝜌, a fact which makes it easy to use for numerical codes. The hyperbolic

tangent tanh 𝜉𝜉 is suitable for expressing the sharp change in the profile by adjusting theposition of 𝜉𝜉 = 0 to that of transport barriers. Hence, the model consists of two

hyperbolic tangent functions: one is for the ITB and the other is for the edge transportbarrier (ETB). The position of the transport barrier almost corresponds to 𝑥𝑥7 and 𝑥𝑥< ,respectively. The steepness of a barrier is controlled by 𝑥𝑥6 and 𝑥𝑥; . The factor 𝑥𝑥= is

imperative for fitting. As readily seen from the form of 𝑔𝑔, solving the simultaneous
equations after taking the derivative of 𝑔𝑔 with respect to 𝑥𝑥 may not be performed
analytically and is very difficult even numerically. A GA-based optimizer is a powerfultool to tackle this sort of problems. The fitness function may be given as 𝑓𝑓 𝒙𝒙 = 𝜒𝜒 D6

where 𝜒𝜒 6 =F

𝑦𝑦F − 𝑔𝑔(𝜌𝜌F ; 𝒙𝒙 )/𝜎𝜎F 6 . Here, 𝑖𝑖 denotes the index of a discrete set of themeasurement points, 𝑦𝑦F is a measured quantity and 𝜎𝜎F denotes the measurement error.

We now attempt to fit the measurement data of the electron and ion temperatureswith the function 𝑔𝑔. An evolutionary run has been executed until the 500th generation
and requires about 0.2 sec. It is obvious in Fig. 2 that thehyperbolic tangent functions play an important role inexpressing the steep temperature gradient in both thecore and edge of the plasma. The GA-based optimizer
has revealed its capability in that one single function isable to fit the measurement points in a wide variety ofdischarges. The numerical implementation for fitting israther simple and there is no need to compute theJacobian by differentiating the function unlike theGauss-Newton algorithm. Any function is of course
applicable and one may seek the form of a model
function, depending upon the situation.

References[1] P. Charbonneau, Astrophys. J. Suppl. Ser. 101 (1995)309.
[2] http://www.hao.ucar.edu/modeling/pikaia/pikaia.php

Fig. 2: Profile fitting of (a)the electron temperature and(b) the ion temperature forthe box type ITB plasma.

Pattern selection in drift wave turbulence with parallel shearflow
Y. Kosuga[1,2][1] Research Institute for Applied Mechanics, Kyushu University, Fukuoka, Japan[2] Research Center for Plasma Turbulence, Kyushu University, Fukuoka, Japan

Convective cells - nonlinearly generated ExB vorticities - play important roles in the
nonlinear dynamics of turbulent plasmas. As they are linearly unstable, convective cells
are nonlinearly driven from primary drift wave turbulence. There are two important
limits of convective cells. One is a radially elongated cell, which is called a streamer. The
other is a poloidally elongated cells, and is called a zonal flow. Since streamer and zonal
flow have different impact on turbulent transport, it is a critical issue to understand theirselection rule[1].
In this work, we discuss the role of parallel shear flow on the selection of zonal flow
and streamer in drift wave turbulence[2]. Employing a modulational instability model forsecondary flow generation, we show that parallel shear flow enters as a critical parameter.
Parallel shear flow modifies dispersion of underlying drift waves. This in turn leads to theenhancement/suppression of secondary flow response. As a result, shearing feedback on
the background fluctuation envelope becomes a function of parallel flow shear. Though
the modulational instability analysis, we show that there is a parameter window of parallel
flow shear for streamers to be selected in the system.
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 Grants-inAid for Scientific Research of JSPF of Japan (JP15H02155, JP15K17799, JP17H06089),
the Joint research project in RIAM.

References[1] S. Champeaux and P.H. Diamond, Phys. Lett. A 288 214 (2001)[2] Y. Kosuga, Phys. Plasmas 24 122305 (2017)

Enhancement and suppression of turbulence byenergetic particle driven geodesic acoustic modes
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

Problems on interactions between micro-turbulence and macro scale flows, such as zonal flows,
are important subjects in magnetically confined plasmas [1]. In toroidal plasmas, geodesic acoustic modes (GAMs), which are oscillatory zonal flows driven by turbulence and/or by energetic
particles, have attracted much attention, because GAMs are expected to suppress turbulence by
their velocity shear. Actually, the suppression of turbulence and transport have been reported inturbulence simulations [2]. However, recently, the enhancement of turbulence by GAMs driven
by energetic particles (which are called EGAMs) has been observed in gyrokinetic simulations,with the subsequent destruction of a transport barrier [3]. In this way, EGAMs can either mitigate or enhance the turbulence. This dual effect of the GAMs on turbulence requires theoreticalinvestigation. In this study, we investigate the phase-space dynamics of spatially inhomogeneous
turbulence in the presence of EGAMs. We show that EGAMs can carry clump of turbulence
spatially, which leads to the enhancement of turbulence in the regions where the turbulence isstable, while the turbulence is suppressed in the region where the turbulence is unstable [4].

2

Model

In this study, we consider the dynamics of spatially inhomogeneous turbulence with a transportbarrier in the presence of EGAMs. Details of the problem settings are written in [4]. We
consider the drift wave type turbulence, which is governed by the wave-kinetic equation, andfocus on the dynamics in the phase space, where the phase space consists of the real space, x(the radial direction in toroidal plasma) and the wavenumber kx . The spatio-temporal behaviorof the EGAM is given, and the time evolution of the turbulence in the phase space is numericallycalculated based on the wave-kinetic equation.

3

Results

The snapshot of the turbulence action, Nk , in the phase space is shown in Fig. 1. In the casewithout EGAM, the intensity of the turbulence is large in the region where the turbulence isunstable, and only the turbulence with small kx leaks into the stable region across the shear
layer (transport barrier). On the other hand, clumps of turbulence can be seen in the stable
1

L

γ
0

x

−1

0

0.01
0
2 (b)

k

γL

kx

0.01
0
2 (a)

−2

−2
1000

0
−1

0

−2

−2
1000

−3
0

(c)

−3
0

(d)

800

800

−1

t

t

−1

600

600
−2

400
0

20

x

40

−2

400
0

20

x

40

Fig. 1: Snapshots of Nk in the case of (a) VG = 0, and (b) VG = 0.1 are shown, where VG is theamplitude of the EGAM. The time evolution of turbulence intensity is shown for the case of (c)VG = 0, and (d) VG = 0.1 [4].
region, when the EGAM amplitude becomes finite. The clump of turbulence propagates withthe phase velocity of the EGAMs. This is due to the turbulence trapping [5] by the EGAM. We
found that the trapped turbulence wave-packets leak across the transport barrier. As a result,turbulence is enhanced by EGAMs in the stable region, while turbulence suppression is obtainedin the unstable region. The propagation of the turbulence is ballistic, with the phase velocity of
the EGAM. Hence, there appear a new global characteristic velocity for turbulence dynamics, inaddition to the local group velocity and that of the turbulence spreading [6]. This paper firstly
shows that the turbulence trapping in the phase space strongly affects the turbulence profile in
the real space, and that this effect can be a new mechanism of the turbulence propagation.

References[1] P. H. Diamond, et. al., Nucl. Fusion, 49, 045002 (2009).
[2] K. Hallatschek, D. Biskamp, Phys. Rev. Lett., 86, 1223 (2001).
[3] D. Zarzoso, et. al., Phys. Rev. Lett., 110, 125002 (2013).
[4] M. Sasaki, et. al., Scientific Reports, 7, 16767 (2017).
[5] P. Kaw, R. Singh, P. H. Diamond, Plasma Phys. Control. Fusion, 44, 51 (2002).
[6] O. D. Gurcan, et. al., Phys. Plasmas, 12, 032303 (2005).

2

Global mode analysis of ion-temperature-gradient instabilities using the gyro-fluid model inlinear devices
Tomotsugu OHNO1, Naohiro KASUYA2, Makoto SASAKI2, Masatoshi YAGI3
1. Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, 6-1 Kasuga-Koen, Kasuga, Fukuoka
816-8580, Japan
2. Research Institute for Applied Mechanics, Kyushu University, 6-1 Kasuga-Koen, Kasuga, Fukuoka 816-8580, Japan
3. National Institutes for Quantum and Radiological Science and Technology, 2-166 Omotedate, Obuchi,
Rokkasho-mura, Aomori 039-3212, Japan

1.

Introduction

It is important to understand turbulent phenomena in magnetically confined plasmas. A level of turbulent transport isdetermined by competition of several instabilities and formation of turbulent structures [1]. Fundamental mechanisms ofstructural formation in turbulent plasmas can be studied in simple cylindrical configurations [2]. One of the candidates tocause anomalous transport in fusion plasmas is the microscopic instability due to the ion temperature gradient (ITG) [3].
Studies of excitation conditions of the ITG instability have been carried out in basic experiments [4]. The ion temperature
can be measured using diagnostics as an ion-sensitive probe [5]. On linear device PANTA [6], ion temperature
measurements and numerical analyses of the ITG instability have been progressing. The purpose of our research is topredict excitation conditions of the ITG instability in PANTA. The ITG instability is predicted to be excited when the ratioof the ion temperature gradient to the density gradient exceeds the threshold value near the unity [7]. Numericalsimulations by using the fluid model in PANTA shows that a mode with k⊥ρs ~ 1 can become unstable even in low ion
temperature plasmas as in linear devices. Here, k⊥ and ρs are the wavenumber in the perpendicular direction to themagnetic field and the effective Larmor radius, respectively. Therefore, a detailed analysis including the finite Larmor
radius (FLR) effect is required. Analyses using a set of localized gyro-fluid equations have been carried out [8]. The ion
Larmor radius is not much smaller than the plasma radius in the linear device, so the global mode structure must be
calculated. In this research, a numerical simulation code is developed to solve global eigen-functions of the ITG instability
with the gyro-fluid model. Here we report mode structures and parameter dependences of the ITG instability in PANTA byusing the global code.

2. Global Mode Analysis
2.1 Gyro-Fluid Model
A set of gyro-fluid equations is derived by taking the moments of the electrostatic gyro-kinetic equation in the velocity
spaces [9]. The target plasma has a cylindrical configuration with a homogeneous magnetic field parallel to the axial
direction, so the magnetic curvature terms can be eliminated. Applying the gyro-kinetic ordering, the linear forms of theequations in the cylinder are given to be⎛ˆ 2 ⎞ ρ ∂Ψ∇dn
+ ∇ // u// + ⎜⎜1+ η⊥ ⊥ ⎟⎟ s
= 0,
dt
2 ⎠ Ln ∂y⎝

1 dT//ρ ∂Ψ
2ν+ ∇ // ( 2u// + q// ) + η// s= − ii (T// − T⊥ ),
τ dt
Ln ∂y
3τ

(1)

du//
+ ∇ // τ n + T// + Ψ = 0,
dt

(2)

(3)

⎡∇
⎤
ˆ2
1 dT⊥
ˆˆ 2 ⎥ ρs ∂Ψ = ν ii (T − T ),+ ∇ // q⊥ + ⎢ ⊥ + η⊥ 1+ ∇
//
⊥
τ dt
2
⎣
⎦ Ln ∂y 3τ

(4)

(

)

(

)

where n is the ion density, u// is the ion velocity, T is the ion temperature, τ = Ti0 / Te0, Ti0 and Te0 are ion and electron
temperatures at the plasma center, η = Ln / LT, Ln is the density gradient length, LT is the ion temperature gradient length, q
is the heat flux, and νii is the collision frequency between the ions. The subscripts // and ⊥ represent the quantities in theparallel and perpendicular direction to the magnetic field, respectively. The time and spatial length are normalized with Ωci−1and ρs, where Ωci is the ion cyclotron frequency. Ψ is the gyro-averaged potential Ψ ≡ Γ01/2Φ, where Γ1/2
0 = (1− bτ / 2 )and b = −∇ ⊥2 . Quantity b represents the magnitude of the FLR effect. The FLR effect is included in Ψ, ∇ˆ ⊥2 and ∇ˆˆ 2 terms,
where
∇ˆ ⊥2
bτ / 2Ψ≡−Ψ = Ψ1,2(1+ bτ / 2)

The quasi-neutrality relation is given to be

(5)

ˆˆ 2
∇
Ψ ≡ 2Ψ1 +b 2τ 2 / 2(1+ bτ / 2)

2

Ψ.

(6)

⎛
⎞
b2
Ψ
Γ0 ⎜ n −T⊥ ⎟ − (1− Γ0 ) = Ψ
1+
b
τ
2
τ
⎝
⎠

(7)

to determine the relation between the density and potential. Collisions are dominant in this system, and higher ordermoments of the gyro-kinetic equation give simplified forms of the heat flux as follows:3
1q// = −
∇ //T// ,q⊥ = −∇ //T⊥.
(8)
(9)
ν iiτ
ν iiτ

2.2 Numerical Scheme
To obtain the linear eigen-mode and -frequency, a spectral code is developed. The Bessel expansion in the r direction andFourier expansions in the θ, z directions are applied on Eqs. (1) - (4) and (7) by using
M
J L
⎛π
⎞
φ ( r,θ , z ) = ∑ ∑∑φmjl J m λmj r / a e imθ cos ⎜ ( 2l −1) z / Lz ⎟,⎝2⎠
m=−M j=1 l=1

(

)

(10)

M
J L
⎛π
⎞
u// r,θ , z = ∑ ∑∑ u//mjl J m λmj r / a e imθ sin ⎜ 2l −1 z / Lz ⎟.
⎝2⎠
m=−M j=1 l=1

(

)

(

)

(

)

(11)

Density n and temperature T are also expanded in the same way as in Eq. (11). The equations are solved in the region
between r = 0 (center of the plasma) and r = a (outer boundary of the plasma). The boundary condition in the radial
direction are set to f = 0 at r = 0, a when m ≠ 0, and ∂f / ∂r = 0 at r = 0, f = 0 at r = a when m = 0. Functions Jm(λmj r / a) areused for the expansion to satisfy the boundary conditions, where Jm(r) is the Bessel function, λmj is the j-th point withJm(λmj) = 0, and m is the azimuthal mode-number. Periodic boundary conditions are used in the axial direction, and theeffect of the end plate is not considered. We solve the matrix of the equations with the spectral expansion using the
eigenvalue method. In the eigenvalue method, the eigenvalues of the matrix are obtained using the mathematic library
MKL to give the linear growth rate.
For the linearization, the differential operator d / dt and ∇// are replaced to iω and ikz, where kz is the wavenumberin the parallel direction. The real and imaginary part of ω give the frequency and growthrate, respectively. The axial mode
number is set to l, which gives kz = 2π l ρs / Lz, where Lz is the device length. In the local analysis in Ref. 13, the radial andazimuthal wavenumbers are assumed to be same to give k⊥2 = 2kθ2 = 2(mρs / a)2, but in the global model k⊥ are given to be
λ̂mj = λmj ρ s / a for function Jm(λmj r / a).

2.3 Target Plasma
The geometry of the plasma is a simple cylindrical shape,
and the magnetic field has only the axial component with a uniform intensity. For the simulations, experimental
parameters in PANTA are used; device length Lz = 4.0 m, plasma radius a = 0.07 m, density n0 = 1.0×1019 m-3, Ln = 0.07 m,
νii = 350 s-1, magnetic field B = 0.1 T, temperatures Te0 = 3 eV and Ti0 = 0.3 eV. The temperatures and magnetic field giveρs = 1.1 cm, ρi = 3.5 mm and Ωci/2π = 3.8×104 Hz for argon plasmas. With these parameters, ρi/Ln = 0.05 ≪ 1.0, kzρi =5.5×10−3 ≪ 1.0 with l = 1 and k⊥ρi = 0.02 ≪ 1.0 with m = 1, so the gyro-kinetic ordering is satisfied. The otherparameters for the analysis are τ and ηi.

3.

Linear Instability

Numerical analyses are carried out using the developed gyro-fluid code. Figure 1 shows the dependences of thegrowthrate on the azimuthal and axial mode number. Here simplification with η// = η⊥ ( = η ) is applied. The mode
numbers of the most unstable modes in this parameter are m = 2 and l = 1. The critical value ηc of the instability in thiscondition is about 1.2, and the value of the growthrate increases as η increases. The global code can give the radial
distribution of T⊥, T//, n, u//, Φ, and Ψ in addition to the growthrate and eigen-frequency. Figure 2 shows the radial profilesof T⊥, T//, n and u// with m = 2, l = 1, η = 1.2 and τ = 1 of an argon plasma. The growthrate and frequency are given to be ωi= 2.9×10-3 and ωr = −9.6×10-3. Similar results were obtained by the analysis with the local model introducing the kr value
as a parameter.
The ion temperature is the most important parameter for the ITG instability, so the dependency on the magnitude andgradient length of the ion temperature is evaluated. Figure 3(a) shows the contour plot of the growthrate in the τ and η
space. The critical value ηc for the ITG instability changes depending on the magnitude of τ. The minimum of ηc is ηc = 0.8,
and increases as τ increases (ηc = 1.7, when τ = 2.0, for example). These are the same tendencies as in the local model
analyses as in Fig. 3(b). The cross in Fig. 3 indicates one of the experimental conditions in PANTA (τ = 0.1, η = 0.2 withan argon plasma). It suggests that a higher temperature gradient is needed to observe the excitation of the ITG mode in
PANTA.

Fig. 1: Dependency of growthrate on m and l with τ = 1, an argon
plasma.

Fig. 2: Radial profiles of eigenmodes of (a) T⊥, (b) T//, (c) nand (d) u// with m = 2, l = 1, η = 1.2, τ = 1, an argon
plasma.

4.

Fig. 3: Contour plot of the growthrate in the τ and η space with
m = 2, l = 1, an argon plasma. The cases with the (a) global and(b) local analysis are plotted.

Parameter Dependency of Radial Structures

Dependencies on the other experimental parameters are also evaluated. The ion mass depends on the using discharge gas,
which makes the Larmor radius different. The effective plasma radius changes, because the spatial length is normalized by
the Larmor radius. Figure 4 shows the ion mass dependences of critical ηc with m = 1-5 in the cases with the global andlocal analyses. The cases of helium, neon, and argon discharges are plotted. The difference in the dependency on the mode
number is partly due to the selection of kr. In the global model, kr of the most unstable mode is determined by calculating
the eigenmode, though it is only a parameter in the local model. In the case with helium, as shown in Fig. 4, the mdependence is greatly different. The local analysis gives critical value ηc = 1.9, which is larger than that ηc = 1.2 in the casewith the global analysis, when m = 1. This is because the value of k⊥ of the most unstable mode is larger than the value
obtained from kr = kθ assumed in the local analysis. Figure 5 shows the relative amplitudes of the k⊥ components in the
eigenmodes of the density. The maximum value is normalized to be the unity. λ̂mk represents the normalized wave length of
the base Bessel functions. As is seen in Fig. 5, there is a peak at λ̂mk = 0.7, but in the local analysis the smaller value k⊥ =2 (2m / a) ρs = 0.14 is used. Figure 5 also shows that the value of λ̂mk giving the maximum amplitude is almost the samein spite of the ion mass variation.
In addition, the difference in the treatment of kθ is the other cause. The value of kθ is given as m / r, and is evaluated at
one point, typically at r = a / 2 in the local model. On the other hand, in the global model, since 1 / r terms are expanded,
so the following relation is used for the azimuthal derivative;
K
1
∂
2mAk
Ai = ikθ A i = i∑ 2
dr* J m (λm,k r* )J m (λm,i r* ).
∫0∂y
k=1 J m+1 ( λm,i )(12)This calculation gives connection between the various Bessel functions. The radial structure also affects the evaluation ofkθ. Figure 6 shows the radial profiles of the amplitudes of n when m = 2, l = 1, η = 1.2 and τ = 1. In the case of helium, thestructure localized near the center of the plasma is formed as compared with the cases of argon and neon. This is because
the characteristic spatial scale ρs is changed. From this result, it is not appropriate under this condition to set the evaluation
position at r = a / 2 in the local analysis. Since the differences in the mode structures greatly affect the growth rates of the
ITG instability in this way, it can be said that it is meaningful to analyze the detailed mode structures by the global
analysis.

Fig. 5: Relative amplitudes of the k⊥ components in the eigenmodesof the density n, when τ = 1, η = 1.2, m = 1 and l = 1.
Dependences on, normalized wave length, are plotted.

Fig. 4: Ion mass dependences of the critical value ηc, when τ= 1. Those of modes with m = 1 - 5 are shown in thecase with the (a) global and (b) local analyses.

Fig. 6: Radial profiles of amplitudes of the density n with m = 2, l =1, η = 1.2, τ = 1.

5. Summary and Discussion
To investigate the excitation conditions of ITG modes in PANTA linear device, the linear instability analyses were carriedout. A gyro-fluid code, which solves global eigen-functions including the FLR effect, was developed and used for the
analysis. The same tendency from the local model was obtained, which gives the threshold for the instability in the linear
device. Wavenumber k⊥ of the most unstable mode was obtained with the radial mode structure analysis. The radial
structure is affected by the relative size between the plasma radius (to give the boundary condition) and the Larmor radius(to give the typical spatial length of the instability). The mode is more localized near the center of the plasmas in the caseof a discharge gas with a smaller mass number as helium.
From the gyro-fluid model, it was found that the experimental parameters of PANTA are present in the stable region,and for excitation of the ITG instability 4 times or more larger η is required. In linear machine CLM, the ITG instability
has been observed in the case when the radial positions of the steep temperature and density gradient are different fromeach other to satisfy the excitation condition on η [10]. The condition of destabilization should be explored by the
distribution controls of the temperature and density in the future.

Acknowledgements
This work is supported by JSPS KAKENHI Grant Number JP16K06938, and the collaboration program of NIFS
(NIFS17KNST112, NIFS16KNXN323) and RIAM of Kyushu University.

Reference[1] P.H. Diamond, et al., Plasma Phys. Control. Fusion 47, R35 (2005).
[2] G. R. Tynan, A. Fujisawa and G. McKee, Plasma Phys. Control. Fusion 51, 113001 (2009).
[3] W. Horton, Rev. Mod. Phys. 71, 735 (1999).
[4] A.K. Sen, et al., Phys. Rev. Lett. 66, 429 (1991).
[5] N. Ezumi, et al., J. Nucl. Mater. 337, 1106 (2005).
[6] S. Inagaki, et al., Sci. Reps. 6, 22189 (2016).
[7] Y. Miwa, et al., Plasma Fusion Res. 8, 2403133 (2013).
[8] G. Hattori, et al., Plasma Fusion Res. 10, 3401060 (2015).
[9] W. Dorland and W. Hammet, Phys. Fluids B 5, 812 (1993).

磁場閉じ込めプラズマにおける主要素以外の効果に関する検討:
不純物と中性粒子
On Effects of Impurities and Neutrals in Magnetically Confined Plasmas
N. Kasuya M. Yagi 1, A. Fukuyama 2, T. Ohno 3, M. Sasaki, S. Abe 3, T. Kobayashi 4 and
S. Inagaki
Research Institute for Applied Mechanics, Kyushu University, Kasuga, Fukuoka 816-8580,
Japan
1
National Institute for Quantum and Radiological Science and Technology, Obuchi,
Rokkasho, Aomori 039-3212, Japan2
Department of Nuclear Engineering, Kyoto University, Kyoto 606-8501, Japan
3
Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, Kasuga,
Fukuoka 816-8580, Japan
4
National Institute for Fusion Science, Toki, Gifu 509-5292, Japan
E-mail: kasuya@riam.kyushu-u.ac.jp
Magnetically confined fusion plasmas consist of several elements. Not only main plasma but
also their isotopes, impurities and even neutral particles are important players for the plasma
confinement anywhere in the core, edge and SOL regions. These synergetic effects must be
clarified. In this research, two kinds of topics are considered regarding this object: i)development of a simulation tool for investigating the impurity control in fusion reactor byusing the integrated simulation code TASK, and ii) turbulence simulation taking account of
inhomogeneity of neutral density in linear devices by using the resistive drift-wave turbulencecode NLD.
The first topic is on the development of a simulation tool for investigating the impurity
control. High Z impurities comes from the vessel wall as tungsten used in ITER, soestablishment of the control method is crucial for the plasma confinement. There are someexperimental results that electron cyclotron resonant heating (ECH) is effective to spill out theaccumulated impurities in the core region [1]. To clarify the appropriate condition for the
impurity control, we are developing an integrated simulation code to evaluate transport ofmulti species ions in tokamaks. The transport code TASK/TR [2] is being modified for that
purpose. To calculate transport of tungsten, atomic processes as ionization and recombination
between multiple charge states must be considered. Near the center of plasma, where theplasma pressure gradient is rather small, the neoclassical transport process is the key todetermine the impurity transport [3]. For the first step, response of the background plasmas iscalculated by applying an additional ECH, whose power deposition position is localized near
the plasma center. A monotonic response is obtained in accordance with the input power in theparameter region of the ASDEX experiment, which is possible to affect the impurity transport.
In the present version of TASK/TR, several kinds of transport models are installed. For theneoclassical transport, Hinton-Hazeltine model [4], Chang-Hinton model [5], and NCLASS
model [6] can be selected. Three cases using these models are compared, and show no
qualitative difference.
The second topic is on the effect of inhomogeneity of neutral density in linear devices. It is
important to consider combination of inhomogeneities, which drive and/or damp instabilities

in magnetized plasmas. Inhomogeneity of neutral particle density is taken into account foranalyzing turbulent states in a simple cylindrical configuration by combining two kinds ofsimulations in this research. One is a simulation for 2-D profiles of neutral particles, and theother is that for resistive drift wave turbulence. Radial and axial profiles of neutral particles
are calculated by the Monte Carlo method [7] with a configuration of linear device PANTA[8]. The ratio for ionization of neutral particles depends on the electron temperature, and thehigher temperature in the operation range makes the neutral density smaller near the center ofthe plasma. Figure 1 shows the radial profile of the neutral density. Then the radial profile ofthe ion-neutral collision frequency in is set by introducing these neutral density profiles into
resistive drift wave turbulence simulations [9]. Smaller ion-neutral collision frequency makes
the instabilities more unstable, and its inhomogeneity enhances mean electrostatic potential
formation [10]. Figure 2 shows the dependency of the potential on the dip of the in profile,where the values of in at the center of the plasma is related to (1 - d1). The potential at r = 0
shows rapid increase when d1 < 0.5 as the in dip is increased. Plasma flows can be driven by
turbulent stresses from excited modes in magnetized plasmas, and convective derivative
nonlinearity generates the azimuthal flow near the center of the plasma in this case. The radialinhomogeneity of the neutral density affects azimuthal flow generation by changing the phasestructure of the most unstable eigenmodes [11]. Eigenmode analyses show that the mode
structure has a complex Bessel-type function shape in the central region of the plasma, andthe imaginary part arises from the radial inhomogeneity of the damping term caused byion-neutral collisions. The amplitude of turbulent stress is proportional to the inhomogeneity
under a marginally stable condition. Two roles of the neutral particles on the turbulence, directchange of linear instabilities and formation of background flows, become clear. Global
structural formation is also an important factor for determining the plasma turbulent state, andthis result clearly shows that several kinds of radial background distributions, the plasma andneutral densities in this case, can influence the global structures.

FIG. 1: Radial profiles of the neutral density,
averaged in the axial direction. The caseswith Te = 3, 5 [eV] and M = 0.1, 0.3, 0.5 areshown.

FIG 2. Dependency on the in dip of theformed potential at the plasma center (r = 0).
d1 represents the value of in at r = 0
normalized with max(in), and the spatiallyconstant in profile is given with d1 = 1.

Acknowledgements
This work was partly performed under the Joint Special Design Team Collaborative Research
Program (29K019), and was partly supported by JSPS KAKENHI Grant Number
JP16K06938, by the collaboration program of NIFS (NIFS17KNST112, NIFS16KNXN323)

and RIAM of Kyushu University. Some numerical simulations were carried out on "Plasma
Simulator" (FUJITSU FX100) of NIFS.
References[1] NEU, R., et al., Plasma Phys. Control. Fusion 44 (2002) 811.
[2] FUKUYAMA, A., http://bpsi.nucleng.kyoto-u.ac.jp/task/.
[3] ANGIONI, C. et al., Nucl. Fusion 54 (2014) 083028.
[4] HINTON, F. L. and HAZELTINE, R. D., Rev. Mod. Phys. 48 (1976) 239.
[5] CHANG, C. S. and HINTON, F. L., Phys. Fluids 29 (1986) 3314.
[6] HOULBERG, W. A., et al., Phys Plasmas 4 (1997) 3230.
[7] IGNATENKO, M., et al., Jpn. J. Appl. Phys. 46 (2007) 1680.
[8] INAGAKI, S., et al., Sci. Rep. 6 (2016) 22189.
[9] KASUYA, N., et al., Phys. Plasmas 15 (2008) 052302.
[10] KASUYA, N., et al., Phys. Plasmas 25 (2018) 012314.
[11] KASUYA, N., et al., J. Phys. Soc. Jpn. 87 (2018) 024501.

13th Burning Plasma Simulation Initiative (BPSI) Meeting

Simulation studies on transport of multi-species
plasmas including impurity ion
Masanori Nunami
National Institute for Fusion Science, Toki 509-5292, Japan
Anomalous and neoclassical transport of multi-species plasmas including impurity ion areinvestigated by gyrokinetic and drift-kinetic simulations. In a high ion temperature plasmas
heated by neutral beam injection in Large Helical Device, the gyrokinetic simulations indicatethat the particle ﬂuxes can be varied according to changes of plasma proﬁles and the constitution of particle species, and the turbulent contribution of the impurity particle ﬂux remains
radially inward directed within certain changes of the temperature gradients. On the other
hand, the drift-kinetic simulations show that the neoclassical particle ﬂuxes of ions can change
due to generated radial electric ﬁeld and the external momentum sources.

1 Introduction
In order to design the fusion reactors, the transport phenomena of plasma heat and particles
should be accurately predicted by theoretical approaches, and numerical simulation approaches
based on the kinetic frameworks are powerful for it. Simulation studies on the transport ofthe plasma consists of multi-ion-species are strongly needed to understand the physics of theburning plasma in the ITER, future reactors, and also the plasma in stellarators such as theLarge Helical Device (LHD) [1], because there exists quite complicated transport properties
through the multi-transport channels produced in the multi-species plasmas. In addition, in themulti-ion-species LHD plasmas heated by neutral beam injection (NBI), the extremely hollow
impurity density proﬁles called impurity hole are often observed, and understandings of theirgeneration mechanism is a critical issue for magnetically conﬁned plasma researches. In suchplasmas, the transport ﬂuxes of heat and particles should be calculated by kinetic simulations
for the neoclassical and anomalous contributions to the transport, where the transport should bedetermined by the total transport ﬂuxes including the both contributions. In this work, we try
to calculate the transport ﬂuxes of multi-ion-species LHD plasmas with impurity hole structure
by gyrokinetic and drift-kinetic simulations.

author’s e-mail: nunami.masanori@nifs.ac.jp

13th Burning Plasma Simulation Initiative (BPSI) Meeting

2 Turbulent transport
Turbulent transport has been considered to be one of the most critical issues in the magneticallyconﬁned plasmas researches. In recent our linear gyrokinetic analyses of the micro-instability inthe LHD impurity hole plasma, which consists of four species (e, H, He, and C), it was foundthat the ion temperature gradient (ITG) mode is a dominant instability [2]. Here, we performthe nonlinear gyrokinetic turbulence simulations by the GKV code [3, 4] in order to evaluate
turbulent contributions of the plasma transport. In Fig. 1, ITG-dependences of the turbulent
contributions of the radial transport ﬂuxes of the heat and particle for all species satisfying theambipolar conditions. All heat ﬂuxes are increased with increasing ITG, which is one of thetypical characters of the ITG-driven turbulent transport. On the other hand, we ﬁnd that theparticle transport ﬂuxes have quite diﬀerent ITG-dependencies among diﬀerent particle species.
Furthermore, the impurity carbon particle ﬂux remains the negative (radially inward directed)
within wide ranges of the ion temperature gradients from their nominal values.

(trb)

(trb)

(right) of each species(left), and particle ﬂuxes ΓsFig. 1 GKV results for turbulent heat ﬂuxes Qs
as functions of the normalized ion temperature gradients (LT i ≡ (−d ln Ti /dr)−1 ) in the LHD
impurity hole plasma. Here LT i0 are nominal temperature gradient lengths of ion, and each ﬂuxes
are normalized by the arbitrary unit.

3 Neoclassical contributions
The neoclassical contributions of the transport of the LHD plasma are evaluated by using the
DKES/PENTA code [5]. The code solves the drift-kinetic equation and evaluate the ambipolar
radial electric ﬁeld Er . Since the neoclassical particle ﬂuxes are quite sensitive to Er and theexternal momentum torque by NBI heating. We evaluate the eﬀects of the external torque13th Burning Plasma Simulation Initiative (BPSI) Meetingon the neoclassical ﬂuxes, and it is found that the co-injected external torque can cause not
only ion-root with negative Er but also electron-root with positive Er [6]. In the case of theelectron-root, the neoclassical particle ﬂux of carbon can be outward directed which is opposite
directed to the turbulent ﬂuxes. On the other hand, the ion-root Er is often observed in thehigh-Ti LHD plasma [7], it should be discussed the possibilities of outward-directed neoclassicalimpurity carbon particle ﬂux in the case of the ion-root. In the simulations performed here,
we have not yet included other important eﬀects, e.g., the poloidal variation of the electrostatic
potential, and the direct contributions of the radial beam ﬂux Γb to the ambi-polar condition.
The studies on the neoclassical impurity particle ﬂux with above eﬀects will be performed as
future works.

Acknowledgment
This work is supported in part by the RIAM collaborative Research Program (29FP-21), theJapanese Ministry of Education, Culture, Sports, Science, and Technology, Grant (26820398,
16K06941, 17K14899), NIFS collaborative Research Programs (KNST095, KNTT043, KNTT035),
and the FLAGSHIP2020, MEXT within the priority study 6. The simulations are performed on”Plasma Simulator” (FUJITSU FX100) of NIFS.

[1] Y. Takeiri, et al., Nucl. Fusion 57, 102023 (2017).
[2] M. Nunami, et al., 26th IAEA FEC, TH/P2-3 (2016).
[3] T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).
[4] M. Nunami, T.-H. Watanabe, and H. Sugama, Plasma Fusion Res. 5, 016 (2010).
[5] D.A. Spong, Phys. Plasmas 12, 056114 (2005).
[6] M. Nunami, et al., Plasma Fusion Res. 12, 1203039 (2017).
[7] T. Ido, et al., Plasma Phys. Controlled Fusion 52, 124025 (2010).

Simulation study on internal transport barrier formation using gyrofluid model(II)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 ITBTOK 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 recent simulation results on ITB formation by the improved transport model.
2. Model equations
GF-ITG code consists of three field equations (ion continuity equation, parallel momentum
equation and ion temperature evolution equation)[1,2]. Here, we have improved ion
temperature evolution equation as
⎛ dT
⎛ dn
3
1 ∂Φ ⎞
1 ∂Φ ⎞
n0 (r) ⎜+ κT⎟ − Ti (r) ⎜ + κ n
⎟=
⎝ dt
⎝ dt
2
r ∂ϑ ⎠r ∂ϑ ⎠592εω̂ d T −A ∇ // T + A∇ //V + χ ⊥ ∇ ⊥2 T + ST
2τ
5
5 π

(1)

where the quasi-linear modification of ion temperature is taken into account, namely
Ti (r) = T0 (r) + T!0,0 (r, t).

The Gaussian type of heat source is introduced as,

ST (r) = P exp(−(r − Sr )2 / Δ 2 )(2)

The total power is estimated by⎛⎛ Sr2 ⎞⎛ (1− Sr )2 ⎞⎞ VPSr Δ (1−Sr )/Δ
PT = 4π Ra ∫ r drST (r) = VPΔ ⎜ exp ⎜ − 2 ⎟ − exp ⎜ −exp(−ξ 2 )d ξ
⎟⎟ +
∫
2
Δ
2
⎝ Δ ⎠
⎝
⎠⎠
⎝
0
−Sr /Δ
1

2

2

2

(3)

where V=2p2Ra2 is the volume. For Sr=0, D=0.2, PT=0.04VP and for Sr=0.4, D=0.05,
PT=0.0177VP. In the following simulation, P=10, Sr=0.4, D=0.05 are adopted.
3. Simulation results
Figure 1 shows mode spectrum in Fourier space (Type I). Here, resonant modes are located in
the area where the relation qa < n / m < qmin is satisfied. n is the toroidal mode number and mis the poloidal mode number. The safety factor profile q is shown in figure2, where qa=3,
qmin=1.35. The ion temperature profile and ion density profile are also shown. In Type I mode
selection, low n off-resonant modes are also taken into account.

Fig.1 Mode spectrum in Fourier space (Type I).

Fig.2 Safety factor profile, ion temperature
profile and ion density profile.

Figure 3 shows temporal evolution of ion temperature profile.

It is seen that clear ITB is notformed in this model with quasi-linear profile modification. The contour plot of fluctuating
ion temperature at t=20 is shown figure 4. The strong convective transport is induced by the
local heating.
Next, we investigate the effect of off-resonant modes on ITB formation. Figure 5 shows mode
spectrum in Fourier space (Type II). The red dots show the resonant modes and low n
off-resonant modes and blue dots show the medium and high n off-resonant modes. The ion
temperature profiles for cases with Type I and Type II are shown in figure 6. It is found thatthese off-resonant modes weakly affect ion temperature profile evolution.

Fig.3 Temporal evolution of ion temperature

Fig.4 Contour plot of fluctuating ion

profile.

temperature at t=20.

Fig.5 Mode spectrum in Fourier space (Type

Fig.6

II).

profiles for Type I and Type II mode

Comparison

of

ion

temperature

selections.

4. Discussion
We have performed the flux driven simulation using gyro-fluid code taking account of thequasi-linear modification of ion temperature on ITB formation and investigate the time
evolution of ion temperature profile using mode spectrum pattern Type I. It is found that thiseffect inhibits ITB formation. The effect of off-resonant modes on ITB formation is also
investigated changing mode spectrum pattern. It is shown that Type II off-resonant modes
very weakly affect ITB formation. Since quasi-linear modification of ion temperature profilestrongly affects ITB formation, so that we should carefully reconstruct ion transport equation
as well as Shafranov shift effect. It is left as a future work.
Acknowledgments
The author would like to thank Drs. N. Miyato, A. Matsuyama and T. Takizuka for usefulcomments. The computation was carried out using the SELENE supercomputer system at
IFERC-CSC.
References[1] S. Tokunaga et al., Nucl. Fusion 49 075023(2009).
[2] M. Yagi et al., 14th BPSI meeting, RIAM, Kyushu Univ. (2016)
https://www.riam.kyushu-u.ac.jp/sosei/bpsi/2016/BPSI2016_proc.pdf