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

第 16 回核燃焼プラズマ統合コード研究会
16th Burning Plasma Simulation Initiative (BPSI) Meeting(Ver.1.1)
日時：2018 年 11 月 29 日(木)－30 日(金)場所：九州大学筑紫キャンパス 応用力学研究所 2 階大会議室(18 min talk+7 min discuss or (15 or 10) min talk+5 min discuss)
11 月 29 日(木)9:10 – 12:30 核融合エネルギーフォーラムサブクラスターとの合同会合
12:30 – 13:30 昼休み
13:30 – 13:40 はじめに
村上 (京大)(座長：村上)13:40 – 14:05 講演 1-1
林 (量研)
Predictive integrated modelling of plasmas and their operation scenarios
towards exploitation of JT-60SA experiment
14:05 – 14:30 講演 1-2
矢本 (量研)Extension of SONIC code toward mixed-impurity seeding capability
14:30 – 14:55 講演 1-3
本多 (量研)
Development of the integrated model with the iterative solver GOTRESS
14:55 – 15:20 講演 1-4
福山 (京大)Progress in kinetic full wave analyses in fusion plasmas
15:20 – 15:40 休憩(座長：林)
15:40 – 16:05 講演 1-5
登田 (核融合研)
Modeling of turbulent particle and heat transport in helical plasmas based on
gyrokinetic analysis
16:05 – 16:30 講演 1-6
成田 (量研)Particle transport modeling based on gyrokinetic analyses of JT-60U plasmas
16:30 – 16:55 講演 1-7
村上 (京大)Simulation study of toroidal flow generation by ECH in non-axisymmetrictokamak plasmas
16:55 – 17:15 講演 1-8
山本 (京大)Effects of electron cyclotron heating on the toroidal flow in HSX plasmas
17:15 散会
18:30 – 21:00 懇親会 (炙り炉端 山尾 博多駅前にて)1

11 月 30 日(金)9:30 – 9:35 事務連絡
(座長：糟谷)9:35 – 10:00 講演 2-1
Yagi 矢木 (量研)
Revisit ion-mixing mode
10:00 – 10:25 講演 2-2
Kosuga 小菅 (九大)How pattern is selected in drift wave turbulence: role of parallel flow shear(tentative)
10:25 – 10:45 講演 2-3
Gyung Jin Choi (Seoul Univ.)
Gyrokinetic simulation study of parity dependence of magnetic transport
10:45 – 11:05 休憩(座長：登田)11:05 – 11:30 講演 2-4
大澤 (九大)Stable structure of hydrogen in tungsten di-vacancy and its isotope effect
11:30 – 11:55 講演 2-5
沼波 (核融合研)
Kinetic simulations for particle transport of multi-species plasmas in LHD
11:55 – 12:20 講演 2-6
奴賀 (核融合研)
Analysis of energetic particle confinement in LHD using neutron measurement
and Fokker-Planck codes
12:20 – 13:30 昼休み(座長：矢木)
13:30 – 13:50 講演 3-1
前田 (京大)Modelling of heat transport in LHD using neural network with
non-dimensional input parameters
13:50 – 14:05 講演 3-2
凌 (京大)
Integrated simulation study of LHD type fusion reactor by TASK3D
森下 (京大)14:05 – 14:20 講演 3-3
Integrated transport simulation of LHD plasma using data assimilation
14:20 – 14:45 講演 3-4
佐々木 (九大)Chirality of helical flows in plasma turbulence
14:45 – 15:10 講演 3-5
糟谷 (九大)Study of plasma instability by numerical simulations in basic experimental
devices
15:10 – 15:30 まとめ糟谷 (九大)15:30 散会2

JT-60SA実験に向けたプラズマと運転シナリオの統合予測モデリング

Predictive integrated modelling of plasmas and their operation scenarios
towards exploitation of JT-60SA experiment林伸彦 1)、ガルシア ジェロニモ 2)、JT-60SAチーム
N. Hayashi 1), J. Garcia 2), and the JT-60SA Team1) 量研、2) フランス原子力庁 磁場閉じ込め核融合研究所
1) QST, 2) CEA IRFM
Modelling strategy towards exploitation of JT-60SA experiment
JT-60SA is a large superconducting tokamak and will start plasma operations in 2020 well in advance of
ITER operation. The JT-60SA mission is to contribute to early realization of fusion energy by supporting the
exploitation of ITER and by complementing ITER in resolving key issues for DEMO reactors. For the
mission, JT-60SA is equipped with flexibly applicable actuators, i.e., negative- and positive-ion-based
neutral beams (NBs) and electron cyclotron waves, to control heating, current and torque profiles (total input
power Pin=41 MW). With the actuators, JT-60SA explores inductive, advanced inductive and steady-statehigh-beta plasma operation scenarios [1]. Integrated modelling and prediction of plasmas and their operation
scenarios plays a key role to plan experiments and prepare the model validation through the prediction
confirmation in real experiments. For the purpose, not only JT-60SA modelling itself but also the validationof models for JT-60U and JET experiments to develop the best modelling framework to predict, and themodel verification between integrated codes to improve the prediction reliability have been carried out in theclose collaboration between Japan and EU [1,2]. These modelling activities are carried out intensely towards
the start of JT-60SA experiments and the following results were obtained.
Prediction results
The integrated code TOPICS with the CDBM anomalous heat transport model validated by JT-60U/JET
experiments [2] is used for the plasma prediction of high-beta steady-state scenario, which is the mostchallenging scenario due to simultaneous achievement of full current drive, low divertor heat load and so on,with their strong interactions. The previous prediction [2] used a simple NB model without a finite-orbit
effect of fast ions, however, the prediction depends on the accuracy of NB models in the verification between
integrated codes. The prediction is improved using the Orbit-Following Monte-Carlo code OFMC, which isalso verified by the ASCOT code. Pedestal profiles are determined on the basis of EPED1 model.
Temperature profiles are solved, while the electron density profile is prescribed. By the power control with
Pin~26 MW, a steady-state plasma (βN=4.3, HH=1.6) is obtained with an internal transport barrier (ITB) as
shown in Fig.1(a). Even with small power perturbations (+/− 0.2 MW) to the above reference case, ITB
continuously moves outward/inward, βN increases/decreases (Fig.1(b)) and the current becomes
over/inductively driven. Although the plasma is unstable against the perturbations, the long time scale ofdeviations is about the current diffusion time in the ITB region (order of 10 s in Fig.1(b)) and enables tosustain the ITB position and target performance by enhancing/reducing powers every several to ten seconds
without their runaways.
Other integrated modelling for inductive scenarios reveals that the rotation with the neoclassical toroidal
viscosity (NTV) due to the toroidal magnetic field ripple degrades the pedestal height, but it is high enoughto achieve target parameters, and error field correction coils in JT-60SA have the potential to control therotation by changing NTV. The obtained predictions clarify the JT-60SA capability to explore the plasma
scenarios indispensable to ITER and DEMO.
[1] G. Giruzzi et al., Nucl. Fusion 57(2017)085001[2] N. Hayashi et al., Nucl. Fusion 57(2017)126037

Fig.1 (a) Profiles of ion/electron temperatures, safety factor
and prescribed electron density in a JT-60SA steady-state
plasma with βN=4.3. (b) Time evolution of βN in a referencecase (solid line, profiles are shown in (a), input power isadjusted at t=0 s and kept constant to 40 s) and two cases with
±0.2 MW perturbations from t=0 s.

Extension of SONIC code toward mixed-impurity seeding capability
Shohei YAMOTO1), Kazuo HOSHINO2), Tomohide NAKANO1), Nobuhiko HAYASHI1)
1)

QST, 2)Keio Univ.

1. Introduction
The JT-60SA tokamak is now being constructed and is expected to ignite the first plasma at 2020 andto perform ITER-supporting and -complementing experiments from 2023 under various operation
scenarios. One of the underlying operation scenarios is the high-beta operation under full noninductive current drive conditions (steady-state high-beta operation) [1]. The operation will take placewith impurity seeding to reduce the heat load towards the divertors and with keeping the separatrix
density low enough to match the core operation condition. On the other hand, the concentration of the
impurities in the core causes a harmful effect to sustain the high-performance plasma by radiation
cooling and dilution. Therefore, it is important to establish a method to control the impurity transport
in the core and SOL/divertor regions. So far, a favourable spatial distribution of radiation power has
been achieved by Ne + Ar (mixed gas) seeding experiment in JT-60U, i.e. the radiated power is almost
localized in the divertor region and as a consequence, better energy confinement than Ne-only and Aronly cases has been achieved [2]. However, the interpretative and predictive simulation of such mixed
gas seeding operation has not been performed. This is partly because the number of impurity species
that the previous version of integrated divertor plasma transport code SONIC [3,4] could kineticallysolve was limited to one. The issue has been essentially resolved by restructuring SONIC code with
Multiple-Program Multiple-Data (MPMD) framework [5], which allows SONIC to calculate transport
processes of two impurity species by a kinetic impurity transport code IMPMC. The extended SONIC
code was applied to the analysis of the JT-60SA divertor plasma with two impurity species, i.e., theintrinsic C and seeded Ar gas impurities, and demonstrated the radiative divertor plasma scenario [5].
The purpose of this paper is to examine the effects of different impurity seeding species in the JT60SA divertor plasmas step-by-step in order to study the potential impurity seeding operation regime.
Aiming for this purpose, the SONIC code has been further extended to handle three or more impurity
species kinetically based on the MPMD framework mentioned above. Now the SONIC code is capableof calculating the mixed seeding impurities Ne + Ar and intrinsic C transport by IMPMC. The
impurity-impurity interaction such as the physical sputtering of C by Ne and Ar bombardment has
been also implemented. The effects of seeding mixture of Ne and Ar in JT-60SA steady-state highbeta operation scenario is demonstrated by means of the extended version of SONIC.
2. Calculation condition
Figure 1 shows the computational grid for
JT-60SA. The locations of the injection portsof each impurity species are also shown in Fig.
1. The intrinsic C impurities are also
considered in the calculation.
The following steps are set to discuss theeffects of additional Ne seeding taking the Ar
seeding case [5] as a reference; (i) the
parametric survey of the Ne seeding rate (0.01– 0.1 Pa m3/s) with other parameters kept thesame as ref. [5], and (ii) the exploration of theD2 gas puff and the seeding rate of impurities
which satisfies the target parameters of the
scenario (the electron density at the separatrix
ne,sep~1.7×1019 m-3, and the divertor heat load
qdiv<10 MW/m2). The Ar seeding rate is
automatically adjusted to keep the total
radiation power to be 13MW in the plasma.

Fig. 1 Computational grid for JT-60SA. The
locations of the injection ports of each impurity
species are also indicated.

3. Results
3.1 Step (i)Figure 2 shows the fractional radiation power in thedivertor region as a function of Ne seeding rate. As
shown in Fig. 2, the radiation power is more localized
in the divertor region as the Ne seeding rate increases.
This is due to the additional seeding of Ne impurities,
which are mostly radiate in the divertor region.
Furthermore, the Ar radiation power in the divertor
region is also increased because the lower charge state
densities of Ar has been increased in the Ne+Ar seedingcase.
In addition, the electron and D+ densities at the outer
midplane are reduced compared to the Ar seeding caseas shown in Fig. 3. The possible reason is that theradiation power localized in the divertor region leads tothe cooler electron temperature in the divertor region.
In that case, the recycling is enhanced and the neutralpressure becomes high. Therefore, the pumping flux
becomes high and the D+ density is decreased.
Therefore, the Ne+Ar seeding is effective to obtainradiative divertor operation with keeping the upstream
electron density lower than the Ar-only seeding case.
However, the mechanism is not fully understood andthus we will investigate in detail in near future.
3.2 Step (ii)
In order to obtain the Ne+Ar seeding operation
conditions which satisfies the target parameters of the
scenario, D2 gas seeding rate has been increased with
fixed seeding rate of Ne 0.1 Pa m3/s. A set ofparameters (Ar seeding: 0.1 Pa m3/s, Ne seeding: 0.1
Pa m3/s, and D2 seeding: 10.0 Pa m3/s) has been found,which satisfies the target parameters. The radial
distribution of the electron density at the outer
midplane becomes ~1.7×1019 m-3, and the target heatload along the outer divertor plate becomes less than
10 MW/m2 as seen from Fig. 4. The spatial(a)Fig. 2 Fractional radiation power inthe divertor region as a function of Ne
seeding rate.

Fig. 3 Electron and D+ separatrix
densities at outer midplane as afunction of Ne seeding rate.

(b)

Fig. 4 Radial distributions of the plasma parameters; (a) the electron density at the outer
midplane, and (b) the target heat load along the outer divertor platedistributions of the radiation power density for each impurity species are shown in Fig. 5. We havealso compared the obtained result with Ne-only and Ar-only seeding cases with respect to the core
dilution and core radiation. In Ne+Ar case, the fuel purity at the edge (𝑛"# /𝑛%) become highest among
Ar-only, Ne-only, and Ar+Ne cases (𝑛"# /𝑛%~ 0.79, 0.77, and 0.83 in each case, respectively), whichis good for the core performance improvement. In addition, the Ar radiation power in the edge isreduced from 1.1MW (Ar-only case) to 0.6 MW (Ne+Ar case). Therefore, the Ar radiation in the coreis possibly reduced in Ne+Ar case.

4. Summary
The purpose of this paper is to examine the effects of different impurity seeding species in the JT60SA divertor plasmas step-by-step in order to study the potential impurity seeding operation regime.
Aiming for this purpose, the SONIC code has been further extended to handle three or more impurity
species kinetically based on the MPMD framework. Now the SONIC code is capable of calculating
the mixed seeding impurities Ne + Ar and intrinsic C transport by IMPMC. The effects of seeding
mixture of Ne and Ar in JT-60SA steady-state high-beta operation scenario is demonstrated by meansof the extended version of SONIC. The following steps are set to discuss the effects of additional Ne
seeding taking the Ar seeding case [5] as a reference; (i) the parametric survey of the Ne seeding rate
(0.01 – 0.1 Pa m3/s) with other parameters kept the same as ref. [5], and (ii) the exploration of the D2
gas puff and the seeding rate of impurities which satisfies the target parameters of the scenario. In step(i), the radiation more localized in the divertor region can be seen as the Ne seeding rate increases. In
addition, the electron and D+ densities at the outer midplane are reduced compared to the Ar seedingcase. The Ne+Ar seeding is effective to obtain radiative divertor operation with keeping the upstream
electron density lower than the Ar-only seeding case. In step(ii), a set of parameters which satisfies
the target parameters of the scenario have been achieved (Ar seeding rate: 0.1 Pa m3/s, Ne seeding
rate: 0.1 Pa m3/s, and D2 puff rate: 10.0 Pa m3/s). In Ne+Ar case, the fuel purity at the edge (𝑛"# /𝑛%)become highest among Ar-only, Ne-only, and Ar+Ne cases (𝑛"# /𝑛%~ 0.79, 0.77, and 0.83 in each case,respectively), which is good for the core performance improvement. In addition, the Ar radiation
power in the edge is reduced from 1.1MW (Ar-only case) to 0.6 MW (Ne+Ar case). Therefore, the Arradiation in the core is possibly reduced in Ne+Ar case.
References[1] JT-60SA Research Plan
(http://www.jt60sa.org/pdfs/JT-60SA_Res_Plan.pdf)[2] N. Asakura, et al., Nucl. Fusion 49 (2009) 115010.
[3] H. Kawashima, et al., Plasma Fusion Res. 1 (2006) 031.
[4] K. Shimizu, et al., Nucl. Fusion 49 (2009) 065028.
[5] K. Hoshino, et al., Contrib. Plasma Phys. 56 (2018) 638.

[MW/m3]Ar 3.5MW

Ne 2.3MW

C 7.0MWFig. 4 The spatial distributions of radiation power density for each impurity species. The total
radiation power of each impurity species is also indicated.

Development of the integrated model with the iterative solver
GOTRESSM. Honda, E. Narita and N. Hayashi
National Institutes for Quantum and Radiological Science and Technology, Naka,
Ibaraki 311-0193 Japan

Introduction
The steady-state transport equation solver, GOTRESS, has been developed, which
directly finds the solution where a transport flux matches an integral of sources and sinks
using global optimization techniques such as a genetic algorithm and Nelder-Meadmethod [1]. So far, GOTRESS has been used in a manner that an equilibrium and source
profiles are given in advance. For the purpose that GOTRESS is used to examine thefeasibility of prescribed kinetic profiles used in an operation scenario and to develop a
new scenario, however, a consistent solution among predicted kinetic profiles, anequilibrium and source and sink profiles must be necessary. This fact motivates us todevelop a novel integrated model with GOTRESS as a kernel of the model.

GOTRESSIn the light of the aim of the code, GOTRESS is very similar to TGYRO [2]. Bothcodes solve a steady-state transport equation and calculate the kinetic profiles that giverise to the transport flux coinciding with the integral of source and sink profiles. However,they take a completely different approach to get a solution. While TGYRO discretizes thevolume-integral form of the equation and solves it with Newton method and fixed point
iteration, GOTRESS prepares the two kinds of the volume-integral equations and directlytries to find a solution of each equation, such as temperature and its radial derivative,
using global optimization techniques. In doing so, GOTRESS can avoid possible
oscillations originating from the discretization of the equation and the estimation ofderivatives when a stiff transport model applies. Further details of the methodology canbe consulted in [1].

GOTRESS+TGYRO is characterized as a profile solver in the OMFIT framework [3], which is acomprehensive integrated modeling framework that has been developed to enable physics
codes to interact in complicated workflows. As explained above, GOTRESS can naturally
play a same role as TGYRO in an integrated modeling framework. Accordingly, theintegrated model GOTRESS+ has been developed, mainly consisting of the equilibrium
and current profile solver ACCOME [4] and the neutral beam (NB) heating code OFMC[5], other than GOTRESS. The settings of auxiliary heating systems such as NB heatingand electron cyclotron heating (ECH) and the coils that support an equilibrium aredescribed in an input file of ACCOME. The prescribed density and temperature profiles
are also given as initial guess. After the iterative calculations in ACCOME, a consistentsolution is obtained between an equilibrium and a current profile according to the given
kinetic profiles. ACCOME also estimates ECH and EC current drive (ECCD). The
profiles required for OFMC are sent to OFMC, and then the NB heating profiles and NB
current drive (NBCD) are estimated. Now the current density profiles, the equilibrium
and heating profiles have been computed, which are all required for GOTRESS.
GOTRESS predicts the temperature profiles based on them. At this point, the first

Fig. 1: Workflow of GOTRESS+. The predicted profiles are illustrated and each color of a line
corresponds to the number of times of iteration.

iteration of GOTRESS+ is completed, as shown in Fig. 1. The predicted temperature
profiles are sent to ACCOME and replaces the prescribed ones used in the first iteration.
The second iteration then commences. Note that optionally ACCOME can make use of
the NBCD profile calculated by OFMC instead of its own Fokker-Planck solver, but the
choice is up to a user. This iteration continues until the profiles are well converged.

Comparison of the predicted profiles with TOPICSTOPICS has already carried out the feasibility test of the high-bN scenario in JT60SA, so-called #5-1: BT/Ip=1.72T/2.3MA, PNBI/PECH=18.85MW/7MW [6]. GOTRESS+now applies to this case. As shown in Fig. 2(a), bothresults are very good agreement in terms of ion
temperature, even though they come from thoroughly
different integrated models. Good agreement holds forelectron temperature as well. When it comes to theresidual ion power dPi, which is the difference between
the transport flux and the integral of sources and sinks,

dPi of GOTRESS+ is much smaller than that of
TOPICS, indicating that the profiles predicted by
GOTRESS+ approach to a steady state more closely.
Time consumption of GOTRESS+ is 6,812 sec, which ismuch faster than that of TOPICS being 42,909 sec.

References[1] M. Honda, Comput. Phys. Commun. 231 (2018) 94.

Fig. 2: Comparison of (a) Ti
and (b)the residual ion
power between GOTRESS+and TOPICS.

[2] J. Candy et al., Phys. Plasmas 16 (2009) 060704.
[3] O. Meneghini et al., Nucl. Fusion 55 (2015) 083008.
[4] K. Tani, M. Azumi and R. S. Devoto, J. Comput. Phys. 98 (1992) 332.
[5] K. Tani, M. Azumi, H. Kishimoto and S. Tamura, J. Phys. Soc. Jpn. 50 (1981) 1726.
[6] N. Hayashi et al., Proc. Int. Conf. on Fusion Energy (Ahmadabad, India, 2018), PPC/1-1 (2018).

Modeling of turbulent particle and heat transport in helical plasmas
based on gyrokinetic analysis
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

A quantitative prediction of turbulent transport is one of the most critical issues for realizing magnetic fusion energy. Recently, a large number of gyrokinetic simulations of the turbulent transport
in toroidal plasmas have been performed. The gyrokinetic analysis results in tokamak and helicalplasmas have been compared with the experimental observation results. The gyrokinetic simulation for helical plasmas consumes much larger computer resources than those for tokamak plasmas,
because the former requires a large number of mesh points along field lines in order to resolve helical ripple structures. Since it is still not easy to couple the nonlinear gyrokinetic simulation withan integrated transport simulation code, especially for helical plasmas, the predictive model, which
can quickly reproduce the nonlinear simulation results, is highly demanded. The GKV code [1] hasbeen widely used to investigate the ion temperature gradient (ITG) mode and zonal flows in the
Large Helical Device (LHD) for studying the turbulent transport. Gyrokinetic simulations using
the adiabatic electron assumption are performed for the high and the low ion temperature LHD
cases in shot number 88343 [2]. The reduced model for the ion heat diﬀusivity was proposed [3]to reproduce the nonlinear simulation results given by the turbulence simulation with adiabatic
electrons. This reduced model is the function of the linear growth rate for the ITG mode andthe zonal flow decay time. The ion energy flux by this reduced model is in good agreement withthe experimental results for the high-Ti plasmas at t = 2.2s [3]. How to apply the reduced modelof the turbulent ion heat diﬀusivity in the adiabatic electron condition to the transport code has
been shown for helical plasmas [4]. The simulation in the kinetic electron condition shows thelarger ion energy flux than the experimental results for the high-Ti plasmas [5]. On the other
hand, the electron and ion energy fluxes obtained from the simulation with kinetic electrons areclose to those of the experimental results in the low-Ti plasmas [5] at t = 1.8s. The simulation
result with adiabatic electrons in the low-Ti plasmas shows that the ITG mode becomes stablearound ρ(= r/a) = 0.5. This work [6] presents the predictive transport model for the particle andelectron heat diﬀusivities in addition to the ion heat diﬀusivity, including the eﬀect of kinetic electrons on the plasma instability. For the purpose, the gyrokinetic equations for both electrons andions are solved to evaluate the diﬀusivities and fluxes of the heat and particle transport. Whenthe small number for the wavelength was taken, the reduced model for the ion heat diﬀusivity
was proposed [7]. The larger number of the wavelength in the wider wavelength region is taken
than the previous simulation [7], to show the models for the particle and electron heat transport.
The electron and ion heat diﬀusivity models are presented, to reproduce the nonlinear simulation
results in terms of the linear growth rates and the linear response of zonal flows. The quasilinear
flux models are presented to reproduce the particle transport in addition to the heat transportobtained by the gyrokinetic simulations.
The turbulence driven by the microinstabilities in LHD plasmas is studied, using the gyrokinetic
local flux tube code GKV [1]. The electromagnetic gyrokinetic equations are solved for bothelectrons and ions. The temperature and density profiles, and field configuration obtained fromthe LHD experimental results of the high-Ti plasmas at t = 2.2s and of the low-Ti plasmas att = 1.8s and 1.9s [2] are used. The major radii of the LHD plasmas are given by R = 3.75m for thehigh-Ti plasmas and R = 3.6m for the low-Ti plasmas. In the low-Ti plasmas, the magnetic field
configuration is shifted more inward than in the high-Ti plasmas. The generation of zonal flows
can be enhanced in the inward shifted configuration [8]. The electron and ion heat diﬀusivities are
represented in terms of the linear simulation results, L and τ̃ZF asandχmodel
A1e LB1e
e
=
,
B2e
χGB
A2e + τ̃ZF
/LB3e
i

(1)

A1i LB1i
χmodel
i
=
,
B2i
χGB
A2i + τ̃ZF
/LB3i
i

(2)

where the coeﬃcients are given by A1e = C1e CT e
Cz−ξ/b = 1.3 × 10, A2e = C2e CT
Cz−ξ/b =αi +1−c/(2b) −1/(2b)1−c/(2b)2.0, A1i = C1i CT
Cz= 2.6 × 102 and A2i = C2i CT
Cz−1/(2b) = 1.8 × 10. The
exponents are given by B1e = αe a = 0.30, B2e = ξ/b = 0.62, 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. Here, the variables are shown as
αe = 0.19, C1e = 6.8 × 10−2 , C2e = 2.1 × 10−2 , ξ = 0.10, αi = 0.41, C1i = 0.13, C2i = 4.9 × 10−2 ,CT = 6.6 × 10, a = 1.6, Cz = 0.19, b = ∫0.16 and c = 0.27. The quantity L related with themixing length estimate is defined as L ≡ (γ̃k̃y /k̃y2 )dk̃y . The zonal flow decay time is defined by
α +1−cξ/b

⟩

⟩ ⟨

⟨

∫τ

1−cξ/b

τZF ≡ 0 f dtRk̃x (t), where Rk̃x (t) ≡ ϕ̃k̃x ,k̃y =0 (t) / ϕ̃k̃x ,k̃y =0 (t = 0) , the upper limit τf in the
integral is set to τf = 30R/vti and τ̃ZF = τZF /(R/vti ). The normalized electron and ion heatobtained from the nonlinear simulation are compared with theand χ̄i /χGB
diﬀusivities χ̄e /χGB
i
i
GB
model
/χGB
/χi and χmodel
model predictions χei . The heat diﬀusivity models reproduce the nonlineariGBwithin the relative error 0.20.
simulation results χ¯e /χi within the relative error 0.21 and χ̄i /χGBi
The quasilinear models are constructed for both the particle and the energy fluxes. In the quasilinear flux formulation, the particle and energy fluxes are written as
∫
QL

Γ̃

= CΓ

and
Q̃QL
= CQj
j

∫

⟨

Γ̃lin
k̃y
|ϕ̃lin
|2
k̃y

⟨

⟩

Q̃lin
j,k̃y
|ϕ̃lin
|2
k̃y

⟨

⟩

⟩

L 2
|ϕ̃N
| dk̃y
k̃y

⟨

⟩

L 2
| dk̃y
|ϕ̃N
k̃y

(3)

(4)

for the species j, where the quantities with the superscripts lin and N L represent the linear
and nonlinear simulation results. Here, the tilde ˜ represents the normalization of the energy andintensity
particle fluxes by the values of nTi vti ρ2i /R2 and nvti ρ2i /R2 , respectively. The saturated⟩⟨
NL 2of the electrostatic potential fluctuation obtained from the nonlinear simulation, |ϕ̃k̃y | , is wellfitted with the model function of γ̃k̃y /k̃y2 and the zonal flow decay time τ̃ZF ,⟨⟩2 model|ϕ̃k̃y |=Cq1 (γ̃k̃y /k̃y2 )αq1
αZF
Cq2 + τ̃ZF/(γ̃k̃y /k̃y2 )αq2

(5)at each k̃y , where the parameters are Cq1 = 1.0 × 102 , Cq2 = 9.2 × 10−4 , αq1 = 0.54, αq2 = 0.12
and αZF = 1.6. To give the quasilinear flux models Γ̃modeland Q̃modelfor the species j, the model
ql
j,ql

⟨

⟩

L 2
function (5) is substituted into |ϕ̃N
| in Eqs. (3) and (4). When the relative errors of thek̃yfluxes at each k̃y and the total fluxes integrated over the k̃y space are minimized between the
nonlinear simulation results and the quasilinear flux models, the fitting parameters are determinedL
as CQe = 0.78, CQi = 0.58 and CΓ = 0.73. The fluxes from the nonlinear simulation Γ̃N L , Q̃N
and
e
NL
model
model
model
Q̃i are compared with the quasilinear flux models Γ̃ql , Q̃e,ql and Q̃i,ql , where the relative
errors are given by 2.3, 0.24 and 0.24, respectively. The relative error shown above for the particle
flux model is larger than those for the energy fluxes, because the quasilinear particle flux becomesclose to zero at some radial points.
The gyrokinetic equations for both electrons and ions are solved by numerical simulations to modelthe diﬀusivities and the fluxes for the particle and heat transport. First, the electron and ion heatdiﬀusivities are evaluated from the nonlinear simulations for high-Ti and low-Ti plasmas in the
LHD, where the ITG mode is destabilized. The use of the linear simulation results enables us toreproduce the nonlinear simulation results for the electron and ion turbulent diﬀusivities by theheat diﬀusivity models shown. Since the density gradient is close to zero in some radial regionsof the LHD plasmas, the reliable diﬀusivity model for the particle transport can not be shown.
The quasilinear flux models for the electron and ion energy transport are proposed to reproduce
the nonlinear simulation results. The nonlinear simulation results of the electron and ion energy
fluxes are reproduced by the quasilinear flux models at the accuracy similar to that of the heatdiﬀusivity models. In addition, the quasilinear particle flux model, which can be applicable evenfor the flattened density profiles in the LHD, is presented. Thus, the promising transport models,such as the heat diﬀusivity models and the quasilinear flux models for helical plasmas, are proposed
based on the gyrokinetic simulation results. How to apply the ion heat diﬀusivity model to the
dynamical transport code was already reported [4] and the dynamical transport simulation result
for the ion temperature profile will be compared with the experimental result in the LHD. The
study on how to install the electron heat diﬀusivity model and the quasilinear flux models to the
dynamical transport code is in progress and will be reported elsewhere.
This work is partly supported by the NIFS Collaboration Research Programs, NIFS18KNST129,
NIFS18KNXN363 (Plasma Simulator), NIFS18KNTT045, the JSPS KAKENHI Grant Number,
16K06941 and the Collaborative Research Program of Research Institute for Applied Mechanics,
Kyushu University, 30FP-4. This work is carried out using the JFRS-1 supercomputer system at
Computational Simulation Centre of International Fusion Energy Research Centre (IFERC-CSC)
in Rokkasho Fusion Institute of QST (Project code: MTFHP).
[1] T. -H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).
[2] K. Tanaka et al., Plasma Fusion Res. 5, S2053 (2010).
[3] M. Nunami, T. -H. Watanabe and H. Sugama, Phys. Plasmas 20, 092307 (2013).
[4] S. Toda et al., J. Phys.: Conf. Ser. 561, 012020 (2014).
[5] A. Ishizawa et al., Nucl. Fusion 57, 066010 (2017).
[6] S. Toda et al., Phys. Plasmas accepted (2019).
[7] S. Toda et al., Plasma and Fusion Research 12, 1303035 (2017).
[8] T. -H. Watanabe, H. Sugama and S. Ferrando-Margalet, Phys. Rev. Lett. 100, 195002 (2008).

Particle transport modeling based on gyrokinetic analyses of
JT-60U plasmas
E. Narita1 , M. Honda1 , M. Nakata2,3 , M. Yoshida1 , N. Hayashi1 , H. Takenaga1
11

National Institutes for Quantum and Radiological Science and Technology, Naka,
Ibaraki 311-0193, Japan
2
National Institute for Fusion Science, Toki, Gifu 509-5292, Japan
3
The Graduate University for Advanced Studies, Toki, Gifu 509-5292, Japan

Introduction

Fusion power strongly depends on the density and its profile, which are governed by turbulent
particle transport in tokamak plasmas. Peaked density profiles are favorable for an increase infusion power, and they tend to be observed in a low collisionality region [1, 2]. In addition tocollisionality, particle sources due to external fueling such as the neutral-beam (NB) injection
have been found to provide a contribution to density peaking [3, 4], but previous gyrokinetic
studies on particle transport have shown a relationship between density peaking and dominant
instabilities without finite particle sources [5]. Another issue on particle transport modeling is
the computational cost required to predict the density profiles. When one uses the turbulent
transport models based on the first-principles to calculate the turbulent particle flux, it takes
from several hours to days to obtain a stationary-state density profile. In order to tackle the issue,
a neural-network (NN) based approach has recently been undertaken, and NN-based transportmodels are now proposed [6, 7]. This paper presents a new one, which enables us to not only
predict density profiles fast but also investigate eﬀects of the turbulent transport mechanisms,
considering the realistic particle sources.

2

Semi-empirical modeling of the quasilinear particle flux

The turbulent particle flux Γ̄ for species a in a quasilinear limit is given as Γ̄a = D̄(R/Lna +CT R/LTa + CP ), where D̄, R/Lna and R/LTa are the particle diﬀusivity in proportion to thefluctuation amplitudes, the density and temperature gradients, respectively, and all of them are
non-dimensional. The diagonal diﬀusive term is proportional to R/Lna , and the oﬀ-diagonalcomponent consists of the thermodiﬀusive term (CT R/LTa ) and another pinch one (CP ). Applying a semi-empirical method [8], we have calculated CT , CP and D̄, in the electron particle
flux for a dataset of JT-60U H-mode plasmas [2], which are dominated by the ion temperature
gradient (ITG) mode or the trapped electron mode (TEM). First, CT and CP are obtained by
the linear calculations performed with the local flux-tube gyrokinetic code GKW [9, 10]. Next,
assuming that the particle source from the NB injection Γ̄e,NB is balanced with the neoclassicalparticle flux Γ̄e,neo and the turbulent one Γ̄e,turb : Γ̄e,NB + Γ̄e,neo + Γ̄e,turb = 0, we evaluate Γ̄e,turb
referred to as Γ̄e,exp , and estimate D̄ to match Γ̄e,exp . In this way, CT , CP and D̄ are obtained
and they give the magnitude of the particle diﬀusion and pinch terms quantitatively.

3

Analyses of particle diﬀusion and pinch eﬀects

The JT-60U experimental dataset referred to in the previous section has a variety of density peakedness. We haveselected two plasmas from the dataset. The diﬀerence in the
profile of the density gradient between the two cases is shownin figure 1(a). Hereafter, the peaked case and the other one
are referred to as cases A and B, respectively. They have acomparable Γ̄e,exp in the core region due to the approximately
equal external particle source from the NB injection. For thetwo cases, the diﬀusion and pinch terms are calculated using
CT , CP and D̄ obtained in section 2. As shown in figure 1(b),case A has the diﬀusion larger than case B especially aroundρ = 0.5 due to the steeper density gradient. Since the sumof the diﬀusion and pinch terms is balanced with Γ̄e,exp , itis found that the total pinch term of case A works to oﬀset
the large diﬀusion as shown in figure 1(c). Therefore, for caseA, the inward pinch allows the peaked density profile thatdrives the large diﬀusion around ρ = 0.5. Both thermodiﬀusive and another pinch terms of case A contribute to a
reduction in the particle flux in comparison with those of case
B. The thermodiﬀusive term of case A is larger in the inward
flux direction than that of case B due to the steeper electron temperature gradient for case A. On the other hand, theanother-pinch terms of both cases are positive, corresponding the outward direction, and that of case A is closer toFigure 1. Radial profiles of (a)the density gradients with unzero around ρ = 0.5. It has been found that CP is almost
certainty bands, (b) the diﬀuproportional to the ion temperature gradient. Therefore, thesion and (c) total pinch terms
diﬀerence in the another-pinch term between cases A and B isof cases A and B.
mainly due to the ion temperature gradient. These eﬀects ofthe electron and ion temperature gradients on the pinch terms show that the ITG mode for case
A is more suppressed than that for case B, and case A has a more TEM-like instability, resulting in the steeper density gradient for case A. The tendency agrees with the previous study [5],which shows that the transition from the ITG mode to the TEM in the ITG-dominated region
accompanies the increase in the density gradient.

4

Development of a neural network particle transport model

In order to incorporate the understanding of the diﬀusion and pinch mechanisms into prediction of density profiles, CT , CP and D̄ obtained in section 2 are used to train NNs. The
constructed NN predicts CT , CP and D̄ with 11 plasma parameters including R/Ln , R/LTe andR/LTi , in ∼ 0.02 seconds with a single CPU, and it has been validated using a subset that is
not used for NN training, as shown in figure 2. Hereafter, the NN is used as the semi-empiricalparticle transport model to predict density profiles.

Figure 2. Comparison between the value calculated using the semi-empirical method and the
NN regression for (a) CT , (b) CP and (c) D̄. Here, σ is the root mean square error.
Prediction of density profiles isperformed with the integrated code
TOPICS [11], that solves the electron particle transport equation using the turbulent particle flux givenby the semi-empirical particle transport model. The reproducibility ofthe density profile with the developedmodel has been checked for a plasmathat is not included in the dataset for
Figure 3. The predicted density profile for the plasmathat is not included in the experimental dataset used
NN construction. As shown in figfor training of the NNs in comparison with the experure 3, the predicted ne almost falls
imental profile.
within the uncertainty band. In addition to density profile prediction, the diﬀusion and pinch terms are calculated individually. It
has been confirmed that the diﬀerence in the diﬀusion and pinch terms between cases A and B
shown in section 3 can be reproduced by the developed model.

References[1] C. Angioni et al., Phys. Rev. Lett. 90 205003 (2003).
[2] H. Takenaga et al., Nucl. Fusion 48 075004 (2008).
[3] C. Angioni et al., Nucl. Fusion 47 1326 (2007).
[4] T. Tala et al., Proc. Int. Conf. on Fusion Energy (Ahmadabad, India, 2018), EX/4-4 (2018).
[5] E. Fable et al., Plasma Phys. Control. Fusion 52, 015007 (2010).
[6] J. Citrin et al., Nucl. Fusion 55 092001 (2015).
[7] O. Meneghini et al., Nucl. Fusion 57 086034 (2017).
[8] E. Narita et al., Plasma Phys. Controlled Fusion 60, 025027 (2018).
[9] A.G. Peeters and D. Strintzi Phys. Plasmas 11 3748 (2004).
[10] A.G. Peeters et al., Comput. Phys. Commun. 180 2650 (2009).
[11] N. Hayashi and the JT-60U team Phys. Plasmas 17 056112 (2010).

Eﬀects of Electron Cyclotron Heating on the Toroidal Flow in
Helical Plasmas
Yasuhiro YAMAMOTO, Sadayoshi MURAKAMI, Ching-Chieh CHANG, Santhosh T.A. KUMAR1) ,Joseph N. TALMADGE1) , Konstantin LIKIN1) and David T. ANDERSON1)
1)

Kyoto University, Kyoto 615-8450, Japan
University of Wisconsin-Madison, Madison, Wisconsin 53706, USA

The toroidal torque related to the electron cyclotron heating (ECH) is investigated in HSX plasma. Radialdiﬀusion of energetic electrons by ECH makes the ion cancelling current and this ion current generates the j × B
torque, which may play important role in the toroidal rotation observation during ECH in the HSX plasma. We
investigate the energetic electron distribution by ECH applying GNET code, which can solve the 5D drift kinetic
equation for the energetic electron. We evaluate the j × B torque and the the collisional torque due to the frictionof the toroidal drift motion of energetic electrons. As a result, we obtained a significant torque due to the ECH
and that the larger torque is obtained in Mirror configuration than that in QHS one. These results are consistent
with the experiment observations.

1 Introduction
An important role of toroidal flow in the turbulence transport is suggested in many experiments. Recently, spontaneous toroidal flows in the electron cyclotron heating(ECH) plasma have been observed in many tokamak andhelical devices, e.g. JT-60U, LHD and HSX. It has beenrequired to make the mechanism clear and many studies
have been done experimentally[1] and theoretically[2].
The Helically Symmetric Experiment (HSX) is thefirst quasi-symmetric stellarator device. There are two typical configurations of HSX. One is Quasi Helical Symmetry (QHS) configuration, which has the helical direction ofsymmetry in |B|. The (m, n) = (1, 4) mode in Boozer spectrum is dominant in QHS configuration. Another one is
Mirror configuration, where a set of auxiliary coils makes
toroidal mirror terms, (0, 4) and (0, 8) modes, to the magnetic field spectrum in order to break the symmetry.
The flow measurement experiments have been donein HSX with the charge exchange recombination spectroscopy (CXRS)[3]. QHS configuration makes neoclassical viscosity smaller than that of Mirror configuration
because of the helical symmetry, and so we can expectthat toroidal flow velocity in QHS configuration would belarger than that in Mirror configuration. However, smaller
toroidal flow has been observed in the QHS configuration.
It has not been understood well yet.
In this study, we evaluate the the direction and strengthof torques by ECH using GNET code[4], which can solvea linearlized drift kinetic equation in the 5D phase-space,and discuss what makes the diﬀerence.

Fig. 1 The quasi-linear diﬀusion term of ECH.

ing Monte Carlo method. We split the gyrophase averaged
electron distribution function, f , into a stationary part,fMax , and oscillating part by ECH, δ f , as f = fMax + δ f ,where we consider the stationary part is Maxwellian. The
drift kinetic equation for δ f is given by
∂δ f
∂δ f
∂δ f
+ (⃗vd + ⃗v∥ ) ·+ v̇ ·− C(δ f ) − L(δ f )
∂t
∂⃗r
∂⃗v
= S ql ( fMax ) (1)where ⃗v∥ and ⃗vd are the parallel velocity to the magnetic
field and the drift velosity. Also, C(δ f ), L(δ f ) and S ( fMax )are the collision operator, the orbit loss term, and the quasilinear diﬀusion operator as the source term for the absorption of ECH, respectively.

2 Theory and simulation model
In order to study about ECH, we apply GNET code, which
can solve the drift kinetic equation in 5-D phase space us-

1

The ECH driving term is discribed by the quasi-linear
diﬀusion theory. We consider that interaction between the
EC waves and particles are so long that the plasma canreach the steady state. Here we ignored the quasi-linear effect, S ql (δ f ), in the quasi-linear diﬀusion source term for
simplicity. Since we evaluate the torque by ECH in thesteady state, we assume that the density and the temperature are constant. Under these conditions, the source term

v⊥ / vthe

42.0e-113

1.0e-112

0.0e+00

1-1.0e-110-2.0e-11-4(a) QHS-2(b) Mirror-10
1v|| / vthe

2

3

4

(a) ρ ∼ 0.1

Fig. 2 The experimental data of the density and the electron
temperature. The plots are experimental data, and thelines are fitting of them. In simulations, the fitting profiles are used.

v⊥ / vthe

4

S ql is given by[6]S ql ( fMax ) = −δ(⃗r − ⃗r0 ) ·-3

4.0e-1232.0e-1220.0e+00

1-2.0e-12-4.0e-120∂ ql ∂ fMax
D∂vi i j ∂v j

-4

-3

-2

-1

(2)

2

3

4

(b) ρ ∼ 0.3v⊥ / vthe

where Dql
r0 is thei j is the quasi-linear diﬀusion tensor, and ⃗heating position obtained by ray tracing code. The fundamental O-mode ECH is applied in HSX. One typical caseof the quasi-linear source term is shown in Fig.1. It shows
the heating from blue region to red region.
ECH can drive the radial electron current je . The net
current in steady state should vanish to maintain the quasineutrality, so the return current, jb (= − je ) , must flow inthe bulk plasma. Therefore, the bulk plasma feels jb × B
torque.
On the other hand, the supra-thermal electrons drift
toroidal due to the precession motion in the direction opposite to the je × B torque. During the slowing down ofsupra-thermal electrons, they transfer their momentum tothe bulk plasma due to collisions. If we consider the isotopic source, the torque of the passing particles in Co direction cancels that of the passing particles in Counter direction. The trapped particles, however, have the precession motion, and it can make the net collisional torque.
The two torques should cancel in the completely symmetric configuration[5]. However, non-symmetric magnetic
modes enhance the radial electron flux and break the cancellation of the two torques.

42.0e-113

1.0e-112

0.0e+00

1-1.0e-110-2.0e-11
-4

-3

-2

-1

0
1
v|| / vthe

2

3

4

(c) The δ f in the total δ f .

Fig. 3 The velocity distribution of δ f and the radial flux. Theyshow the deviation from Maxwellian.

3 Simulation Results
We perform the simulation assuming typical temperature,density and radial electric field of HSX ECH plasma as
shown in Fig.2. The plasma parameters are as follows: themagnetic axis major radius Rax ∼ 1.2 [m]; the averaged
minor radius a ∼ 0.15 [m]; toroidal magnetic field strength
BT = 1.0 [T]; the absorption power of ECH Pabs = 100[kW].
Applying GNET code, we evaluate the velocity distribution of δ f . Fig.3 (a) and (b) shows the velocity distribution at normarized minor radius ρ ∼ 0.1 and 0.3 sur-01
v|| / vthe

2face. The velocity distribution of total δ f is shown in Fig.3(c). They are the deviation from Maxwellian distribution,
where the red region means the increasing and blue is decreasing. Energetic electrons can be found outer regions
apart from the heating point. This result means that thereis the radial electron flux as shown in Fig.4 .
HSX has the helical symmetry of (1,4) mode, and so
we have to consider 2-dimensional torque as Fig.5. We
cannot compare the size of vectors among three configuration because it is emphasized in order to make it easy
to see. Also, we use the toroidal-poloidal angle as theaxis, not rθ and Rϕ, for simplicity. As a result, the total torque in completely symmetric configuration is almost
perpendicular to the symmetry direction. However, evenin QHS configuration, the j × B torque is much larger thanthe collisional torque, and the total torque has the component parallel to the helically symmetry direction. The direction of the parallel component is the same direction as
the observed flow in the experiments. We consider that thesymmetric conponent is important for the flow. Therefore,
we evaluate the helical force, corresponding to the helicalcomponent of the torque, below.

Fig. 4 The radial flux profile of supra-thermal electrons and theabsorption power density.

Fig.6 shows the helical force in each configuration. In
the completely helically symmetric configuration, the total
force is quite small. Even QHS configuration has the net
force of the symmetry direction. As noted above, the j × B
force is dominant in QHS and Mirror configurations. Also,
the torque in Mirror is larger than that in QHS configuration. This is consistent with the experiments.
In order to investigate the diﬀerence of the helicalforce and the radial orbit among these configurations, we
calculate the collisionless drift orbit of the energetic electron with the energy E = 5keV and the pitch angle λ ∼ 20◦(passing) or 80◦ (trapped). The orbit in each configurationis shown in Fig.7. As seen from Fig.7 (a), the orbits of apassing particle are the same among the three configurations, so the three orbits look like one line. However this is
not the case for trapped particle. In the completely symmetric configuration, the helically trapped particle goesaround and doesn’t move radially. To the contrary, the particle goes radially in QHS configuration and the orbit in
Mirror configuration is much larger than that in QHS configuration. The radial flux can be roughly understood as theradial mean free path, which is determined by the collision
frequency and the radial drift velocity. When the collisionality is low enough for electrons to move along the drift
orbit, larger orbit makes more flux. Therefore, the torqueis much smaller in completely helically symmetric configuration, and that in Mirror configuration is larger than thatin QHS configuration.
QHS configuration has not only the helically symmetric mode but also other non-symmetric modes. Also, Mirror configuration has two large non-symmetric magnetic
modes additionally. The eﬀect of non-symmetric modes isshown in Fig.8. In QHS configuration, 7 non-symmetricmodes are 10 ∼ 20% of (1,4) mode near the axis (ρ ∼ 0.1)
and enhance the radial flux. When we ignore each magnetic mode one by one, the helical force decreases little bylittle. When all of the 7 modes are not including, the profile

Fig. 5 The direction of each torque. (a) is that of the completely
helically symmetric configuration, (b) is that of QHS configuration, and (c) is that of Mirror configuration. The
green dashed vector is the component of helical symmetry direction, and the blue dashed vector is its perpendicular component. The background contour shows themagnetic field strength pattern.

is similar to that of the helical symmetriy case. Therefore
we can consider that there is no especially dominant effective mode among them. In Mirror configuration, Mirror
terms, (0,4) and (0,8) mode, dominate the enhancement of
the radial flux. Except for the two Mirror modes, the helical force has very similar profile to that of QHS case.

43Conclusions

We have evaluated the collisional and j × B torques by
ECH, using GNET code. In axisymmetric and helicallysymmetric configuration, the collisional and j × B torques
almost cancel each other in the direction of symmetry.
QHS and Mirror configuration has the component of thesymmetry direction. This symmetry direction is the same(a) QHS config.

(b) Mirror config.

Fig. 8 The helical force profiles not including several modes (a)and (b).

Fig. 6 The helical total force. They include j × B and collisional
force.

[3] A. Briesemeister , et al., Plasma Phys. Control. Fusion 55(2013) 014002.
[4] S. Murakami, et al., Nucl. Fusion 40 (2000) 693.
[5] M. N. Rosenbluth and F. L. Hinton, Nucl. Fusion 36 (1996)
55.
[6] M.Brambilla, Kinetic Theory of Plasma Waves Homogeneous Plasmas, Oxford Science Publications, 1998(a) Passing orbit.

(b) Trapped orbit in the completely
helically symmetric config.

(c) Trapped orbit in QHS config. (d) Trapped orbit in Mirror config.

Fig. 7 The collisionless orbit of a supra-thermal electron. Fig.
(a) shows the passing orbits of the three configurations
and Fig.(b)-(d) shows each trapped orbit.

with the experimental flow direction, so we consider that
helical force is important. As a result, we obtained thelarger torque in Mirror configuration than that in QHS configuration. This result is consistent with the experimental
observation. It is also found that what makes the dependence is the orbit enhanced by the non-symmetric magnetic mode.
In order to compare the simulation with the experiment, we must solve the momentum balance equations
with viscosity and torque. Now we are tackling the flow
calculation. In the neoclassical prediction, the flow should
be small in the core region where the gradients are small.
However large flow was observed in Mirror configuration.
As a preliminary calculation, we solved them with some
limitations and assumptions. The torque makes the parallel flow to the magnetic field line in the expected direction.
Therefore ECH torque would explain the large flow in Mirror configuration.
[1] M. Yoshida, et al., Nucl. Fusion 49 (2009) 115028.
[2] P. H. Diamond, et al., Nucl. Fusion 53 (2013) 104019.

4

Ion Mixing Mode Revisited
M. Yagi1 , K. Miki1 , N. Miyato1 , and B. Scott†1

QST, Rokkasho Fusion Institute†
IPP Garhing

Introduction
The particle transport is an important research topic for burning control in DEMO reactors. For
example, the hollow density profile is often seen after pellet injection or gas-puﬀ. In such a case,the inverted density gradient appears in the edge region which produces the particle pinch. To understand particle pinch mechanism, we have performed local gyrokinetic simulation using DEFEFI
code[1, 2] with the positive density gradient (dn/dr > 0). We have observed inward particle flux forsuch condition. To understand the mechanism of particle pinch, we have revisited the ion-mixing
mode theory proposed by B. Coppi[3]. It is pointed out that the electron drift wave contributes tothe particle pinch rather than the ion-mixing mode for the inverted density gradient.

Simulation Results
Simulation parameters relevant to plasma edge in ASDEX Upgrade are used: R0 = 165cm, ne =2.0 × 1013 cm, Te = 100eV, Ln = −7cm, LTi = LTe = LT = 3.5cm, Ti /Te = 1, B = 2.5T,
a/R = 0.303, L⊥ /R = LT /R = 0.0212, q = 3.5, β = 6.44 × 10−5 , νi (L⊥ /cs ) = 9.56 × 10−3 ,νe (L⊥ /cs ) = 8.23 × 10−1 . Here we assume Lx = Ly . Grid numbers and time step for simulation
are given as (nx , ny , nz , nw ) = (32, 128, 32, 32, 16) and ∆t = 5 × 10−3 . Figure 1 shows the temporal
evolution of particle flux for the cases with the normal density gradient Ln = 7cm and the inverted
density gradient Ln = −7cm. It is seen that inward particle flux appears for Ln = −7cm. Figure 2
shows the dependence of mode frequency on wavenumber with Ln = −7cm. Here the negative sign
corresponds to the electron diamagnetic frequency, since ω∗e < 0 for Ln < 0.

Revisit of Ion Mixing Mode
A.

Local Analysus

We assume the non-adiabatic electron response is written by iδ model as
ñe
eφ̃
=
(1 + iδ),
n0
Te

1

(1)

Fig.2 Dependence of mode frequency onwavenumber.
Fig.1 Temporal evolution of particle flux.

where

)
(
k
3
δ=
ω − ω∗e + ω∗e ηe ,
ωχ
2

(2)

2
k = 1 + αT , αT the thermal force coeﬃcient and ωχ = k∥2 vthe
/νe . The quasi-linear particle flux isgiven by
(
)
2
∑
∑
2c
eφ̃
e
kθ φ̃kθ ñkθ = −2n0 DBδkθΓ = ⟨ñṽE×B,r ⟩ = − ℑ
,
(3)
B
Te
kθ >0

kθ >0

e
where DB= cTe /eB. If δ > 0, the Γ < 0 which means the inward particle flux.

we obtained the local dispersion relation of the ion-mixing mode including the toroidal ITG
eﬀect as( 2 2)()
k∥ cs 112ϵnK
1 + iδ − + − 2+ bs +
1+
= 0,
(4)
Ω
ω∗e Ω2ΩΩ
where K = (1 + ηi )/τ , bs = (kθ ρs )2 , ϵn = Ln /R, and Ω = ω/ω∗e . If we set bs = 0 and ϵn = 0 andassume K ≫ Ω, then the original model is obtained.

B.

Nonlocal Analysis

Using the ballooning transformation and taking the strong ballooning limit, we obtain the eigenvalue
equation as
[ 2
(
(
))
]
1 ∂2ϕ
∂2
cs
2ϵn
1
2
2
+
Bϕ
+
b
s
+
s
−
χ
ϕ
s
2 Ω2 q 2 R2 ∂χ2∂χ2 ω∗eΩ2k νe ω∗e−i(η̂e + Ω)q 2 R2 ϕ = 0,
(5)
A ve2
where
K
η̂e ≡ 1.5(ηe − 1.5), A ≡ 1 + , B ≡
Ω
2

(

1
1−
Ω

)

1
2ϵn
+ bs +,
AΩ(6)

and −∞ ≪ χ ≪ +∞. This equation is solved by the Hermite expansion and the dispersion relation
is given by( )2 {(
)
}2 (
)
qΩ
1 1
2ϵn
ϵn (2s − 1)2bs
1−
+ bs ++ bs s +
ϵn
Ω A
Ω
Ω
{
}2
k νe ω∗e 2 2
η̂e + Ω
2
× 2n + 1 + iqR= 0.
2n + 1 A ve2
bs s2 + ϵn (2s − 1)/Ω
(7)
For the simulation parameter, K = −1, so that we can derive the simplified dispersion relation. We
consider two limits. In the case of bs Ω ≪ 2ϵn ,()2ν Ω2 (Ω + η̂e )Ω(Ω + 2ϵn )2 + s 2n + 1 + i= 0,(8)
2n + 1 Ω − 1
where
s=

νe ω∗e 2 2ϵ3n2k(2s − 1), ν =
q R .
2
bs q
ϵn (2s − 1) ve2(9)For ν = 0 and bs ≪ 1,
Ω3 + 4ϵn Ω2 + 4ϵ2n Ω + s∗ = 0,(10)where s∗ = s(2n + 1) . This dispersion describes ITG mode with the special limit with K = −1.
If we set s∗ = 0, we have Ω = −2ϵn which implies the toroidal ITG mode is marginal. On theotherhand, if we set ϵn = 0, the slab ITG mode is found to be unstable. In the case of Ω ≫ 2ϵnand bs ≪ 1,
2(Ω + ϵ∗n )(Ω − 1)2 − νe∗2 Ω(η̂e + Ω)2 = 0,

(11)

where
ϵ∗n ≡

ϵn (2s − 1),bs s2

νe∗ ≡ 2k

νe LT /cs me q
, ŝ = (2n + 1)s.
ω∗e LT /cs mi ϵn ŝ

(12)This dispersion implies that the electron drift wave is destabilized by the non adiabatic electron
response though the parallel electron heat conductivity for the inverted density gradient. This is
another branch which can drive the inward particle flux in this limit. Figure 3 and 4 show mode
frequency and linear growth rate given by the general dispersion relation Eq.(7). Here we use kθ ρs =0.3, s = 1, other parameters are same as those used in the nonlinear simulation. For comparison,
those given by Eqs.(10) and (11) are also shown by the label wr(ITG) and wr(EDW), respectively.
Good agreement between Eq.(7) and Eqs.(10), (11) is obtained. Here Ωi < 0 corresponds tounstable mode since ω∗e < 0. From Figures (3) and (4), we conclude the electron drift wave plays acrucial role for the particle pinch for the inverted density gradient limit rather than the ion-mixingmode.
For summary, The ion-mixing mode is revisited in the case of inverted density gradient taking
account of toroidai curvature. It is shown that there exists two unstable modes, namely, the ionmixing mode and the electron drift wave. The ion-mixing mode is driven by negative compression(slab ITG mode) and non-adiabatic electron response though parallel electron heat conductivity.
The toroidal eﬀect is subdominant. For the case with the inverted density gradient, the electron
drift wave is dominant which produces the particle pinch. As future work, ion mode should be
identified which is observed in high kθ regime.
3

Fig.3 Dependence of mode frequency on radialquantum number n.

Fig.4 Dependence of growth rate on radialquantum number n.

References[1] B. Scott, Phys. Plasmas 17 (2010) 102306.
[2] K. Miki, N. Miyato, B. Scott, M. Yagi, Proc. of M&C 2017, Korean Nuclear Society, 2017-01.
[3] B. Coppi and C. Spight, Phys. Rev. Lett. 41 (1978) 551.

4

Role of density modulation in driving nonlinear streamer flows
in drift wave turbulence
Y. Kosuga[1,2][1] Research Institute for Applied Mechanics, Kyushu University, Fukuoka, Japan[2] Research Center for Plasma Turbulence, Kyushu University, Fukuoka, Japan

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,2]. In order to selectively excite one flow or another, we must understandtheir generation mechanisms in detail.
In this work, we discuss the role of density modulation in streamer flows[3]. In theprevious studies, modulational feedback from shearing was discussed. This mechanism isapplicable to both zonal flows and streamers. In this process, shearing flows modulate
underlying turbulence, which modulate the wave number of underlying turbulence. This
leads to the amplification of the modulating flows. For the case of streamers, in addition
to this process, we show that density modulation opens a novel feedback mechanism toamply the original flow. Unlike zonal flows, streamer has a finite poloidal mode number.
As a consequence, streamer flows modulate density field. With this effect, the frequencyof the underlying drift wave turbulence is modulated, which leads to the modulationof the poloidal wave number, very much like the wave number modulation in Langmuir
turbulence. The relative importance of the shearing modulation and density modulation isevaluated, and for typical plasma parameters, the density modulation becomes dominant
effect. With this density modulation feedback, we have revisited the selection rules of
streamers and zonal flows. We find that streamer flows are selected i.) for the medium
value of ρ∗ ∼ 1/200 − 1/300 and small fluctuation amplitude ñ/n0 < 1% or ii.) large levelof fluctuation amplitudes ñ/n0 < 10% and ρ∗ ∼ 1/10. Thus streamer flows can be driven

in the core region of the medium size tokamaks, or the linear machine such as PANTA.
Stimulating discussion with S. Inagaki, F. Kin, M. Sasaki, P.H. Diamond, K. Itoh, S.-I.
Itoh, is acknowledged. This work was supported by Grants-in-Aid for Scientific Researchof JSPF of Japan (JP15H02155, JP17H06089, JP18K03578), the Joint research projectin RIAM.

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

Gyrokinetic simulation study of parity dependence of magnetictransport
G.J. Choi1, B.J. Kang1, T.-H. Watanabe2 and T.S. Hahm1
1

Department of Nuclear Engineering, Seoul National University, Seoul, Korea
2
Department of Physics, Nagoya University, Nagoya, Aichi, Japan

1. Introduction
In the past, physics study of influence of finite plasma 𝛽 on transport of magnetic fusion
plasmas has been focused on linear stabilization of ion temperature gradient (ITG) modes [1,2]and onset of kinetic ballooning mode (KBM) [3] above critical 𝛽. The finite-𝛽 stabilizationof ITG modes is well known to originate from magnetic field line bending [4]. More recently,
considerable attention has been paid to finite-𝛽 effect on transport in fully nonlinear stage [5,6].
Nonlinear electromagnetic gyrokinetic simulations have showed that magnetic heat transportinduced by electromagnetic ITG turbulence make non-negligible modification to heat flux.
Parity analyses of the electromagnetic ITG turbulence have showed that nonlinearly generated
tearing parity modes play an important role in magnetic heat transport [7] by making full
stochasticity of fluctuating magnetic field [8], and thus significantly change 𝛽-dependence ofelectron magnetic heat flux. In this work, we perform a further detailed study of the parity
dependence of magnetic heat transport of finite-𝛽 tokamak plasmas in ITG mode regime,
analyzing results of nonlinear electromagnetic gyrokinetic simulations.

2. Nonlinear magnetic heat flux
Previously performed parity analysis of tokamak transport in ITG mode regime [7] showedthat tearing parity modes, nonlinearly generated from linearly unstable twisting parity ITG
modes, have significant contributions to electron magnetic heat flux. Following previous works,
we analyzed results of nonlinear electromagnetic gyrokinetic simulations with cyclone basecase parameters performed using GKV code [9]. Full-kinetic electron with real mass ratio wasassumed, and 𝛽𝑒 = 0.6% was chosen for the simulations. We obtained nonlinear electron andion magnetic heat fluxes by averaging their amplitudes over the fully nonlinear stage. Fig. 1
shows 𝑘𝑦 components of electron and ion nonlinear magnetic heat fluxes, each decomposed

Fig. 1. GKV simulation results for 𝑘𝑦 components of electron (left) and ion (right) magnetic heat
fluxes (for 𝑘𝑥 = 0 the most linearly unstable 𝑘𝑥 ) in nonlinear stage.

by parity of fluctuating radial magnetic field. First of all, the most linearly unstable mode (𝑘𝑥 =0, 𝑘𝑦 𝜌𝑖 = 0.3) has twisting parity and it induces inward magnetic heat flux, which partiallycancels with electrostatic part and therefore results in reduction of net outward heat flux. Other
twisting parity modes also induce inward magnetic flux. Meanwhile, nonlinearly generated
tearing parity modes induce outward magnetic heat flux, which is unfavorable to tokamak
confinement. Also, note on obvious disparity between electron and ion magnetic heat transport.
The tearing parity modes give significant contributions to electron magnetic heat flux, but theircontributions are negligible for the case of ion magnetic heat flux. Further quantitative analysis

Fig. 2. GKV simulation results for 𝑘𝑦 components of electron (left) and ion (right) magnetic heat
fluxes (after summation over 𝑘𝑥 ) in nonlinear stage.

revealed that tearing parity modes make magnetic heat flux more efficiently than twisting parity
modes for electron transport, while it is the opposite for ion transport. Fig. 2 shows heat fluxes
containing contributions from all 𝑘𝑥 components. It also shows significant contributions from
nonlinearly generated modes on electron magnetic heat flux and negligible contributions forion transport. However, the disparity between twisting and tearing parity significantly weakens.
Further analysis and physical interpretation on this is ongoing.

3. Time evolution of magnetic heat flux

Fig. 3. GKV simulation results for time-evolution of electron (upper) and ion (lower) magnetic heat
fluxes (after summation over 𝑘𝑥 and 𝑘𝑦 ), normalized by square of radial magnetic field.

We also compared nonlinear magnetic heat fluxes to quasi-linear ones. Fig. 3 shows time
evolution of electron and ion magnetic heat fluxes. It is clear that normalized nonlinear heat
flux is more than normalized quasi-linear heat flux for electron magnetic transport. Also, thedirection of electron nonlinear heat flux is opposite to that of quasi-linear heat flux. The trend
is the reverse for ion magnetic transport; normalized nonlinear heat flux is less than normalized
quasi-linear heat flux. Once again, clear disparity between electron and ion magnetic transport
is revealed.

References[1] W. Horton et al, Phys. Fluids 28, 3050 (1985).
[2] B.G. Hong and W. Horton, Phys. Fluids B 1, 1589 (1989).
[3] W.M. Tang, J.W. Connor and R.J. Hastie, Nucl. Fusion 20, 1439 (1980).
[4] F. Jenko and W. Dorland, Plasma Phys. Control. Fusion 43, A141 (2001)[5] S.E. Parker et al., Phys. Plasmas 11, 2594 (2004)[6] M.J. Pueschel and F. Jenko, Phys. Plasmas 17, 062307 (2010)[7] D.R. Hatch et al., Phys. Rev. Lett. 108, 235002 (2012)[8] W.M. Nevins et al., Phys. Rev. Lett. 106, 065003 (2011)[9] T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006)

Stable structure of hydrogen in tungsten di-vacancy
and its isotope effectKazuhito Ohsawa1, Takeshi 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, the divertor armor tiles are covered by tungsten (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 mechanical toughness etc. However, a large
amount of H isotopes (D and T) retention has been reported in the tungsten materials.
In particular, tritium inventory is an important subject associated with safety of fusion
reactors. H isotopes are not retained in perfect crystal lattice of W because of its low H
solubility. However, H isotopes can be accommodated in vacancy-type lattice defects
nucleated by plasma particle irradiation. Stable structures of multiple H atoms trappedin a W mono-vacancy have been reported in our previous work [1]. In the present work,
we pay attention to stable structures of H atoms trapped in W di-vacancy. Besides, theisotope effects are investigated by taking into account zero point energy (ZPE) of H atoms.
2. Methodology
Stable structures of multiple H atoms trapped in the di-vacancy are investigated interms of first-principle calculations on the basis of density functional theory. The
calculations are performed by using Vienna ab initio simulation package (VASP) with
PBE type potential. The cut-off energy of plane wave is 350eV. We use large supercell
composed of 432 lattice points ( 6 × 6 × 6 BCC lattice) in the present calculations in orderto reduce the effects of periodic boundary condition imposed on the supercell. Lattice
relaxations are iterated until resultant force of each atom is lowered than 0.003 eV/Å.

The binding energies of H atoms to the di-vacancy not including ZPE corrections aredefined by(1)
Eb = E [V2 ] + kE H I − E [V2 H k ],I
where E is the cohesive energy of the corresponding supercell; V2, H , and V2Hk indicate
di-vacancy, interstitial H atom, and di-vacancy trapping k H atoms, respectively. The[ ]ZPEs of H atoms are calculated by harmonic oscillation approximation, that is, the forcesacting on the H atoms in stable structures are assumed to be proportional to thedisplacement from the equilibrium positions. The binding energies including ZPE

correction are[ ]EbZ = Eb + kZ H I − Z [V2 H k ] ,(2)where Z means the total ZPE of H atoms in the corresponding supercell.
3. Results
We investigate stable structure of di-vacancy
in W lattice. First nearest neighbor configuration,that is, two vacancies line up in the <111>direction, is most stable. H atoms trapped in the
di-vacancy are located in the vicinity of octahedral
interstitial sites (O-sites) on the inner surface ofthe di-vacancy, as shown in Fig. 1. There are 12 Osites next to the di-vacancy. Besides, the O-sites
are classified into two types. Six O-sites located intheinterface

area

of

the

di-vacancy

arerepresented by red numbers. While, six O-sites
located in the periphery area are expressed by

Fig. 1 Di-vacancy in the firstnearest neighbor configuration andO-sites next to the di-vacancy.

Fig.2 Stable structures of H atoms in W di-vacancy. k indicates the number of H
atoms trapped in the di-vacancy. Each H atoms occupying O-site is represented by
red or blue number corresponding to the O-site. The arrows indicate thedisplacement directions of H atoms from the O-sites.

blue numbers, as shown in Fig. 1. Stable structures of H atoms trapped in W di-vacancy
not including ZPE corrections are shown in Fig. 2. The H atoms preferentially occupy
the O-sites in the interface area when the number of H atoms in the di-vacancy is lessthan five. While, H atoms occupy the O-sites located in the periphery area as H atoms
increase. Stable structures of H atoms not including and including ZPE corrections are
not always identical, as shown in Fig. 3. The dependence of the binding energy on thenumber of H atoms trapped in the di-vacancy is exhibited in Fig. 4. The binding energies
including ZPE corrections E bZ are larger than those not including ZPE corrections E b .

Fig. 3: Stable structures of four H atoms inthe di-vacancy (a) not including and (b)including ZPE corrections, respectively.

Fig. 4: Binding energies of H atoms notincluding E b and including E bZ ZPE
corrections.

4. Discussion
The affiliation of H atoms to the O-sites in the interface area is larger than that in
the periphery area. However, mutual repulsive interaction acts on the H atoms. Actuallyalmost all H atoms are displaced from O-sites so as to keep distance each other, as shownin Fig. 2. The stable structures of H atoms are supposed to be determined by the balanceof the affiliation and repulsive interaction. Furthermore, the stable structures of H
atoms and the binding energies are affected by the ZPEs of H atoms. According to Fig. 4,
we find that E b < E bZ , i.e. ZPE causes the increment of the binding energies. The lighter
H isotopes have larger ZPE. Therefore, the lighter H isotopes are supposed to be trappedin the vacancy-type lattice defects more than the heavier ones.
References[1] K. ohsawa et al., Phys. Rev. B 85 (2010) 184117.

16th Burning Plasma Simulation Initiative (BPSI) Meeting

Kinetic simulations for particle transport ofmulti-species plasmas in LHD
Masanori Nunami
National Institute for Fusion Science, Toki 509-5292, Japan
Using gyrokinetic and drift-kinetic simulations, we investigate the anomalous and neoclassical particle transport of multi-species plasma in Large Helical Devise. In a high ion
temperature plasmas including the carbon impurity ion with hollow density proﬁle, the gyrokinetic simulations within wide ranges of radial gradients of temperatures and densities represent
that the turbulent contributions of particle ﬂuxes are varied satisfying the ambi-polar condition between multi-ion-species and electron. And the turbulent contribution of the impurity
carbon particle ﬂux remains radially inward directed. For the neoclassical contributions, thedrift-kinetic simulations indicate that neoclassical particle ﬂuxes of ions can change due toambi-polar radial electric ﬁeld and the external momentum sources. The external torques of
injected heating beam cause not only the ion-root, but also the electron-root. In the case ofelectron-root, the neoclassical contributions of the particle ﬂuxes may be balanced with the
turbulent contributions obtained from gyrokinetic simulations.

1 Introduction
For designing of fusion reactors, physics of the plasma heat and particle transport is one of
the critical issues in fusion plasma researches. Numerical approaches via gyrokinetics and driftkinetics frameworks are useful to make progress in the researches. In the multi species Large
Helical Devise (LHD) [1] plasmas heated by neutral beam injection (NBI), we often observe theextremely hollow density proﬁle of impurity ions, which is called impurity hole. It is one of
the advantages for fusion reactors because the impurity hole structure may avoid the impurity
accumulations which exacerbate core plasma performances. Therefore, the understandings of the
impurity hole generations are strongly demanded for magnetically conﬁned plasma researches.
In the multi-species plasmas, the particle transport ﬂuxes should satisfy the balance relations
for the species s,()∂ns (ρ)1 ∂(Trb)+ 0 V 0 Γ(NC)+Γ= Ss .
s
s
∂t
V ∂ρ

(1)

If the system is in the steady state, ∂ns /∂t = 0, and the auxiliary particle sources and sinks
can be negligibly small, Ss ∼ 0, the neoclassical contribution of the particle transport ﬂuxes
author’s e-mail: nunami.masanori@nifs.ac.jp

16th Burning Plasma Simulation Initiative (BPSI) Meeting
(NC)

Γs

(Trb)

should be balanced with the turbulent contributions Γs

(NC)
(Trb)
Γs
∼ Γs
.

for each species s, namely

Therefore, in this work, the neoclassical and turbulent contributions to theparticle transport of the plasma are investigated based on the gyrokinetic and drift-kinetic
frameworks.

2 Turbulent particle fluxes
From here, we consider the impurity hole LHD plasma which consists of four species (e, H, He, C).
We perform the gyrokinetic turbulence simulations using the GKV code [3, 4] in order to evaluate
turbulent contributions of the plasma transport. In the previous linear estimates in the plasma,
it was found that the ion temperature gradient (ITG) mode is a dominant instability [5]. From
the nonlinear simulations with wide ranges of temperature gradients of ions and electron within
50% changes from their nominal values, it is found that the particle transport ﬂuxes have quite
diﬀerent temperature gradient dependences between diﬀerent particle species. Furthermore, theimpurity carbon particle ﬂux remains negative, namely radially inward-directed.
In addition to the temperature gradients, the density gradient lengths for several particle
species also aﬀect the particle transport ﬂuxes. Figure 1 shows the density gradient dependencesof the turbulent particle ﬂuxes for each species. In the results, it is found that there are quite
diﬀerent dependences. The density gradients of bulk ions (H and He) have little impacts oneach particle ﬂux. In contrast, the density gradients of electron and carbon ion strongly aﬀect
the turbulent particle ﬂuxes of electron and hydrogen. However, for the carbon particle ﬂux, itcannot be radially outward-directed unless the impurity hole structure is disappeared, in which
the carbon density proﬁle becomes to be peaked.
Therefore, at least in the results from the temperature and density gradient dependences, theturbulent contribution of the carbon particle ﬂux is inward-directed, and the neoclassical particle
ﬂux of carbon is demanded to be outward-directed to balance with the turbulent contribution.

3 Neoclassical particle fluxes
For the simulations to evaluate the neoclassical contributions of the particle transport in the
impurity hole LHD plasma, we employ the DKES/PENTA code [6]. In our previous work [7], due
to the external torque from the NBI heating on the neoclassical ﬂuxes, the co-injected external
torque can cause not only ion-root with negative ambi-polar radial electric ﬁeld Er but also
electron-root with positive Er . In the case of the electron-root, each neoclassical particle ﬂux ischanged from the ion-root case. Especially, the intensities of the particle ﬂuxes are moderated,(NC)and −Γs(Trb)

close to the turbulent contributions Γs

(NC)
(Trb)
Γs
+ Γs
∼ 0.

obtained by gyrokinetic simulations, i.e.,

Therefore, if the impurity hole LHD plasma is in the electron-root, theparticle balance relation, Eq.(1) may be satisﬁed within the kinetic frameworks.

16th Burning Plasma Simulation Initiative (BPSI) Meeting

Fig. 1 Turbulent particle ﬂuxes for changing (a) electron density gradient, (b) hydrogen density gradient,(c) helium density gradient, and (d) carbon density gradient. Here, each ﬂuxes are evaluated at
r/a = 0.6 and normalized by the arbitrary unit.

Acknowledgment
This work is supported in part by the RIAM collaborative Research Program (30FP-12), JSPS
KAKENHI (16K06941 and 18H01202), NIFS Collaborative Research Program (KNTT043, KNTT045
and KNST095), and by FLAGSHIP2020, MEXT within the priority study 6. The results wereobtained by using the K computer at the RIKEN (hp170260, hp180200 and hp180083), JFRS-1
at IFERC CSC (KML3D, GKVISO and GT5DSTL), and ”Plasma Simulator” at 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] M. Nunami, et al., 27th IAEA FEC, TH/P6-8 (2018).
[6] D.A. Spong, Phys. Plasmas 12, 056114 (2005).
[7] M. Nunami, et al., Plasma Fusion Res. 12, 1203039 (2017).

Analysis of energetic particle confinement in LHD using neutron diagnostics
and Fokker-Planck codes
H. NUGA1 , R. SEKI1,2 , K. OGAWA1,2 , S. KAMIO1 , Y. FUJIWARA1 , M. OSAKABE1,2 , M.
ISOBE1,2 , T. NISHITANI1 , M. YOKOYAMA1,2 , and LHD Experiment Group1

2

National Institute for Fusion Science, Toki, JapanSOKENDAI (The Graduate University for Advanced Studies), Toki, Japan

Introduction The physics of the fast ion transport is the one of the important issue for theachievement of the magnetically confined fusion reactor. This is because the plasma temperature is sustained by the kinetic energy transfer from fast ions generated by external heatings andfusion born alpha particles. For this reason, the fast ion transport mechanism has been investigatedin several magnetic confinement devices and simulations. Since the fast ion transport is governedby several mechanisms, however, the accurate prediction of the fast ion transport is diﬃcult. In out
research, we estimate the fast ion confinement time from the measurements and simulations in the
Large Helical Device (LHD) instead of the theoretical approach.
The neutron measurements in deuterium plasmas play an important role to estimate the fast ion
confinement time. In present fusion devices, such as LHD, the fusion reaction between thermal
and non-thermal ions is dominant instead of thermal-thermal fusion reaction. For this reason,
we can survey the fast ion density profile and the fast ion velocity distribution from the neutron
measurements.
To investigate the fast ion confinement time, a series of the short-pulse neutral beam (NB) injection experiment has been performed in LHD. The decay time of the fusion born neutron emission
rate has been observed after the NBs have been turned oﬀ. There are two factors which cause theneutron emission rate decaying. One is the reduction of the fusion cross-section due to the fast ion
slowing down. The other factor is the reduction of the fast ion density due to the fast ion transport. Since the contribution of the fast ion slowing down can be estimated by simulation codes,
the diﬀerences between experiment and simulation results shows the contribution of the fast ion
transport. In this paper, we investigate the fast ion confinement time by using LHD experiment
data and simulation codes.
Experiment scenario Figure 1 shows the typical waveform of this series of experiment. NB
port-through power and energy, electron density and temperature on axis, and the neutron emission
rate are displayed. In this experiment, three tangential NBs (NB#1-#3) and a perpendicular NB(NB#4) are used. Four NBs (NB#1-#4) have the pulse width of 40 ms. Plasmas are sustained by
2.5 MW electron cyclotron heating system.
Analysis tools The neutron emission rate is simulated by using FIT3D code[1, 2] and CONV FIT3D
code. FIT3D calculates the fast ion birth profile including the finite orbit width eﬀect. CONV FIT3D
calculates the slowing down of NB fast ions. This code solves the fast ion slowing down equation
based on the classical Coulomb collision theory. Since the energy distribution of NB fast ions can
1

NB#3

ENB~165keV

NB#4

ENB~ 60keV

8
6
4
2
0

e15

neutron rate[1/s]

port

PNB[MW]
Te0[keV]neutron rate[e14/s]

NB#2

NB#1

ENB~180keV ENB~150keV

ne0[e18m-3]

8
6
4
2
0
4
3
2
1
0
3
2

e14

e13
3.5

4

4.5
time[sec]

5

5.5

1
0
3.5

4

4.5
time[sec]

5

Figure 2: Time evolutions of measured andsimulated neutron emission rate (SN139599).

5.5

Figure 1: Typical waveform of this series ofexperiment (SN139605).
be obtained by CONV FIT3D, the time evolution of the neutron emission rate can be also estimated. The deuterium-deuterium fusion cross section is approximated by the beam-thermal fusion
dominant model derived by Ref. [3]. In addition, the coeﬃcients of the astrophysical S function ischosen from Ref. [4].
As noted above, the decay of the neutron emission rate comes from two factors; the fast ion
slowing down and the fast ion transport. Therefore, the measured neutron decay time τexp
n can bedescribed as:
−1
−1
(τexp
= (τsim
+ τ−1
n )
n )
c ,

(1)

where, τsimand τc indicate the simulated neutron decay time and the fast ion confinement time,nrespectively. Since the present simulation ignores the fast ion transport, the fast ion confinement
time can be estimated from the diﬀerences between measured and simulated neutron decay time.
In the following calculation, it is assumed that the plasma species is pure deuteron, the eﬀective
charge Zeﬀ = 1 and T e = T i . The neutron decay time is not sensitive to the eﬀective charge
though the neutron emission rate is sensitive to the plasma eﬀective charge. This is because thecontribution from the fast deuteron-electron Coulomb collision is dominant to the slowing down of
NB fast ions. From our previous investigation[5], it can be considered that these assumptions arevalid to the present calculation.
Analysis Result Figure 2 is the example of the time evolution of the measured and simulated
neutron emission rate. Since the simulation ignores the fast ion transport, it is found that thesimulated neutron decay time is longer than the measured neutron decay time. Figure 3 shows therelation between measured and simulated neutron decay time. It is found that as simulated neutron
decay time increases, the measured neutron decay time separates from the simulated value due tothe presence of τc . The fast ion confinement time in each beam, which are displayed in fig. 3, cansimsim
be obtained by least-square fitting method in τexpn = τn τc /(τn + τc ). From fig. 3, it is found thatthe perpendicular NB fast ions have a poor confinement than the tangential NB fast ions. This canbe explained by the strong ripple induced particle diﬀusion in LHD.
2

neutron decay time (exp) [sec]

0.4
0.3

SN139591-SN139603: Raxis=3.6, Bt=2.75(CCW)eff

NB1 τc =0.58[sec]
=0.55
NB2
=0.44
NB3
=0.088
NB4

0.2
0.1
00

0.1

0.2
0.3
0.4
0.5
0.6
0.7
neutron decay time (simulation) [sec]0.8

Figure 3: The measured neutron decay times are plotted against the simulated neutron decay times.

Conclusion In this paper, to investigate the fast ion confinement in LHD, a series of experiments
and simulations have been performed. We estimate the fast ion confinement time from the differences between measured and simulated neutron decay time. The fast ion confinement time isobtained by the least-square fitting method. From this estimation, the confinement time of tangential NBs have approximately τc ∼ 0.5 sec and the confinement time of perpendicular NB has
approximately τc ∼ 0.09 sec. In the present paper, we have estimated the fast ion confinement timein single magnetic configuration. The investigation of the dependence of the confinement time tothe magnetic configuration is future works.

[1] S. Murakami, et al., Fusion Technology 27 Suppl., (1995).
[2] P. Vincenzi, et al., Plasma Physics and Controlled Fusion 58, 125008, (2016).
[3] D. Mikkelsen, Nucl. Fusion 29, 1113, (1989).
[4] H. -S. Bosch, G. Hale, Nucl. Fusion 32, 611, (1992).
[5] H. Nuga, et al., Nucl. Fusion 59, 016007, (2019).

3Modelling of heat transport in LHD using neural network with
non-dimensional input parameters
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].
In the previous study, TASK3D is assumed that the transport coeﬃcients are given
as the sum of neoclassical (NC) and turbulent (TB) transport coeﬃcients. The integrated heat and particle transports simulation of NBI heating plasma of LHD had beenpreviously performed and this simulation result had been compared with experimental
result. The result shows that simulation result is good agreement some experimental
result when the turbulence transport model of ion is gyro-Bohm(gradT) model and ofelectron is gyro-Bohm model[1, 2]. However, more experiments are performed in LHD,and simulation result of TASK3D is not good agreement with other some experimental
result. It is considered that there are regions of plasma parameters that can not be reproduced in the turbulence transport models. Therefore, a heat transport model that is ingood agreement with many experimental results of diﬀerent plasma parameters in LHD
is necessary.
In LHD, more experiments are performed than those shots, and a database of experiments using TASK3D-a has been constructed. Studies to model transportation frommany LHD experimental data have been done. By modeling the thermal diﬀusion coeﬃcient χ from the experimental data by statistical method, good agreement with χ ofthe experiment was obtained[3]. In addition, the heat flux Q using neural network(NN)is also predicted for the experimental data of DIII-D[4].
In this paper, we construct NN that predicts direct heat flux Q without considering
neoclassical transport and turbulent transport components. We also verified the parameters aﬀecting the result of the constructed NN using the sensitivity analysis method.

1Modeling of heat flux using Neural Network
Simulation model
we use database in the range of r/a < 0.9
and the database accumulated by TASK3D-a.
This research is a collaborative research with
Dr. Yokoyama and we share the data set[3].
We selected 18 non-dimentional input parameters(Table1). In this study, we evaluated the performance of NN using RMSE and RMSPE for error. We divided the data into training, test, validation data and evaluated NN by cross validation method. First, we learned NN using training data so that RMSE between QNN predictedfrom neural network (NN) and experimental value
Qexp (= Qi , Qe ) is minimized. In order to prevent over fitting to training data, we stored theweights and biases of NN when RMSE becomesminimum on validation data. Finally, we calculated RMSE for test data which we did not use
for training and evaluated the performance of themodel against unknown data. We constructed NN
for electron and ion respectively. nality normalized by Banana-Plateau boundary

Table 1: input parameter of NN
Expression
R/ar/aT e /T i
ne a3
1 dT iT i dρ
1 dT eT e dρ
1 dne
ne dρ

q
dq
dρ
ϵe1.5
ff
∗
νep
∗
νip
∗
νeh
∗
νih
∗
νeb
∗
νib

ρ∗e
ρ∗i

Name
Normalized major radius
Normalized mainor radius
Ion to electron temperature ratio
Normalized electron density
Ion temperature scale length
Electron temperature scale length
Electron temperature scale length
Safety factor
Normalized safety factor shear
Electron collisionality normalized by
Plateau-PS boundary
Ion collisionality normalized by PlateauPS boundary
Electron collisionality normalized by 1/νPlateau boundary
Ion collisionality normalized by 1/νPlateau boundary
Electron collisionality normalized by
Banana-Plateau boundary
Ion collisionality normalized by BananaPlateau boundary
Normalized electron Larmor radius
Normalized ion Larmor radius

Simulation result
Fig.1 on the left is the regression plots of comparing heat flux of NN results and experimental results and shows a good agreement with experimental results. Fig.1 on theright is histograms and shows good statistical agreement of the NN model for train andvalidation data but histograms for test data has larger error than other data. The RMSE
of ion and electron for test data are 0.029 and 0.031. RMSE for training and validation
data is smaller than test data, and the value of NN model prediction is good agreementwith experimental results for all data(training, validation, test).
Fig.2 shows the RMSPE values of the sensitivity analysis results. We have also verified the main variables for predicting heat flux in NN using sensitivity analysis. As
sensitivity analysis, we have verified how the output changes by increasing the value ofthe input parameter by 10 to 40%. As a result of the sensitivity analysis of the ions,
RMSPE is in the order of largeness: R/a, T e /T i , ne a3 , which indicates that this is themain component of NN. The sensitivity analysis of our NN model has shown no strong2

(a) ion(b) electron

Figure 1: On the left is regression plots of
Figure 2: The results of electron and ion sen-comparing heat flux of NN results and experimental results. On the right is histogram of ionand electron

sitivity analysis

dependency of the variable in the gyro-Bohm model and OLS analysis. As a result of
the sensitivity analysis of the electron, RMSPE is in the order of largeness: R/a, ne a3 ,q, which indicates that this is the main component of NN.

Summary
We have studied the heat flux of ion and electron in LHD using NN. The developed NN
model of ion and electron has shown a good agreement with experimental results. The
RMSE for test data of ion is 0.029 and of electron is 0.031. We have also verified themain variables for predicting heat flux in NN using sensitivity analysis. The sensitivity
analysis of our NN model has shown no strong dependency of the variable in the gyroBohm model. We will analyze important variables aﬀecting heat transport by decreasing
input variables and analyzing sensitivity.

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] M. Yokoyama, Plasma and Fusion Research 9 (2014) 1302137.
[4] O. Meneghini,Physics of Plasmas 21 (2014) 060702.

3

Integrated simulation study of LHD type fusion reactor
by TASK3D
T. Shinogi, S. Murakami, S. Maeta
Department of Nuclear Eng., Kyoto Univ., Nishikyo, Kyoto 615-8540, Japan

1

Introduction

Design research of a helical fusion reactor is underway based on the present large helical device(LHD). In the previous studies, they verify the plasma operation scenario for the LHD-type helical
reactor FFHR-d1 in view of MHD equilibrium/stability, neoclassical transport, α-particle energy
loss and impurity eﬀect. Moreover, they identify a self-consistent solution of the plasma operation
scenario which achieves steady-state sustainment of the burning plasma with a fusion gain ofQ ∼ 10 is found in the operation regime that has been already confirmed in LHD experiment.
They study FFHR-d1 with simple turbulent diﬀusion coeﬃcients which is several times neoclassicaldiﬀusion coeﬃcients. It is necessary to study with more accurate turbulent transport models.
In this study, we investigate the possibility of a steady state only by α-particle heating, that is
Q = ∞, in the LHD type reactor. We assume LHD type fusion reactor, LHR4, whose size is fourtimes larger than LHD and whose magnetic field strength is 5[T]. Moreover, we assume turbulent
diﬀusion models which are agree with results of LHD experiment in TASK3D.

2

Simulation model

We assume the LHD type fusion reactor, LHR4, which is 4 times larger in the size scale than thoseof LHD. Table 2 shows the comparisons of the configuration parameters for LHD and LHR4.
We apply the integrated transport simulation code, TASK3D, to study the stated sate of
LHDR4. TASK3D is the integrated transport code for helical plasma based on TASK code andhas been developed in collaboration between Kyoto Univ. and NIFS. TASK3D is composed ofmultiple modules each modeling one physical phenomenon. One dimensional diﬀusion transportmodule (TR module) mainly used in this research calculates heat and particle transports of plasma
by solving one dimensional diﬀusion equations. The diﬀusion coeﬃcient in TASK3D is defined as
TB
the sum of neoclassical diﬀusion coeﬃcient and turbulent diﬀusion coeﬃcient as χs = χNC
s + χs ,where subscript s represents particles. Turbulent transport model has been studied to reproduce
results of LHD experiments. One of models is the gyro-Bohm type for electron, which is given byχe = CgB

LHD
LHR4

Major radius [m]
3.6
14.4

Te ρi
eB a

Minor radius [m]
0.632
2.528

Magnetic field strength [T]
5.0
2.75

Table 1: Comparisons of the configuration parameters of LHD and LHR4

1

and is the gyro-Bohm-gradT type for ion, which isTe ρi
χi = CgBgT
eB a

(

)
∆T
a .
T

These turbulent transport model for electron and ion reproduce well the results of LHD experiments. The other model is the one proportional to the inverse of plasma density, which is givenbyχe,i = C1/n

1019
,
ne,i

where CgB , CgBgT and C1/n are constant which show good match with results of experiments inprevious research. In this study, we investigate the steady state of the LHR4 assuming these two
turbulent models.
For simplicity, we solve only the evolution of the ion and electron temperature, and the plasma
density is fixed in this study. We assume three type of density profiles shown in figure 1.

Figure 1: Plasma density profile A, B and C.（n(0) = 5.0 × 1020 /m3 ）We introduce constant coeﬃcient FNC and FTB for neoclassical heat diﬀusion coeﬃcient andturbulent diﬀusion coeﬃcient. Thus, we assume heat diﬀusion coeﬃcient as
TBχs = FNC · χNC
s + FTB · χs .

FNC and FTB show how much transportations are reduced from the current models. In the future
operation of the nuclear fusion reactor, it is heated up by neutral beam incidence (NBI) heating
and raised, and the temperature is raised by plasma heating. However, in this study, we startsimulation with enough high temperature (11 ∼ 21keV) to sustain nuclear fusion reaction. We do
not consider external heating and only α particle heating is assumed as a heat source (Q = ∞) in
this simulation. We change FNC from 10−3 to 1, FTB from 10−3 to 1, central density n(0) from
1.0 to 6.0 ×1020 [m−3 ], and initial central temperature from 11 to 21keV.

3Results

First, we assume the gyro-Bohm-gradT model. We vary initial temperatures of ion and electron toinvestigate the dependence of central temperature at steady state : Ts(e,i) with initial temperature.
The time evolution of the central temperature reaching the steady state is shown in figure 2.
We assume following parameters; FNC = 1,FTB = 0.1, n(0) = 5.0 × 1020 [m−3 ] and B density
profile. We show the typical case where the steady state is not achieved with same condition as
2

Figure 2: Time evolution of electron
central tmperature changing initial temperatures withFTB = 0.1 (Steady state)

Figure 3: Time evolution of electron andion central temperature withFTB = 0.3

above except FTB . In figure 3, FTB = 0.3 is. We have found that in gyro-Bohm-gradT model,
Ts(e,i) reaches almost the same temperature (about 20keV) with FTB = 0.1 even when the initial
temperature is changed. In the case of FTB = 0.3, we do not achieve steady state against any
initial temperature.
We fix the constant coeﬃcient of neoclassical component FNC = 1 and change that of turbulent
component FTB , n(0), density profiles to investigate dependence of Ts(e,i) with them. Figure 4
shows dependence of Ts(e,i) with FTB and n(0) with B density profile. The solid line represents theelectron central temperature at steady state Ts(e) and dashed line represents that of ions. Figure
5 shows dependence of Ts(e,i) with density profiles with FTB = 0.1. The higher the central density

Figure 4: Density dependence of Ts(e,i)
with FNC = 1 changing FTB

Figure 5: Density dependence of Ts(e,i) in several density profiles (A, B and C) with
FNC = 1, FTB = 0.1

becomes, steady state is also achieved for higher FTB value. We found that Ts(e,i) is highest with
C density profile.
We changed the value of FNC in order to consider the case that we could improve the neoclassical transport more than the current LHD. We fix FTB = 0.1 and change FNC from 10−3 ∼ 1
in three density profiles. Figure 6 shows dependence of Ts(e,i) with FNC and n(0) with B density
profile. The eﬀect of lowering the order of FNC is smaller than the influence of FTB . This isbecause heat transport of electrons is dominated by turbulent transport. We show plots of Ts(i)3with various FNC and FTB in B density profile andn(0) = 5.0 × 1020 [m−3 ] in figure 7. In figure 7,
crossmark represents the condition when the steady state is not reached.

Figure 6: Density dependence of Ts(e,i)
with FTB = 0.1 changing
FNC

Figure 7: Dependence of Ts(i) with FNC
and FTB . The color of themarkers show the temperatureof Ts(i) .Cross marks show unsteady state conditions.

Then, we simulate with 1/n model. We fix FNC = 1 and change n(0) and FTB from 10−1 ∼ 1
in three density profiles. Figure 8 shows the dependence of Ts(e,i) with FTB and density with B
density profile and FNC = 1. In 1/n model, we achieve the steady state FNC = 1 and FTB = 1.
We show the dependence of Ts(e,i) with density profiles in figure 9. we found that Changes in thedensity of near the boundary do not aﬀect the central temperature at the steady state in 1/nmodel.

Figure 8: Density dependence of Ts(e,i)
with FNC = 1 changing FTB(1/n model)

4

Figure 9: Density dependence of Ts(e,i) inseveral density profiles (A, B
and C) with FNC = 1, FTB = 1(1/n model)

Conclusion

We have investigated the possibility of a steady state only with α-particle heating, Q = ∞, in the
LHD type reactor.
We assume the LHD type fusion reactor, LHR4, whose device size is four times larger thanthat of LHD and the magnetic field strength is 5[T] at the plasma center. We assume turbulent
4

diﬀusion models which can reproduce well the results of LHD experiment. We have applied theintegrated transport simulation code, TASK3D, to study the steady state of LHDR4. The diﬀusion
coeﬃcient in TASK3D was defined as the sum of neoclassical diﬀusion coeﬃcient and turbulent
diﬀusion coeﬃcient as χs = χNC+ χTBss , where subscript s represented particles. Turbulent
transport model has been studied to reproduce results of LHD experiments. One of models was
gyro-Bohm-gradT model, the other model was 1/n model. In this study, we have investigated thesteady state of the LHR4 assuming these two turbulent models.
We have introduced constant coeﬃcient FNC and FTB and assumed heat diﬀusion coeﬃcient as
TBχs = FNC ·χNC
s +FTB ·χs . FNC and FTB have shown how much transportations has been reduced
from the current models. We did not consider external heating and only α particle heating was
assumed as a heat source (Q = ∞) in this simulation.
We have found that in gyro-Bohm-gradT model, Ts(e,i) has reached almost the same temperature (about 20keV) with FTB = 0.1 even when the initial temperature has been changed. In
gyro-Bohm-gradT model, we achieved steady state with FNC = 1, FTB = 0.28 in the high density.
On the other hand, in 1/n model, we achieved the steady state with FNC = 1, FTB = 1.
From now on, we will study with finite Q value.

References[1] A.Sagara,et al., Fusion Engineering and Design 89 (2014)[2] T.Goto,et al., Nuclear Fusion Volume57 Number6 (2017)[3] S.Murakami,et al., PlasmaPhys.Control.Fusion57 (2015)

5

BPSI 2018

Integrated Transport Simulation of LHD Plasma
Applying Data Assimilation
Yuya Morishita, Sadayoshi Murakami, Masayuki Yokoyama1 , Genta Ueno2
Kyoto University, NIFS1 , The Institute of Statistical Mathematics2

1

Introduction

Data assimilation is one of forms of machine learning and has been developed mainly in meteorology and oceanology. Data assimilation is to seek to find the optimum combination of a numericalmodel and observation by training numerical model parameters using observed data. Thereby, it
would be possible to estimate the state of a system more accurately.
In this study, we apply data assimilation to the integrated transport simulation (TASK3D) of
LHD plasma. We use TASK3D code [1] as the numerical simulation model and estimate electron
and ion temperatures, Te , Ti . We assume the following turbulent heat transport models[2, 3].
Their coeﬃcients, Ce , Ci are also estimated.
( )( )ρeTe
TB
Electron (gyroBom model) : χe = CeeBa( )( )( ′)
aTi
ρi
Ti
Ion (gyroBom+gradT model) : χTB
= Ci
i
eB
a
Ti
Here, B is the magnetic field strength, ρ is the larmor radius, and a is the plasma minor radius.
Time series of observed data Te (t) and Ti (t) are assimilated by the Ensemble Kalman Filter(EnKF)[4], which is one of the data assimilation methods. We use the experment data of NBI
heated plasma in LHD (shot: 114053) as the time series observed data.

2

Simulation models

In data assimilation, optimizing numerical models with observations is allowed by replacing theconstants, variables and initial values of the models with probability distributions (state vector).
Observations are also replaced with probability distribution (observation noise). Data assimilationis the loop of the two steps: prediction and filter. In prediction step, time evolution of the probability distributions is calculated based on the numerical models. The distributions of prediction
are expanded by system noise, which occurs with time evolution using the numerical models. In
filter step, the prediction is compared to the observation data and the state vector is optimized.
The probability distributions become sharp in filter step. For this cause also, it is necessary to
expand the distributions by system noise. The variance-covariance matrixes of system noise andobservation noise are hyperparameters of data assimilation. Although there are several methodsof data assimilation, we use the EnKF, which can be applied for more general numerical models.
The EnKF is Monte Carlo approximation of the Kalman Filter (KF) [4]. In the KF, linear
systems are considered in time evolution, and all probability density functions are assumed to
be Gaussian. The EnKF is an extended KF so that nonlinear system models can be used byrepresenting probability density functions by ensemble. Time evolution of the probability density
functions is calculated by developing each ensemble member with the numerical simulation models.
1

BPSI 2018

Figure 1: NBI heat deposition evaluated by GNET-TD code.
In this study, we fix the variance-covariance matrix of observation noise, and determine thematrix of system noise using a maximum likelihood estimation. Likelihood is the probability
that the observation data is obtained actually from estimated probability distribution. Maximum
likelihood estimation is to seek the parameter θ that gives the maximum likelihood. The log
likelihood formula for data n = 1 ∼ N is


N
N
m
∑
∑
∑
l(θ) =
logp(yn |Yn−1 ) ≃log 
αn(j)  − N logm,
n=1

n=1

j=1

where, for each ensemble member j,1exp
αn(j) = √(2π)l |Rn |

{

}
1
(j)
(j)
(yn − Hxn|n−1 )T Rn−1 (yn − Hxn|n−1 ) ,2

y is the observation data vector (dimension l), x is the state vector, R is the variance-covariancematrix of observation noise, and H is the observation matrix [4].
We simulate NBI heated plasma of LHD (shot 114053). The heat deposition is evaluated by
GNET-TD code [5] as in figure1. We set the cycle of data assimilation to 40 ms and the numberof ensemble members to 100. The standard deviation of observation noise is fixed to 0.2 keV and
we assume that the system noise for Te is 1.2 times that for Ti , because the discrepancy between
the observations and the estimated values by TASK3D in Te is greater than that in Ti .

3

Simulation results

The result of the maximum likelihood estimation for the system noise are shown in figure 2. From
this results, it is assumed that the proper noise variance for Te is 1.2 keV2 (∼1.1 keV) and thatfor Ce and Ci are 4.5 and 1.5. The variance for Te which is larger than 1.5 keV2 (∼1.2 keV) is notproper, because the noise is so strong that the probability distributions diverge.
The results of the prediction for Te using data assimilation with the proper system noise areshown in figure 3. Improvement of the accuracy about Te estimates can not be seen before 1.0
sec, and the estimates take profiles close to the estimates by only TASK3D. This is because theprofiles of Te are predicted more strongly based on the factors other than the variables in the statevector.
The results of the prediction for Ti using data assimilation are shown in figure 4. In this case,improvement of the accuracy about Ti estimates can be seen. Figure 5 is the results of Ce and Ci
2

BPSI 2018

Figure 2: Log likelihood for various system noise variance (N = 14).

Figure 3: Te profiles of the observations, the estimates using data assimilation and the estimatesby only TASK3D.
estimates averaged in radial direction. There are quite diﬀerences between the estimates and theconventional fixed values [5].

4Summary

We have introduced data assimilation method to integrated transport simulation (TASK3D) of
LHD plasma. The strength of the system noise have been determined by maximum likelihood
method. We have seen improvement of the accuracy about Ti estimates, on the other hand, havenot seen that about Te estimates.
In future, we will include the particle density and the NBI head deposition into the state

Figure 4: Ti profiles of the observations, the estimates using data assimilation and the estimatesby only TASK3D.
3

REFERENCESBPSI 2018

Figure 5: Ce (left) and Ci (right) profiles optimized using data assimilation.
vector, and adjust strength of system noise and observation noise more finely.

References[1] S. Murakami, et al., Plasma Phys. Control. Fusion 57, 054009 (2015).
[2] A. Wakasa et al., THC/P4-29, 23rd IAEA Fusion Energy Conf., Daejeon, Korea, (2010).
[3] A. Wakasa et al., Proc. 39th EPS Conf. and 16th Int. Congress on Plasma Physics (2012)
P2.028[4] T. Higuchi, et al., Data douka nyumon [Introduction to data assimilation]. Tokyo: Asakurasyoten, (2017).
[5] H. Yamaguchi, et al., JPS Conf. Proc. 1, 015045 (2014).

4

Chirality of helical flows in plasma turbulence
M. Sasaki1,2,* , N. Kasuya1,2 , K. Itoh2,3 , S. Toda4 , T. Yamada2,5 , Y. Kosuga1,2 , Y.
Nagashima1,2 , T. Kobayashi4 , H. Arakawa6 , K. Yamasaki1 , A. Fujisawa1,2 , S. Inagaki1,2 ,and S.-I. Itoh1,21Research Institute for Applied Mechanics, Kyushu University, 816-8580, Japan
Research Center for Plasma Turbulence, Kyushu University, 816-8580, Japan
3
Institute of Science and Technology Research, Chubu University, 487-8501, Japan
4
National Institute for Fusion Science, 509-5292, Japan
5
Faculty of Arts and Science, Kyushu University, 819-0395, Japan
6
Shimane University, Graduate School of Engineering, 690?8504, Japan
2

1

Introduction

Helical flows are ubiquitous in nature. In magnetized confined plasmas, helical flows that consists of flows in the parallel and perpendicular to the magnetic field directions are driven spontaneously. The parallel flow is related to toroidal intrinsic rotations [1] and is sometimes driven
by externally injected momentum. The perpendicular flow is related to zonal flows and poloidal
mean flows [2]. These flows contribute to the stabilization of the MHD instabilities and the
suppression of the transport. Flow bifurcations have often been observed; the inversions of thetoroidal flow [1], the poloidal flow [3]. The topology of the helical flow is bifurcated in thesecases, where the topology is sometimes characterized by the chilarity. In such a plasma, multiple
free energy sources for instabilities coexist. The temperature and density gradients drive thedrift wave type instabilities. The inhomogeneities of flows excite Kelvin-Helmholtz instability,
the interchange mode, and the D’Angelo mode. The fundamental processes of the formation of
the helical flow in such a system should be understood.

2

Model

We perform a three-dimensional simulation of turbulence in magnetized cylindrical plasmas [4].
The density and the parallel flow are assumed to have their gradients in radial direction. In
order to account for a situation when the density and the parallel flow are inhomogeneous inspace, the particle and parallel momentum sources are introduced. The intensities of the sources
control the stabilities of the resistive drift wave and the D’Angelo mode. In order to change
the dominant instability, the simulations are performed by changing the density source intensity
shot-by-shot with the fixed parallel momentum source.

3

Results

By the change of the particle source intensity, the bifurcation of the turbulent states is found.
The transitions among three states are observed: the D’Angelo / resistive drift wave are domi1

Fig. 1: Energy of each fluctuation as a function of the intensity of the particle source. The
stream lines of the mean helical flows are shown for three different cases.
nant when the particle source intensity is small / large, respectively. In the intermediate region,
the hybrid mode becomes dominant, which has a characteristics both of the D’Angelo modeand the drift wave, but is a different type of mode [5]. The hybrid mode is driven by both thedensity gradient and the parallel flow shear. The flow topology is also bifurcated, associatedwith the transition of the turbulent states. The turbulent fluctuations drive the mean azimuthal
flows and induce the parallel momentum flux, so that the flow with the structure of co-axial
helix is formed. The each fluctuation energy and the associated helical flows are shown in Fig.1.
Depending on the turbulent states, the spatial patterns of the fluctuation change, which leads
to the change of spatial structures of the momentum transport, and the flow formation. The
results illustrate the generic feature of turbulence flow generation in non-equilibrium magnetizedplasmas.

References[1] P. H. Diamond, et. al., Nucl. Fusion, 53, 104019 (2013).
[2] P. H. Diamond, et. al., Plasma Phys. Control. Fusion, 47, R35 (2005).
[3] S. Inagaki, et. al., Sci. Rep., 6, 22189 (2016).
[4] M. Sasaki, et. al., Phys. Plasmas, 24, 112103 (2017).
[5] M. Sasaki, et. al., Plasma Fusion Res. 12, 1401042 (2017).

2

基礎実験装置のプラズマ不安定性に関するシミュレーション研究
Study of plasma instability by numerical simulations in basic experimental
devicesN. Kasuya, T. Ohno 1, M. Sasaki and M. Yagi 2
Research Institute for Applied Mechanics, Kyushu University, Kasuga, Fukuoka 816-8580,
Japan
1
Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, Kasuga,
Fukuoka 816-8580, Japan
2
National Institute for Quantum and Radiological Science and Technology, Obuchi,
Rokkasho, Aomori 039-3212, Japan
E-mail: kasuya@riam.kyushu-u.ac.jp
Clarifying fundamental mechanisms of turbulence is important in magnetic confined fusion
plasmas [1], and researches on structural formation of plasma turbulence has been progressing
by using linear devices [2]. A linear device has a cylindrical shape without magnetic curvatureand can produce high density plasmas, so is suitable for basic experiments to observe plasma
instabilities in details [3]. Numerical simulations of linear devices have been also carried out,
and a reduced fluid model for drift wave turbulence is used to describe the 3-D turbulence [4].
For understanding structural formation mechanism, identification of the spatial structure of
eigenmodes as the drive of turbulence. In basic experimental devices, the typical spatial scalesof the device and gyro-motion of ions are close to each other, so the finite Larmor radius(FLR) effect is not negligible. A mode with the wave length comparable to the effective ion
Larmor radius can become unstable even in low ion temperature plasmas [5]. Here wedevelop a global simulation code for the ion-temperature-gradient (ITG) instability in linear
devices [6]. For the first step, linear characteristic of the ITG instability [7] is investigatedusing the parameter set of PANTA linear device [8]. The parameter dependency is evaluatedwith a gyro-fluid model to identify influence of the FLR effect.
Fluid models are convenient for analyzing global mode structures in collisional plasmas. The
following set of linearized gyro-fluid equations [9] is used, where the magnetic curvatureterms can be eliminated for the analysis of the cylindrical configuration;ˆ 2  1 dn  / / u/ /  1     0,dt
2  Ln y


(1)

du/ /
  / /  n  T/ /    0 ,
dt

(2)

q 
2
1 

  / /  2u/ /  / /   / /
  ii T/ /  T  ,


L
y
3


n
ˆ2
q 
1 dT
ˆˆ 2  1   ii  / /       1   T/ /  T  ,

 dt
 2
 Ln y 3
1 dT/ /

 dt

 

(3)
(4)

where n is the ion density, u// is the ion velocity, T is the ion temperature, τ is the ratio between
ion and electron temperature at the plasma center,  is the ratio between the density andtemperature gradient length, and the subscripts // and  represent the quantities in the paralleland perpendicular directions to the magnetic field, respectively. The other definitions of theparameters are described in Ref. [6]. The gyro-averaged potential Ψ is represented as Ψ ≡Γ01/2Φ with operators 01/2  1  b / 2 1and b    . Operator b gives the square of the2

perpendicular wavenumber k2, which corresponds to the magnitude of the FLR effect. Two
modified Laplacian operators ̂ 2 and ˆˆ 2 are introduced to include the FLR effects. The
quasi -neutrality relation is given to beb20  n T   1 0   (5)
 1 b 2 
to determine the relation between the density and potential. The Fourier-Bessel expansion isused to calculate the eigenmode structures. The global analysis using this set of equations
considers the boundary condition to determines the radial mode structure, which gives thevalues of the perpendicular wavenumber [6]. Figure 1(a) shows the dependency of the growthrate on τ and η. The ITG instability is unstable when the ratio of the ion temperature gradientto the density gradient η exceeds the threshold value c. The local analysis, as in Fig. 1(b),confirms to reproduce the global analysis result by adjusting the wavenumber to that obtained
from the global analysis. The local model with linearization d/dt →-iω, ∇// → ikz, ∂/∂y → ikθgives the following dispersion relation;
2 

 3

 3    A2 F    F  1   b1   *  2
3 
 2

 A1
, (6)
*
2 2  5 1  3
2 2
 kZ     F     kZ 1  F  A2
0
A
A
3
2
2
3



1
1


where

b 2 
3
A1   b 
A2  1  b1    b1  b2    1  1 ,
4 
2


2 2
b  /2k
b2  2 F .
2 , and * Ln
1  b / 2 

,

F

b/2
1  b / 2,

b1   F,(7)

With the parameters of the linear device the condition for the unstable ITG mode can be
reduced to
A2 > 0.
(8)
For small , c is given as
2 5
c   b .
(9)
3 9
The critical value c changes depending on the magnitude of τ. This dependency is related to
the FLR effect. Without the FLR effect, b = 0, no dependency on  appears as in Fig. 1(c). In
the case of the resistive drift wave instability, break of the Boltzmann relation n  Ψ isimportant for the destabilization. In this model, the FLR effect can break the Boltzmann
relation, as Eq. (5), which affects the instability, but is not strong enough, so this effect alone
cannot make the mode unstable.

In this way, numerical analyses using a gyro-fluid model have been performed to investigate
the FLR effect on ITG instability in cylindrical plasmas. A spectrum code with Fourier-Bessel
expansion has been developed for the analysis of global mode structures. The analyticalexpression of c is obtained from the local dispersion relation. Its  dependency comes fromthe FLR effect. Torus devices are used for fusion researches, and for investigation of global
structures and their dynamics in torus, PLATO project is now progressing to establish a
plasma turbulence observatory [10]. Construction of the platform for analyzing 3-D dynamics

of plasma instabilities with a variety of simulation codes is necessary for comparison between
experiments and simulations.

Fig. 1: Dependencies of the growth rate of the ITG instability on τ and η. The mode withazimuthal mode number 2 and axial mode number 1 is calculated with the experimental
parameters in the PANTA device. The dependencies from (a) the global model and the local
model (b) with and (c) without the FLR effect are shown.
Acknowledgements
This work was partly supported by JSPS KAKENHI Grant Number (JP16K06938,
17H06089), by the collaboration program of NIFS (NIFS17KNST112, NIFS18KNXN373)
and RIAM of Kyushu University.
References[1] P. H. Diamond, et al., Plasma Phys. Control. Fusion 47, R35 (2005).
[2] G. R. Tynan, et al., Plasma Phys. Control. Fusion 51, 113001 (2009).
[3] S. Inagaki, et al., Sci. Reps. 6, 22189 (2016).
[4] N. Kasuya, M. Yagi, K. Itoh and S.-I. Itoh, Phys. Plasmas 15, 052302 (2008).
[5] Y. Miwa, et al., Plasma Fusion Res. 8, 2403133 (2013).
[6] T. Ohno, et al., Plasma Fusion Res. 13, 1401081 (2018).
[7] A.K. Sen, et al., Phys. Rev. Lett. 66, 429 (1991).
[8] G. Hattori, et al., Plasma Fusion Res. 10, 3401060 (2015).
[9] W. Dorland and W. Hammet, Phys. Fluids B 5, 812 (1993).
[10] A. Fujisawa, AIP Conf. Proc. 1993, 020011 (2018).