第 18 回
核燃焼プラズマ
統合コード研究会
18th Burning Plasma Simulation
Initiative (BPSI) Meeting
日時：2020 年 12 月 17 日(木)－18 日(金)場所：オンライン(Ver. 1, 2021 年 3 月 2 日発行)

合同会合第 18 回核燃焼プラズマ統合コード研究会
18th Burning Plasma Simulation Initiative (BPSI) Meetingおよび
核融合エネルギーフォーラムサブクラスター会合
シミュレーションクラスター 炉心プラズマモデリングサブクラスター
プラズマ物理クラスター 定常運転・制御サブクラスター(Ver.1)
日時：2020 年 12 月 17 日(木)－18 日(金)場所：オンライン開催

12 月 17 日(木)(13:00 – 17:00 ポスター講演 Poster presentation)
13:00 – 15:00 ポスター講演 1 Poster 1
コアタイム core time 13:00 – 14:00 Group 1A
14:00 – 15:00 Group 1B
15:00 – 17:00 ポスター講演 2 Poster 2
コアタイム core time 15:00 – 16:00 Group 2A
16:00 – 17:00 Group 2B(開室 start 12:30、閉室 end 17:30)
12 月 18 日(金)9:30 – 9:35

事務連絡 Business announcement(9:35 – 10:55 ポスター紹介 Poster talk)(座長：糟谷)9:35 – 9:55 講演 1-1
糟谷 (九大)Development of integrated transport simulation scheme for impurity control intokamak plasmas
9:55 – 10:15 講演 1-2
本多 (量研)Coupling of EPED1 model to the integrated model GOTRESS+ anddevelopment of JT-60SA operation scenario
10:15 – 10:35 講演 1-3
森下 (京大)1

Data assimilation system based on integrated transport simulation applying a
reduced model of neutral beam injection heating
10:35 – 10:55 講演 1-4
登田 (核融合研)Transport simulation directly coupled with gyrokinetic transport models forhelical plasmas
10:55 – 11:15 休憩(Session Leader：村上)
11:15 – 12:00 議論 Discussion
12:00 – 13:00 昼休み
13:00 – 14:00 小グループ議論 Free discussion(核融合エネルギーフォーラムサブクラスターとの合同会合)
14:00 – 14:10 藤田 (名大)サブクラスター関係連絡事項（定常運転・制御サブクラスターの世話人交代について）
14:10 – 14:40 宮戸 (量研)
IFERC 計算機シミュレーションセンターの現状報告
14:40 – 15:10 本多 (量研)
ITER に於ける統合モデリング活動関係の報告
15:10 – 15:40 花田 (九大)日本の球状トカマクにおける EC によるプラズマ立上げ支援の研究（ITPA 統合運転シナリオグループ活動への貢献）
15:40 – 16:00 休憩
16:00 – 16:30 若月 (量研)
ITPA 統合運転シナリオグループ活動報告
16:30 – 17:00 林 (量研)
サブクラスター活動に関する報告と今後の予定
17:00

散会

2

Poster number
18th BPSI meeting (2020)
Group 1A
P1A-1 登⽥慎⼀郎 (核融合研) S. Toda
Transport simulation directly coupled with gyrokinetic transport models for helical plasmas
P1A-2 鈴⽊航介 (京⼤) K. Suzuki
Fast-ion velocity-space reconstruction using FIDA measurements in LHD plasma
P1A-3 沼波政倫 (核融合研) M. Nunami
Visualization and similarity analysis for data of gyrokinetic simulations
P1A-4 成⽥絵美 (量研) E. Narita
Integrated transport simulations with a neural-network transport model and development of
turbulence saturation rules
P1A-5 持永祥汰 (九⼤) S. Mochinaga
Transport analysis in PLATO tokamak using integrated code TASK
P1A-6 Apiwat Wisitsorasak (Mongkut's Univ. Tech. Thonburi)
Development of extended two-point model for asymmetric scrape-off layer
Group 1BP1B-1 ⽮本昌平 (量研) S. Yamoto
Effects of impurity transport on detached divertor plasma transitions by time-dependent
analysis of integrated divertor code SONIC
P1B-2 梅崎⼤介 (九⼤) D. Umezaki
Effect of large-angle elastic scattering between ions and neutral particles on density profile
on divertor plate in divertor plasma
P1B-3 飯尾太那 (京⼤) D. Iio
Effect of magnetic fluctuations on the toroidal torque driven by electron cyclotron heating
in tokamak plasmas
P1B-4 ⼭本泰弘 (京⼤) Y. Yamamoto
Magnetic configuration and heating location dependences of the toroidal torques by ECH in
LHD
P1B-5 Byung Jun Kang (SNU)
Gyrokinetic studies of fast ion precession driven drift instability in reversed shear burning
plasmasP1B-6 Yong Jik Kim (SNU)
Study of non-local transport events in flux-driven ITG turbulenceGroup 2A
P2A-1 本多 充 (量研) M. Honda
Coupling of EPED1 model to the integrated model GOTRESS+ and development of JT60SA operation scenario
P2A-2 福⼭ 淳 (京⼤) A. Fukuyama
Progress of kinetic full wave analysis in plasmas
P2A-3 太⽥佳吾 (京⼤) K. Ota
Influence of the momentum dependence of radial diffusion on the transport of thermal andenergetic particles in tokamak plasmas
P2A-4 佐々⽊真 (九⼤) M. Sasaki
Extraction of turbulent structures and quantification of energy transfer among structures bysingular value decomposition
P2A-5 Taik Soo Hahm (SNU)
Bounce-kinetic studies on barely trapped particles' effect on trapped electron modes
P2A-6 糟⾕直宏 (九⼤) N. Kasuya
Development of integrated transport simulation scheme for impurity control in tokamak
plasmasGroup 2B
P2B-1 ⽯⽥雅信 (九⼤) M. Ishida
Ion mass number dependence of resistant drift wave turbulenceP2B-2 ⽯⿊裕暉 (京⼤) Y. Ishiguro
Effects of toroidal field ripple on toroidal torque by electron cyclotron heating in JT-60U
plasmaP2B-3 柳原洸太 (量研) K. Yanagihara
Development of ECCD calculation module for quasioptical ray tracing code
P2B-4 森下侑哉 (京⼤) Y. Morishita
Data assimilation system based on integrated transport simulation applying a reduced modelof neutral beam injection heating
P2B-5 ⿊⽥ 侑 (九⼤) Y. Kuroda
Evaluation of ballooning and kink instabilities in PLATO tokamak
P2B-6 庄司悠歩 (京⼤) Y. Shouji
Effects of anisotropic triton birth profile on triton burn-up ratio in LHD plasma

Transport simulation directly coupled with gyrokinetic transportmodels for adiabatic electron response for helical plasmas
S. Toda, M. Nunami, and H. Sugama
National Institute for Fusion Science, Oroshi-cho 322-6, Toki, Gifu, 509-5292, JapanThe quantitative prediction of the turbulent transport in toroidal plasmas is one of the most criticalissues to be solved for realization of fusion energy. 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 ﬁeld lines to resolve helical ripple structures. Since
it is still not easy to couple the nonlinear gyrokinetic simulation with an integrated transport simulation code for helical plasmas, the predictive model, which can quickly reproduce the nonlinear
simulation results, is highly demanded. The heat diﬀusivity models for the electron and ion heattransport, and the quasilinear ﬂux models for the particle and heat transport have been proposed
to reproduce the nonlinear gyrokinetic simulation results in the kinetic electron condition for theion temperature gradient (ITG) mode [1] in the Large Helical Device (LHD). These models arebasically the functions of the linear simulation results, which are the linear growth rate and thecharacteristic quantity for the linear response of zonal ﬂows. The heat diﬀusivity models have been
installed into the transport code for simulating evolutions of the plasma proﬁles in the LHD when
the additional modeling by R/LTi is used [2]. However, the particle diﬀusivity cannot be modeled
for the hollow and ﬂattened density proﬁles, which are typical in the LHD experimental results [2].
Furthermore, the quasilinear ﬂux model is not directly coupled with the transport simulation sofar, because the used quasilinear ﬂux model depends on the cross-phase between the ﬂuctuating
potential and temperature (or the density) and it is considered that this cross-phase term cannot
be modeled by the temperature or density gradients.
In this study, the reduced model [3] for the ion heat diﬀusivity in the adiabatic electron condition is
directly coupled with the transport simulation by performing the GKV [4] linear simulation. The
linear growth rate is evaluated at each time step in the transport simulation. The zonal ﬂow decay
time, which is characteristic for the linear response of zonal ﬂows, is dynamically set as the initialstate in the transport simulation. The neoclassical diﬀusion coeﬃcient is derived from DGN/LHD
database with the low-β limit (β = 0). The reduced model for the ion heat diﬀusivity is shown asχmodelA1,ad Lαad
i
=
,
(1)
A2,ad + τ̃ZF /L1/2
χGB
i
)
( ∫
2
2
)d
k̃
is the quantity related
/
k̃
v
/R)
is
the
gyro-Bohm
diﬀusivity
and
L
≡
(γ̃
(=ρ
where χGB
y
ti
y
i
i
k̃y
with the mixing length estimate. Here, γ̃k̃y (=γk̃y /(vti /R)) is the normalized linear growth rate of−1
the ITG mode. The coeﬃcient is given by A1,ad = 1.8 × 10, ⟨A2,ad = 5.2 ×⟩ The⟩ 10⟨ and αad = 0.38.
linear zonal ﬂow response function is deﬁned by Rk̃x (t) ≡ ϕ̃k̃x ,k̃y =0 (t) / ϕ̃k̃x ,k̃y =0 (t = 0) . Notethat the zonal ﬂow response function for k̃x = 0.25 is used to evaluate the representative values ofthe zonal ﬂow decay time, because there are peaks of the wavenumber spectra around k̃x = 0.25

∫τ
in the nonlinear simulation results. The zonal ﬂow decay time is deﬁned by τZF ≡ 0 f dtRk̃x (t),where the upper limit τf in the integral is set to be τf = 25R/vti in the adiabatic electron condition
for modeling the diﬀusivities or ﬂuxes. The correlation time of the turbulent sources is shorter
than 25R/vti . Therefore, the zonal ﬂow response function for τf > 25R/vti is not considered
to inﬂuence the generated zonal ﬂow level. The normalized zonal ﬂow decay time is deﬁned asτ̃ZF = τZF /(R/vti ).
The simulation result for the adiabatic electron response of the stationary ion temperature proﬁlecan be obtained. The linear growth rate is calculated by the GKV at each time step in the
dynamical transport simulation by TASK3D and the value of the ion heat diﬀusivity is obtained.
In the adiabatic electron condition, the simulation result for the Ti proﬁle show the good agreementwith the experimental result. The stationary proﬁles of the turbulent and neoclassical diﬀusivities
are studied. The ITG mode is destabilized in the radial region 0.20 < ρ < 0.80. The turbulent
transport is dominant compared with the neoclassical transport in the radial region 0.20 < ρ <0.80, where the positive electric ﬁeld is chosen. The transport simulation results for the ion
temperature proﬁles coupled with the reduced model are compared and have the good agreementwith those for the additional modeling of the linear growth rate by R/LTi .
This work is partly supported by the NIFS Collaboration Research Programs, NIFS20KNST168(Plasma Simulator) and the Collaborative Research Program of Research Institute for Applied
Mechanics, Kyushu University, 2020FP-11.

[1] S. Toda et al., Phys. Plasmas 26, 012510 (2019).
[2] S. Toda et al., Journal of Physics: Conference Series 561, 012020 (2014).
[3] M. Nunami, T. -H. Watanabe and H. Sugama, Phys. Plasmas 20, 092307 (2013).
[4] T. -H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).

18th Burning Plasma Simulation Initiative (BPSI) Meeting

Visualization and similarity analysis for data ofgyrokinetic simulations
Shuta Kitazawa1 and Masanori Nunami2
1 Nagoya University, Nagoya 464-8601, Japan,
2 National Institute for Fusion Science, Toki 509-5292, JapanThe gyrokinetic simulation is a powerful tool to analyze the turbulent transport of themagnetically conﬁned plasmas. However, since the data obtained from the gyrokinetic simulations consist of time series of ﬁve-dimensional distribution functions, it is diﬃcult to fully
capture the physics in the turbulent transport phenomena. While the direct visualization of
the distribution function has been quite important for the intuitive analysis for the phenomena,in this work, a novel visualization tool and a novel data analysis way for the GK simulation
data are developed with the aid of the latest visualization and data science techniques. Since
the developed way enables us to visualize the higher-dimensional data in the virtual space andto capture the characteristic structure of the distribution functions, it becomes easy to analyze
the complicated the gyrokinetic simulation data.

1 Introduction
In magnetically conﬁned fusion plasma researches, the gyrokinetic simulation is reliable way
to predict the turbulent transport. In the last decade, it becomes possible to validate thegyrokinetic simulation results against the experimental observations, and to predict the proﬁlesof the plasma temperature and density by performing the simulations with changing the input
parameters of the plasma proﬁles [1]. Since the gyrokinetic simulation can treat the evolutionsof the distribution function in ﬁve-dimensional (5D) phase space, we can obtain full data ofthe time-series of 5D distribution functions. On the other hand, it is quite diﬃcult to fully
capture the physics in the simulation data because of the hugeness of the 5D data. Therefore,
the direct visualization and development of the data analysis scheme for the 5D distributionfunctions obtained from the gyrokinetic simulations are signiﬁcant and strongly demanded forunderstandings of the turbulent transport physics in the magnetized plasmas. In this work, wedevelop a novel visualization tool and a novel data analysis scheme using the latest techniques forthe visualization and data science. Using the developed schemes, it may be possible to analyze
the complicated the gyrokinetic simulation data more intuitively.
author’s e-mail: nunami.masanori@nifs.ac.jp

18th Burning Plasma Simulation Initiative (BPSI) Meeting

2 Mapping of velocity distribution function into wavenumber space
In our previous work [2], we developed the visualization tool for the gyrokinetic simulation datausing the head mounted display (HMD) which is a display device with small display optic infront of each eye, and the Unity [3] which is a cross-platform game engine. The developed tool
enables us to visualize the data of 5D distribution functions from the gyrokinetic simulations in
3D real-space with 2D velocity-space distribution functions obtained by the gyrokinetic code,
GKV [4]. However, if we have performed it without any reduction of the simulation data for
5D distribution functions, it is extremely hard to visualize them for the sake of the hugeness ofthe data. Therefore, in actuality, we cannot fully capture the physics in the simulation data by
the conventional way. In order to reduce the ﬁle size of the distribution function data for thevisualization, we employ the image data for the 2D velocity-space part of the 5D distributionfunction data, because the image data is lighter than the full data of the distribution functions.
Then we collect whole images of the 2D velocity-space part of the distribution function data at
each position in 3D real-space, and map the image into the 3D real-space. Using the mapping
data, we are able to visualize the data of 5D distribution functions with quite light load. The
ﬁgure 1 shows an example of the visualization, and we can plot the mapping of 2D velocity-space
distribution functions into the 2D space (kx , ky ) which is the wavenumber space perpendicular
to the magnetic ﬁeld at ﬁxed the remaining real-space coordinate z. In the ﬁgure, in addition
to the 2D velocity-space structures of the distribution functions, we can easily ﬁnd that there
exists the certain global structures in higher dimensional phase space.

Fig. 1 Mapping of velocity-space distribution functions into wavenumber space. At each position inwavenumber space, the corresponding velocity-space distribution functions obtained by the gyrokinetic simulation are fully mapped.

18th Burning Plasma Simulation Initiative (BPSI) Meeting

3 Similarity analysis
Using the mapping data of the distribution functions obtained in the previous section, we performthe similarity analysis to extract the features of the distribution functions from the gyrokinetic
simulations. In the analysis, we employ the numerical library, OpenCV [5], which is a library of
programming functions mainly aimed at real-time computer vision. Here, we focus on the dataof the turbulence simulations for the “inward-shifted” case of the LHD [6] plasma that has animportant feature, i.e., the enhanced zonal ﬂows which cause the plasma transport reduction.
For the mapping data, we ﬁnd that there exists robust structures in the distribution functionsin the wavenumber region where the zonal ﬂows are most enhanced. In the region, there exists
anti-symmetrical structures for vk -direction, and V-shaped structure for vk -v⊥ space. Therefore,
we regard such characteristic structures as a reference structure of the distribution functions,and ﬁnd the similarities in whole mapping data using similarity analysis functions in OpenCV
library including the feature detection technique, AKAZE [7]. The AKAZE is the accelerated
version of KAZE [8] which is a 2D feature detection and description method that operates
completely in a nonlinear scale space. Using the function, we can detect similarities of structure
to the reference structure in the zonal-ﬂow enhanced region with the anti-symmetrical structuresfor vk -direction and the V-shaped structure for vk -v⊥ space. And we ﬁnd that the similarities
remains for long simulation time. This means that it may be possible to extract the importantregion in the gyrokinetic simulations at an early simulation phase.

4 Summary
The schemes to visualize and analyze the gyrokinetic simulation data of the distribution functions
in 5D phase space have been newly developed. By employing the image data for 2D velocityspace distribution functions instead of the full data of 5D distribution functions, and mapping theimage into the remaining 3D space with the wavenumber space, we can visualize the data of the
distribution functions with quite light load. Using the mapping visualization, furthermore, we
can extract the characteristic structures of the distribution functions by means of the similarity
analysis.

Acknowledgment
This work is supported in part by the RIAM collaborative Research Program (2020FP-1), NIFS
Collaborative Research Program (NIFS20KNST161), and JSPS KAKENHI (20K03907). The
simulations were performed by using JFRS-1 at IFERC CSC (GDKTHEL) and ”Plasma Simulator” at NIFS.
[1] T. Görler, et al., Phys. Plasmas 21, 122307 (2014).
[2] S. Kitazawa, M. Nunami, et al., The 81st national convention of IPSJ, 2M-05, Fukuoka, Japan (2019).

18th Burning Plasma Simulation Initiative (BPSI) Meeting[3] Unity oﬃcial web site, https://unity.com/ .
[4] T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24 (2006).
[5] OpenCV oﬃcial web site, https://opencv.org/ .
[6] Y. Takeiri, et al., Nucl. Fusion 57, 102023 (2017).
[7] P. F. Alcantarilla, J. Nuevo, and A. Bartoli, Fast Explicit Diﬀusion for Accelerated Features in
Nonlinear Scale Spaces, British Machine Vision Conf. (BMVC), Bristol, UK (2013).
[8] P. F. Alcantarilla, A. Bartoli, and A. J. Davison, KAZE Features, Eur. Conf. on Computer Vision(ECCV), Fiorenze, Italy (2012).

Integrated transport simulations with a neural-networktransport model and development of turbulence saturation rules
E. Narita1 , M. Honda1 , M. Nakata2,3 , M. Yoshida1 and N. Hayashi111National 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

Since density and temperature proﬁles severely aﬀect fusion power, accurate predictions ofthese proﬁles are required. The gyrokinetic-based turbulent transport models, such as TGLF [1]have been found to give reasonable agreement with experiments, but introducing the detailed
descriptions of turbulence physics has tended to increase the computational cost. To accelerate
the proﬁle predictions, a neural-network (NN) based approach has been undertaken (see e.g.,[2–4]). Our study is also developing a NN based turbulent transport model DeKANIS basedon a combination of the gyrokinetic code GKW [5] and JT-60U experimental data. DeKANIS
has originally been developed for particle transport [6], and then it has extended to cover bothparticle and heat transport [7]. Since the previous model was hard to perform stable simulations
with integrated codes, the NN model has been improved. Furthermore, the prediction capabilities of two diﬀerent turbulent saturation models adopted in DeKANIS have been checked using
JT-60U and JET plasmas.

2

Improvements on the NN model for integrated simulations

DeKANISﬂuxes for electrons( predicts turbulent) particle and heat
) in a quasilinear limit as
(
R
R
R
R
Γ̄e = D̄ Ln + CT LT + CP and Q̄e = χ̄e CN Ln + LT + CHP , respectively. Here, noneeeedimensional D̄, χ̄, R/Ln and R/LT are the particle and heat diﬀusivities in proportion to
ﬂuctuation amplitudes and the density and temperature gradients, respectively. The ﬂuxes arecomposed of the diagonal (diﬀusion) and oﬀ-diagonal (convection) terms, which are quantitatively determined by calculating the oﬀ-diagonal-term coeﬃcients (CT,P,N,HP ) and the diﬀusivities (D̄ and χ̄e ). Note that the relationship between CT,N , D̄, and χ̄e satisﬁes the Onsager
symmetry. Predicting diagonal and oﬀ-diagonal terms individually could be helpful to understand the eﬀects of the transport processes on density and temperature proﬁle formation. In
addition to the electron ﬂuxes, DeKANIS also predicts the ion heat ﬂux, which is given asχ̄i,eﬀχ̄e,eﬀ LRT nnei TTei . Here, χ̄eﬀ is the normalized eﬀective heat diﬀusivity, which is deﬁned
Q̄i = χ̄e,eﬀ
iwhen the heat ﬂux is expressed with only the diagonal term. If Q̄e is given, Q̄i can be calculated
by determining χ̄i,eﬀ /χ̄e,eﬀ .
DeKANIS estimates CT,P,N,HP , D̄, χ̄e and χ̄i,eﬀ /χ̄e,eﬀ with a fully-connected feed-forward NN
model. The training dataset for the NN construction has been developed with 23 JT-60U H-mode
plasmas. The values of CT,P,N,HP and χ̄i,eﬀ /χ̄e,eﬀ have been given by GKW linear calculations.

On the other hand, D̄ has been determined by the two diﬀerent saturation models. The ﬁrst
model calculates D̄ to match it to the the experimental turbulent particle ﬂux. The case using
the ﬁrst model is referred to as DeKANIS-1. The second model related to DeKANIS-2 will beshown in later. The NN model has been constructed to reproduce the training dataset. Here,the inputs are 12 plasma parameters, such as R/Lne and R/LTe , and the outputs are the six
coeﬃcients: CT,P,N,HP , D̄ and χ̄i,eﬀ /χ̄e,eﬀ . After obtaining the coeﬃcients, DeKANIS computes
χ̄e not to break the Onsager symmetry.
The training dataset has been prepared in the same way as the previous DeKANIS [7], but
the number of the data points has increased from 2,232 to 7,268. As for the hyperparameters of
the NN model, the number of hidden layers and unis has been determined using Bayesian optimization with a Python library Hyperopt [8,9] through a Keras wrapper Hyperas [10], whereas
it was optimized manually for the previous model. Furthermore, the activation function of theoutput layer has changed from the linear function to the softsign one to limit outputs withintree standard deviations of the average in the training data. Setting the upper and lower limits
prevents unreasonable extrapolation, which can cause unstable transport simulations.
Due to theaforementioned
modiﬁcations,
DeKANIS-1 hasrealized the density and temperature proFigure 1. Radial proﬁles of (a) the electron density and the (b) electron and (c)ﬁle predictions
ion temperatures predicted by DeKANIS-1 (blue) and -2 (red) using TOPICS
using the infor a JT-60U plasma outside the training dataset.
tegrated code
TOPICS [11]. Figure 1 demonstrates that DeKANIS-1 works for a JT-60U H-mode plasma thatis not included in the training dataset. It takes about three hours to obtain the steady-state
proﬁles with a single core, which is much faster than the conventional models like TGLF.

3

A theory based turbulence saturation model for applicability
to a wide range of plasma parameters

As the ﬁrst saturation model directly
uses the experimental value, DeKANIS1 is pragmatic in the
known parameter range.
However, it is diﬃcult to use outside
the known parameter
range.
The second
saturation model has

Figure 2. The root-mean-square errors σ of the (a) electron density
and the (b) electron and (c) ion temperatures predicted by DeKANIS1 (blue) and -2 (red) using TOPICS for 14 JT-60U test cases at eachradial point.

been introduced to apply DeKANIS to devices other than JT-60U. It is based on the mixinglength theory, assuming that D̄ can be expressed by the quasilinear diﬀusivity γ̄/k̄θ2 and the()1.5 −2.71
residual zonal ﬂow level LZF [12]: D̄ = 1.86 × 10−4 γ̄/k̄θ2
LZF . To use the scaling formula,the other NN model has been constructed for DeKANIS-2 to reproduce γ̄/k̄θ2 and LZF that havebeen calculated by GKW. As the NN model of DeKANIS-2 does not need to output D̄, it isfree from the experimental values. The modiﬁcations made for the DeKANIS-1’s NN model
have also been performed for the DeKANIS-2’s. The prediction results from DeKANIS-2 for
the JT-60U H-mode plasma are similar to those from DeKANIS-1 (ﬁgure 1), and it has beenfound that the the root-mean-square errors σ of the theirfor the 14 test cases are at the√ results
1 ∑similar level (ﬁgure 2). Here, σ of ξ is given as σξ = N (ξsim /ξexp − 1)2 , where ξsim(exp) is
the simulated (experimental) value of ξ, and N is the number of the data points.
The two DeKANIS models have been
applied to JET plasmas, whose data havebeen taken from the
ITPA International
Multi-Tokamak Proﬁle Database [13].
Figure 3. The root-mean-square errors σ of the (a) electron density
We have estimatedand the (b) electron and (c) ion temperatures predicted by DeKANISΓ̄e , Q̄e and Q̄i for 16
1 (blue) and -2 (red) using TOPICS for 16 JET plasmas at each radial
L- and H-mode plaspoint.
mas, and evaluatedσ, as shown in ﬁgure 3. It can be found that DeKANIS-2 gives better predictions, though
the accuracy is worse than that for JT-60U. Since the accuracy of DeKANIS-2 can be improvedby replacing the mixing-length-like scaling formula with another one, the saturation model will
be updated near future.

References[1] G.M. Staebler et al., Phys. Plasmas 12 102508 (2005).
[2] M. Honda and E. Narita, Phys. Plasmas 26 102307 (2019).
[3] O. Meneghini et al., Nucl. Fusion 57 086034 (2017).
[4] K.L. van de Plassche et al., Phys. Plasmas 27 022310 (2020).
[5] A.G. Peeters et al., Comput. Phys. Commun. 180 2650 (2009).
[6] E. Narita et al., Nucl. Fusion 59 106018 (2019).
[7] E. Narita et al., Proceedings of 17th Burning Plasma Simulation Initiative Meeting (2019).
[8] J. Bergstra et al., Proceedings of the 30th International Conference on Machine Learning,
28 (2013).
[9] http://jaberg.github.io/hyperopt/[10] http://maxpumperla.com/hyperas/[11] N. Hayashi and JT-60 Team, Phys. Plasmas 17 056112 (2010).
[12] E. Narita et al., Plasma Phys. Controlled Fusion 60 025027 (2018).
[13] The ITER 1D Modelling Working Group, Nucl. Fusion 40 1955 (2000).

Transport analysis in PLATO tokamak using integrated code
TASK
S. Mochinaga1), N. Kasuya1,2), A. Fukuyama3),
Y. Nagashima1,2) and A. Fujisawa1,2)
1)

1.

IGES, Kyushu University, 2) RIAM, Kyushu University, 3)Kyoto University

Introduction

PLATO device at Kyushu University is a tokamak to observe detailed spatio-temporal structure of plasma turbulencewith several diagnostics, as tomography systems with more than 1500 detectors and three HIBPs [1,2]. The PLATO
device is designed to have plasma parameters as major radius 𝑅 = 0.7 m, minor radius 𝑎 = 0.2 m, plasma current 𝐼= 40 kA, electron temperature 𝑇 = 0.2 keV, and density 𝑛 = 1.0 1019 m-3. A number of integrated codes have beendeveloped, and are applied to evaluation of performance of several tokamak plasmas [3-5]. This paper presents
transport analyses using the integrated simulation code TASK (Transport Analysis System for tokamaK) [6] to predictexperimental conditions. TASK code has a modular structure, which includes various functions to calculate particle
and heat transport, equilibrium, wave heating, etc. A self-consistent simulation is realized by modules selected inaccordance with the research purpose, which are coupled with a data exchange interface. In the simulation ofmagnetically-confined plasmas, transport calculation is performed under the assumption that plasma is in the
MagnetoHydroDynamics (MHD) equilibrium state. Conventional transport analyses with TASK code have been
carried out by using a fixed-boundary two-dimensional equilibrium code, TASK/EQ, and a one-dimensional transport
code, TASK/TR. It is important to consider the external coil conditions to represent more realistic experimental
situations in PLATO. Therefore, in this study, we have developed a routine to include a free-boundary twodimensional equilibrium code TASK/EQU, which calculates the plasma equilibrium with vertical field coils, toanalyze transport in PLATO plasmas. TASK/EQU code has been also used for the analyses of MHD instabilities in
PLATO [7]. The developed transport analysis routine is used to predict the plasma performance, and dependences onparameters that are externally controllable in experiments are evaluated. PLATO device is no direct heating tools forions, so is investigated how to improve ion temperature by power balance analysis in transport processes.

2.

Transport analysis modelTime evolutions of plasma quantities are obtained by coupling equilibrium and transport calculations. In theequilibrium calculation, the following Grad-Shafranov equation is solved.
𝐼 𝑑𝐼
𝑑𝑝
1
(1)
4𝜋 𝜇
.
𝛻 𝛻𝜓
𝑑𝜓 𝑅 𝑑𝜓𝑅
Here, 𝜓, 𝑝, 𝐼, 𝑅, 𝜇 are the poloidal flux function, the plasma pressure, the poloidal current, the major radius, and thepermeability, respectively. In the free-boundary equilibrium calculation, 𝜓 includes contributions both from the∑ 𝜓 , where 𝑖 is the index number ofas 𝜓 𝜓
plasma current 𝜓 and the external vertical coil current 𝜓
the coils.
In the transport calculation, the following particle and heat transport equations are solved for plasma species s.
1𝜕 𝑉𝛤
1 ∂
(2)
𝑛𝑉
𝑆 ,
𝑉
𝑉 ∂𝑡
𝜕𝜌
1∂ 𝑉𝑄
1 ∂ 3
(3)
𝑛 𝑇𝑉 /𝑃 ,/ ∂𝑡 2𝑉∂ρ
𝑉
Here, 𝑛 , 𝑇 , 𝑆 , 𝑃 , 𝜌 are the plasma density, the plasma temperature, the particle source, the heat source, the𝑑𝑉/𝑑𝜌 is the radial derivative of plasma volume. The heat sources
normalized minor radius, respectively, and 𝑉
for electron and ion are given by
𝑃
𝑃
𝑃 ,
(4)
𝑃 𝑃 ,
(5)
where 𝑃 is the ohmic heating power expressed as
(6)
𝑃
𝜂𝑗 ,
and 𝑃 is the collisional energy transfer from electrons to ions expressed as
3 𝑇 𝑇
𝑛
,
(7)
𝑃
2
𝜏
Here 𝜂 is the plasma resistivity, 𝑗 is the plasma current density, and

.

.

√

𝜏

(8)

is the electron-ion collision time. In the particle transport equation, the quasi-neutrality condition is assumed.
Transport coefficients includes neoclassical and turbulent contributions. For the neoclassical transport model,
NCLASS [8] is used to obtain the diffusion coefficients and convection velocities of particles and heat. The plasma
resistivity is also obtained by using NCLASS. For the turbulent transport model, the CDBM model [9] is used toobtain the thermal diffusivity. The transport coefficients are the sum of the neoclassical and turbulent thermal
diffusion coefficients.
We have developed a scheme for linking equilibrium module EQU and transport module TR. The TASK/TOT
module is used to control the time evolution of the modules. In the equilibrium and transport linkage scheme, each
module exchanges data via data interface BPSD. Quantities regarding to the plasma shape, metric and magneticsurface averages are obtained from the equilibrium calculations, and are used in the transport calculations. On theother hand, in the equilibrium calculation, the pressure profile obtained from the transport calculation is used as
an initial condition, and the poloidal current profile is derived from the safety factor q in the right-hand side of Eq.
(1). In this study, change of the current distribution on the flat-top of PLATO plasma discharges is considered tobe small and is not reconstructed in the evolution.

3.

Transport analysis with PLATO parameter

The developed scheme (EQU-TR combination) is used toevaluate the plasma performance and transport property with
PLATO plasma parameters. Toroidal magnetic field 𝐵 = 0.3 T
and plasma current 𝐼 = 40 kA are kept to be constant, and thedensity at the plasma center is set to be 𝑛 = 1.0 1019 m-3 bysupplying particles from the edge of the plasma as particle
supply with the gas-puffing method. PLATO device has main
three pairs of vertical field coils; D (dipole), Q (quadropole),and H (hexapole) coils, and two pairs of feedback coils; F1 and Fig.1. Equilibrium magnetic configuration of PLATO
F2 coils. Ohmic heating coil in the center solenoid is divided
tokamak obtained by (a)TASK/EQU andinto two pairs; OH and OHC coils as in Fig. 1 [10]. In this(b)TASK\EQ. The position of each vertical fieldcalculation, only D, Q and H coils are considered for control.
coils and the ohmic coils are also shown.
Figure 1 shows the poloidal cross-section of the flux contoursof the equilibrium with EQU calculation. The obtained equilibrium has 𝑅 = 0.7 m, 𝑎 = 0.18 m, ellipticity 𝜅 = 1.9, andtriangularity 𝛿 = 0.4. The plasma profiles in the steady-state are obtained with TR calculation as in Fig. 2. The
obtained temperature has 𝑇 = 130 eV and 𝑇 = 35 eV at the center. The conventional transport analysis method (EQTR combination) is also carried out using the same shaping parameters. Similar results are obtained both with EQ
and EQU, but the electron temperature in the EQU-TR case is lower than the EQ-TR case, because the ohmic heating
power near the center of the plasma is different as in Fig. 2(d). The current profile is radially broader in the EQU case,so the current is smaller near the center with restriction of the same total plasma current. In this calculation, theneoclassical convection directs inward, but contribution of the diffusion is much larger. The turbulent diffusion isdominant except near the center, where the density and temperature profiles are flat.
0.15

0.04

𝑇

0.10

0.05

0.02
0.0

0.2

0.4

𝜌

0.6

0.8

1.0

0.3

(c)

0.0
0.0

0.2

0.4

𝜌

0.6

0.8

1.0

0.15

0.05

0.2
0.1

(b)
0.00

0.0

EQU-TR
EQ-TR

0.20

0.4

𝑇

(a)

0.00

0.25

EQU-TR
EQ-TR

0.5

MWm

0.06

EQU-TR
EQ-TR

𝑗 MAm

0.08

𝑇 keV

𝑛 10 m0.6

EQU-TR
EQ-TR

0.10

𝑃

0.12

0.10

(d)

0.00
0.2

0.4

𝜌

0.6

0.8

1.0

0.0

0.2

0.4

𝜌

0.6

0.8

1.0

Fig.2. Radial profiles of (a) density, (b) electron and ion temperature, (c) current density, (d) ohmic heating power. The caseswith EQU-TR and EQ-TR are shown.

4.

Dependences on externally controllable quantities

The dependence on the particle feeding rate is evaluated, as it is externally controllable in PLATO experiments. The
plasma parameters and coil currents are same with those in the previous section. The particle feeding rate is adjusted
to give the density at the center from 𝑛 = 0.5 to 1.5 1019 m-3. Figure 3(a) indicates that the electron temperature
decreases and the ion temperature increases, as the density increases. This dependency can be explained with thepower balance. From Eq. (3), the power balance of electrons and ions are presented with the ohmic power, energy
transfer from electrons to ions and heat fluxes in the steady state as
Electron

:

Ion

:

𝑃

𝑃
𝑃

0,
0.

(9a)
(9b)

Figure 3(b) shows the values at 𝜌 = 0.3, where the temperature gradient is large, for the explanation. A larger density
leads increase of electron flux and 𝑃 , which makes the electron temperature smaller. Equations (7) and (8) give(10)
𝑃 ∝𝑛 𝑇 . 𝑇 𝑇 .
Electron-ion collision is responsible to the energy exchange between electrons and ions, which tends to makedifference of the temperature between electrons and ions smaller. Increases of 𝑛 and 𝑇 . and decrease of 𝑇 𝑇
determines the magnitude of 𝑃 , as is indicated with the red dashed line in Fig. 3(b). Therefore, the ion temperature
increases, as 𝑃 increases. In addition, the larger plasma resistivity gives 𝑃 increase. Equations (6), and thespitzer's formula for the resistivity,
𝜂(11)give the following relation(12)
𝑃 ∝𝑇 . 𝑗 .
The deviation from this relation indicated by the green dashed line in Fig. 3(b) comes from density dependence ofneoclassical resistivity by NCLASS routine.

Fig.3. (a) Relationship between volume-averaged temperature and density, when the particle feeding rate is changed. (b)Dependency of terms at ρ = 0.3 in the power balance on 〈n〉.

The other transport analysis is performed with variation of an external coil current. The plasma parameters are samewith the previous sections excepts D-coil currents, which are changed in the range from -18.8 to -23.3 kA turn, while
the other coil currents are same. The plasma size 𝑉 is changed in accordance with the D-coil current. In the transportmodel of the TR code, the effect of plasma shape is introduced by the plasma volume, metric 〈|∇𝜌|〉 and 〈|∇𝜌| 〉. As
shown in Fig. 4(a), the density slightly decreases with increase of plasma volume. Particle source S(r) decreases as
the total particle feeding rate is constant, but the density variation is slight, due to decrease of the particle fluxes.
On the other hand, the larger the plasma shape is, the lower the electron temperature and the higher the ion
temperature become. Figure 4(b) shows the values at 𝜌 = 0.3 for explanation of the dependency as well. Since theplasma current is kept constant, the current density decreases in accordance with the increase of the plasma volume.
Power 𝑃 is represented as𝑃 ∝𝑇 . 𝐴 ,13where 𝐴 is the size of the poloidal cross-section of the plasma. Therefore, 𝑃 is decreased depending on 𝐴 as inFig. 4(b), which gives the decrease of the electron temperature. For the ion temperature, 𝑃 decreases slightly as
𝐴 increases, but the energy confinement time becomes longer, so the ion temperature increases due to theimprovement of confinement.

Fig.4. (a) Relationship between volume-averaged density, ion and electron temperature, and the plasma volume. This is the casewhen ID-coil is changed. (b) Dependency of terms at ρ = 0.3 in the power balance on plasma volume.

5.

Conclusion

In this study, the transport analyses were performed with calculation of free-boundary condition equilibria,
considering external vertical coils in PLATO experiments. The plasma performance was predicted in PLATO. The
dependences of plasma performance on the particle feeding rate and vertical field coil current, which are externallycontrollable parameters in experiments, were evaluated. From these evaluations, it was found that higher density andlarger plasma volume give larger ion temperature. Change of energy transfer from electrons to ions, or energy
confinement time is the causes for the ion temperature increase. In the initial phase of PLATO experiments, there isno direct ion heating method, so these indirect methods are possible ways to increase the ion temperature.
In the present code, toroidal current profiles obtained by transport calculation is not reflected on equilibrium
calculation, so the process is necessary to be introduced for self-consistent calculation. In addition, prediction
accuracy of transport will be improved by comparing with experimental results in future works.

References[1] A. Fujisawa, AIP Conf. Proc. 1993, 020011 (2018).
[2] T. Ido, et al., submitted to Rev. Sci. Instrum. (2021).
[3] M. Honda et al., J. Plasma Fusion Res. SERIES 6, 160 (2004).
[4] N. Hayashi and JT-60 Team, Phys. Plasmas 17, 056112 (2010).
[5] Y. Murakami et al., J. Nucl. Mater. 313-316, 1161 (2003).
[6] https://bpsi.nucleng.kyoto-u.ac.jp/task/[7] S. Tomimatsu, et al., Plasma Fusion Res. 15, 1403052 (2020).
[8] W. Houlberg, et al., Phys. Plasmas 4, 3230 (1997).
[9] M. Yagi et al., Phys. Fluids B 5, 3702 (1993).
[10] K. Matsuoka, RIAM Reports, Kyushu Univ. 141, 51 (2011) (in Japanese).

Eﬀect of Magnetic Fluctuations on the Toroidal Torque
Driven by Electron Cyclotron Heating in Tokamak Plasmas
D. Iio, S. Murakami, Y. Yamamoto
Department of Nuclear Engineering, Kyoto University, Nishikyo, Kyoto, 615-8540, Japan

1

Introduction

It is well recognized that toroidal plasma ﬂow and its shear play a signiﬁcant role in magnetohydrodynamic(MHD) mode stabilization and energy conﬁnement improvement. Particularly, self-generated toroidal ﬂow
in Electron Cyclotron Heating (ECH) plasma has been observed and analyzed in many tokamaks, such as
JT-60U and DIII-D. The ECH driving mechanism of the spontaneous plasma toroidal ﬂow needs clarifying.
Although many experimental studies[1,2] have been done, it has not been completely understood yet.

2

Simulation model

It is known that two toroidal torque components contribute to the plasma toroidal ﬂow, i.e., j × B
torque and collisional torque. The j × B torque is a consequence of radial ion current in response to
the outward ﬂux of energetic electrons with ECH. Collisional friction between energetic electrons and bulk
plasma generates the collisional torque. In axisymmetric magnetic ﬁeld conﬁguration such as ideal tokamak,these two components are of the same magnitude but in diﬀerent directions, so that total toroidal torquebecomes almost zero. However, it is shown that the balance of these two components breaks, and the net
torque is generated in non-axisymmetric conﬁguration such as helical devices,[3]. Also, in tokamak devices,
it is predicted that the net torque is generated with a small non-axisymmetric component such as toroidal
ﬁeld ripple and resonance magnetic perturbation (RMP).
In this study, we consider the magnetic ﬂuctuations in which toroidal and poloidal modes of several tens
are assumed, and investigate the behavior of the energetic electron and the toroidal torque generated by
ECH. We use the GNET code, which solves a drift kinetic equation for energetic electrons with ECH in 5D
phase space
) ∂δf
∂δf
∂δf (
+ ⃗vd + ⃗v∥ ·+ v̇ ·− C(δf ) − L(δf ) = S ql (fMax )
∂t
∂⃗r
∂⃗v

with the Monte Carlo method[4]. Magnetic ﬂuctuation is considered in the phase of solving the drift equationof motion. We assume that the ﬂuctuation ﬁeld component δB is perpendicular to the equilibrium ﬁeld
component B0 . The ﬂuctuation ﬁeld is deﬁned as
B = B0 + δB
δB = ∇ × αB0 ,where α is the function at real space, which describes the shape of the perturbation at the real space. We
deﬁne the function α asα(ψ, θ, ϕ) = ε∑mn

bmn (ψ) cos(mθ − nϕ + χmn )[](m − nq(ψ))2.
bmn (ψ) = exp −4w2

3

Simulation results

We calculate the toroidal torques under the perturbations with the mode numbers of (m, n) = (10−30, 10).
In the axisymmetric conﬁguration(ε = 0) these two torques are canceled each other. As seen in the Fig. 1,
as the perturbation strength increases, the j × B becomes dominant and the total torque becomes non-zero.
We also calculate the toroidal torques under 2x, 3x and 4x perturbation mode numbers of (m, n) = (10−30).
The obtained result is shown in Fig. 2. It is conﬁrmed that as mode numbers increase, the toroidal torques
become stronger, especially in the co-direction. According this, it is expected that we can use high mode
perturbations for generating toroidal torque instead of strong perturbations. This result indicates thepossibility of a signiﬁcant role on a net toroidal torque generation by high-mode magnetic microturbulence.

Fig. 1: Total torque at each εFig. 2: Total torque at each modes[1] M. Yoshida, et al., Phys. Rev. Lett. 103, 065003 (2009)[2] J. S. deGrassie, et, al., Phys.Plasmas 11 4323 (2004)[3] Y. Yamamoto, et al., Plasma Fusion Res. 14 3403105 (2019)[4] S. Murakami, et al., Nucl. Fusion 40 693 (2000)

˜

Magnetic Configuration and Heating Location Dependences of the
Toroidal Torques by ECH in LHD
Yasuhiro YAMAMOTO and Sadayoshi MURAKAMI
Department of Nuclear Engineering, Kyoto University, Kyoto, 615-8540, Japan
Toroidal torque by electron cyclotron heating (ECH) is investigated in LHD plasmas, assuming the suprathermal electrons by ECH generate torques on the plasma through jr × B and collisions. The jr × B torque depends
on the radial drift velocity and the fraction of trapped electrons. Therefore, the magnetic configuration andthe heating location aﬀects the toroidal torque. We investigate the magnetic configuration and heating location
dependences of the toroidal torques by ECH in LHD, considering typical three magnetic configurations: theinward shifted, standard, and outward shifted configurations. As a result, magnetic ripple bottom heating makes
larger torque than that of ripple top heating because of the large fraction of trapped electrons. Also, heating at theouter minor radius makes larger toroidal torque than that of heating at the inner radius, and the injection angle
can also change the toroidal torque profile. Moreover, ECH makes the largest toroidal torque in the outward
shifted configuration. Finally, With the obtained toroidal torques, we evaluate the toroidal flow velocities withthe obtained toroidal torques. We obtained the largest flow near the axis in the standard configuration because of
its small viscosity and large toroidal torque.

1 Introduction
Toroidal flow is one of the important issues in turbulence transport and MHD stability. Recently, spontaneous
toroidal flows driven by electron cyclotron heating (ECH)
have been observed in many tokamaks and helical devices
such as JT-60U, LHD, and HSX. Many experimental [1, 2]
and theoretical [3] studies have been undertaken to clarify the underlying mechanism. However, the mechanismof the toroidal flow generation by ECH has not yet been
understood well.
We previously reported that ECH could apply torques
on the plasma through jr × B and collisions[6, 7]. Also,
the jr × B torque overcomes the collisional torque in the
non-symmetric configuration. The jr × B torque would bea candidate for the torque driving the toroidal flow.
In LHD, toroidal flows have been investigated in the
Neutral Beam Injection (NBI) heating and ECH plasmas, where the toroidal flow velocity is measured by
the charge exchange recombination spectroscopy (CXRS)[4, 5]. Also, LHD has the flexibility on the magnetic configuration by shifting the magnetic axis. It is known thatthe confinement of energetic particles of the inward shifted
configuration is better than that of the outward shifted configuration due to the trapped particle orbit improvement.
Therefore, we expect the jr × B torque would be significant in the outward shifted configuration because the radialdiﬀusion of trapped electrons is important for the jr × B
torque. Also, the heating location varies the fraction oftrapped supra-thermal electrons accelerated by ECH. The
heating location dependence on the jr × B torque would
appear through the fraction of trapped electrons. Additionally, including the finite parallel wavenumber, the resauthor’s e-mail: yamamoto y@p-grp.nucleng.kyoto-u.ac.jp

onance condition is shifted due to the Doppler eﬀect, andthe change of the heating source profile in velocity space
would appear.
In this study, we investigate the dependences of magnetic configurations and heating location on the toroidal
torque by ECH in LHD, applying the GNET code, which
can solve a linearized drift kinetic equation for δ f electrons
by ECH in the 5D phase space[8]. Also, we evaluate thetoroidal flow velocities driven by the ECH torques, solving
the diﬀusion equation of toroidal flow.

2

Simulation model

To investigate the electron velocity distribution and thetoroidal torque due to ECH, we apply the GNET code,
which can solve the drift kinetic equation in 5-D phasespace using the Monte Carlo method. We split the gyrophase averaged electron distribution function f into astationary Maxwellian part fMax and an perturbed part byECH δ f as f = fmax + δ f . The drift kinetic equation for δ f
is given by
∂δ f
∂δ f
+ (vd + v∥ ) ·+∂t∂r−C(δ f ) −∂δ f∂v
L(δ f ) = S ql ( fmax ) (1)v̇ ·where v∥ and vd are the velocity parallel to the magnetic
field and the drift velocity, respectively. Also, C(δ f ), L(δ f )and S ql ( fmax ) are the collision operator, the orbit loss term,
and the heating source term of ECH, respectively.
The ECH source term is described by the quasi-linear
diﬀusion theory. We consider only linear eﬀect S ql ( fmax )for simplicity. Then the source term S ql is expressed as
S ql ( fmax )= −∂ ql ∂ fmax
D∂vi i j ∂v j

(2)
-01

˜
ql

where Di j is the quasi-linear diﬀusion tensor and the radial
profile of the absorption power density is set as a Gaussian distribution. We consider that the right-handed electric field of the EC wave is dominant for X-mode.
ECH applies torques on the plasma through j × B andcollisions as below. Since the radial movements of trappedelectrons accelerated by ECH are faster than those of thermal electrons, ECH can drive the radial electron current je .
The net current in the steady state should vanish to maintain the quasi-neutrality, so the return current, jr (= − je ),must flow in the bulk plasma. Therefore, the bulk plasma
feels jr × B torque due to the return current. On the other
hand, the electrons drift toroidally due to the precession
motion. During the slowing down of the energetic electrons, they transfer their momenta to the bulk plasma due
to collisions. If we consider the heating source without initial momentum input, the collisional torque of the particles
passing in the co-direction should be equal to that of theparticles passing in the counter-direction. The trapped particles, however, have precession motion, which can contribute to the net collisional torque. The jr × B and collisional torques driven by the heating source without initialmomentum input should cancel in the completely symmetric configuration in the symmetry direction except for theeﬀect of the transient orbit width[11, 12]. Therefore, theconservation of angular momentum is satisfied and the total toroidal torque should vanish in the axisymmetric configuration. However, non-symmetric magnetic modes enhance the radial electron flux and break the cancellation of
the two torques. Since the jr × B torque is dominant in
LHD, the ECH torque is subject to the radial drift velocityof trapped electrons and the fraction of trapped electrons.

3 Results
Applying GNET code, we solve the linearized drift kinetic
equation for the supra-thermal electrons in 5D phase spaceand evaluate the deviation of the distribution function fromMaxwellian, δ f , in the inner shifted (Rax = 3.6m), thestandard (Rax = 3.75m) and the outward shifted (Rax =3.9m) configurations. We assume an ECH plasma with thecentral electron temperature T e0 = 4keV, the central ion
temperature T i0 = 1keV and the central electron density
ne0 = 1 × 1019 m−3 . Also, the toroidal magnetic field is
Bt = 1.375T.
We investigate the heating location dependence of thetotal torque, which includes the jr × B and collisional
torque. The total toroidal torque profiles with diﬀerent
heating location (r/a = 0.1, 0.3, 0.5, 0.7, and ripple top(θ = 180◦ , ϕ = 0◦ ) or ripple bottom (θ = 0◦ , ϕ = 18◦ ))and the integrated toroidal torque for each heating location are shown in Fig. 1. Here, the EC wave is assumed to
be injected perpendicularly. Ripple bottom heating makes
larger net toroidal torque than top heating because of thefraction of trapped electrons. Also, heating at outer minor

Fig. 1 The toroidal torque profiles with diﬀerent heating location in (a) the inward shifted configuration, (b) the standard configuration and (c) the outward shifted configuration. Figure (d) shows the dependences of the heating location and the configuration on the integrated ECH
torque.

radius makes larger toroidal torque than that by heating at
inner minor radius due to the strong magnetic ripple. The
integrated toroidal torques by ripple top heating vary moderately with diﬀerent minor radii of the heating location
than those by ripple bottom heating. The diﬀerence between the top and bottom heating is comparatively small in
the inward shifted configuration. The outward shifted configuration makes larger toroidal torque in the case of heating at inner minor radii because of the worse confinementof supra-thermal electrons. For the outward shifted configuration, however, the integrated toroidal torque of heating
at outer minor radii decrease because of the boundary.
We evaluate the toroidal flow V(r) driven by ECH
torque, solving the diﬀusion equation
(
)
∂V(r)
T (r)
∂V(r) 1 ∂=rD+
− µ(r)V(r), (3)
∂t
r ∂r
∂r
mi ni R
where T, D and µ, are the ECH torque, the diﬀusion coefficient and the viscosity coeﬃcient, respectively. The viscosity coeﬃcient is evaluated as[13]µ ≈ π1/2 ⟨(n̂ · ∇B)2 /B2 ⟩(R/M)(2eT i /mi ).

(4)

The viscosity coeﬃcients in the three configurations areshown in the Fig. 2 The obtained toroidal torque is thelargest in the outward shifted configuration, while the viscosity coeﬃcient near the axis is the largest in the outward
shifted configuration. The outward shifted configuration
has the strong driving and damping forces. On the other
hand, the inward shifted configuration has the weak driving
and damping forces. Their eﬀects conflict. We investigate-02˜Fig. 2 The neoclassical viscosity coeﬃcients in the three LHD
configurations.

Fig. 3 The toroidal flow velocities obtained by solving the diffusion equation in the three LHD configurations.

the toroidal flows in the three configurations. The steadytoroidal flows driven by heating near the axis r/a ∼ 0.1 areshown in the Fig. 3. As a result, toroidal flow in the standard configuration is the largest with the same input parameter (D = 1m2 /s, P = 1MW) because of its small viscosity and large toroidal torque. The second largest one is
obtained in the inward shifted configuration, which has thesmallest viscosity coeﬃcient around the axis in the three
configurations.

4 Conclusions
We have evaluated the ECH torque, which consists of thejr × B and collisional torques, in LHD plasmas applying
the GNET code to investigate the magnetic configurationand heating location dependences of the toroidal torque.
In a non-symmetric configuration, the jr × B torque canbe larger than the collisional torque and generate the net
torque. The jr × B torque is caused by the radial electron
current driven by ECH, which is subject to the radial drift
velocity of trapped electrons accelerated by ECH and thefraction of trapped electrons.
Heating at the ripple bottom makes more trapped electrons and thus generates larger toroidal torque than heating
at the ripple top. The supra-thermal electrons heated by
ECH need pitch angle scattering before they get trappedin the ripple top heating case. Heating at the outer minor
radius makes larger torque than heating at the inner minor
radius due to the strong magnetic ripple.

The inward shifted configuration has better confinement of supra-thermal electrons than that of the outward
shifted configuration. It indicates that the radial velocity ofsupra-thermal electrons is faster in the outward shifted configuration and generate a larger radial electron current by
ECH. Therefore, the outward shifted configuration makes
the largest toroidal torque when the ECH heating locationis near the axis. When the heating location is set at outerminor radii, the net torque decreases because of the boundary.
Finally, we have evaluated the toroidal flow with obtained toroidal torques, solving the diﬀusion equation. The
driving force, i.e., ECH torque, is the largest in the outward
shifted configuration, although the damping force, i.e., theneoclassical viscosity, is the largest in the outward shifted
configuration, too. As a result, the obtained toroidal flowis the largest in the standard configuration because of its
small viscosity and large toroidal torque. Also, the second
largest one is in the inward shifted configuration, which
has the smallest viscosity coeﬃcient in the three configurations.

Acknowledgments
The authors thank LHD Experiment Group for their
excellent support of this work. This work was supported by JSPS KAKENHI Grant Numbers JP18K03582.
This work is also supported by the National Institute for Fusion Science grant administrative budget(NIFS10ULHH021) and the NIFS Collaboration Research
program (NIFS19KNST150).
[1] M. Yoshida, et al., Nucl. Fusion 49, 115028 (2009)[2] S. T. A. Kumar, et al., Nucl. Fusion, 57, 036030 (2017).
[3] P. H. Diamond, et al., Nucl. Fusion 53, 104019 (2013).
[4] M. Yoshinuma, et al., Fusion Sci. Technol. 58, 375 (2010).
[5] K. Ida, S. Kado, and Y. Liang, Rev. Sci. Instrum. 71, 2360
(2000).
[6] Y. Yamamoto, et al., Nucl. Fusion (submitted)[7] Y. Yamamoto, et al., Plasma Fusion Res. 14, 3403105(2019).
[8] S. Murakami, et al., Nucl. Fusion 40, 693 (2000).
[9] M. Romé et al., Plasma Phys. Control. Fusion 39, 117(1997).
[10] S. Hasegawa, et al., Plasma Fusion Res., 8, 2403083 (2013).
[11] M. N. Rosembluth and F. L. Hinton, Nucl. Fusion, 36, 55(1996).
[12] A. Snicker, et al., Nucl. Fusion 55, 063023 (2015).
[13] K. Ida and N. Nakajima, Phys. Plasmas 4, 310 (1997).

-03

Coupling of EPED1 model to the integrated model
GOTRESS+ and development of JT-60SA operation scenario
M. Honda, N. Aiba, *H. Seto, E. Narita and N. Hayashi
National Institutes for Quantum and Radiological Science and Technology, Naka,
Ibaraki 311-0193 Japan*National Institutes for Quantum and Radiological Science and Technology, Rokkasho,
Aomori 039-3212 Japan

Extension of the integrated model GOTRESS+
A steady-state transport equation solver, GOTRESS [1,2], has been developed, which
benefits from global optimization techniques such as a genetic algorithm and the NelderMead method such that they directly find the solution where a transport flux matches anintegrated source. We have developed a novel integrated model GOTRESS+, mainly
consisting of the equilibrium and current profile solver ACCOME and the neutral beam
(NB) heating code OFMC, other than GOTRESS, for predictions of temperature profiles
consistent with an equilibrium and heating profiles. For the development of operation
scenarios in JT-60SA and ITER, etc., it is indispensable to predict the pedestal, which
accounts for a large proportion of plasma stored energy. In the past GOTRESS+
framework, the pedestal height was artificially given as a boundary condition of thetemperatures in GOTRESS. To predict the pedestal width and height in a reasonablefashion, the EPED1 model [3] has been incorporated into GOTRESS+. The EPED1model is based on the hypothesis that the pedestal is determined by two constraints: thepeeling-ballooning stability, and the empirical scaling of the pedestal height observed in
DIII-D, the latter of which relies on the kinetic ballooning turbulence theory. Our own
EPED1 model works together with ACCOME as an equilibrium solver and the MHD
stability code MARG2D [4]. The shapes of the temperature profiles in the edge barrier
region by EPED1 are characterized as the hyperbolic tangent ones. Multiple interface
programs have been newly developed to join a temperature profile by GOTRESS given
at discrete points and that by EPED1 given by the hyperbolic tangent function at theboundary of GOTRESS simulations. In short, GOTRESS+ is now extended to predict theplasma profiles entirely inside the plasma, provided that the density profiles are given.

Extension of the integrated model GOTRESS+We demonstrate an operation scenario development with GOTRESS+ for the ITERlike operation scenario in JT-60SA, what is called #4-1: 𝐵! = 2.28T, 𝐼" = 4.6MA and𝑃#$ = 24MW + 10MW. The 24MW NB heating is delivered by the positive-ion-source
NBI and the remaining by the negative-ion-source NBI. The CDBM model is used as thecore transport model. As shown in Fig. 1, the temperature profile and the safety factor
profile are converging as iteration continues. It can be seen that the temperature profilein the edge barrier region predicted by EPED1 varied for the first two iterations, but wasalmost unchanged after the third. Followed by the edge temperature profiles, the coretemperature profiles are determined subsequently. The plasma finally obtained virtually
meets the target values: the normalized beta of 2.71(~2.8), the energy confinement
improvement of 1.09(~1.1), where the values in the brackets correspond to the target
ones. The extended GOTRESS+ will be utilized for other scenarios.

References[1] M. Honda, Comput. Phys. Commun. 231 (2018) 94.
[2] M. Honda and E. Narita, Phys. Plasmas 26, 102307 (2019).
[3] P. B. Snyder et al., Phys. Plasmas 16, 056118 (2009).
[4] N. Aiba et al., Comput. Phys. Commun. 175, 269 (2006).

Fig. 1: Predicted profiles of the safety factor and the ion temperature at each iteration by
GOTRESS+. The right figure is the close-up of the ion temperature in the edge barrier region.

Extraction of turbulent structures and quantification of energy
transfer among structures by singular value decomposition
RIAM Kyushu Univ. A, RCPT Kyushu Univ.B, Warwick Univ.C, NIFSD,
IGSES Kyushu Univ.E, CCFEF, Shimane Univ.G
M. SasakiA, B,C, T. KobayashiD, Y. KawachiE,
R. O. DendyB,C,F, H. ArakawaG, S. InagakiA, B

Magnetically confined plasmas are the typical non-equilibrium open system so that thetransport of heat, particles and energy is self-consistently determined by turbulence and its
nonlinear structural formations. Dynamical behavior, such as limit cycle oscillation, is observed incases of zonal flows, toroidal intrinsic rotation, and transport barrier formation, etc [1]. In order tounderstand such phenomena, it is important to extract spatio-temporal structures with a minimal
number of degrees of freedom and to quantify the interaction among them.
In this study, we performed three-dimensional turbulence simulation to demonstrate how to
extract the spatio-temporal structures and to quantify the energy transfer among them by using the
singular value decomposition. We obtained the cases which corresponds to the Kelvin-Helmholtz
turbulence with a dynamical limit cycle oscillation of background plasmas which is accompanied
by flow pulsation and avalanche-like transport [2]. The abruptness of flow change contributes tothe sustainment of the limit cycle.
The singular value decomposition is applied tospatio-temporal data obtained in the simulation, andthe four typical structures are identified from the SVD
modes; the background flow deformation, zonal flow,a coherent fluctuation, and intermittent fluctuation.
The obtained zonal flow is found to be synchronizedwith the limit cycle, which determines the frequency ofthe zonal flow. Next, we formulated the energy of theextracted SVD modes based on Hasegawa-Wakatanimodel, and evaluated the energy transfer functions by
the simulation output. It was clearly shown that theFig. 1: Energy flow among theenergy transfer characteristics dynamically changes
with the limit cycle, and that the nonlinear coupling between the intermittent structure and thecoherent mode plays the significant roles for the flow pulsation [3]. The energy flow pattern issummarized in Fig. 1.
[1] K. Itoh, et. al., 8, 1102168 (2013).
[2] M. Sasaki, et. al., Phys. Plasmas, 26 042305 (2019).
[3] M. Sasaki, et. al., Plasma Phys. Control. Fusion, 63, 025004 (2020).

トカマクプラズマにおける不純物制御のための統合輸送シミュレーションスキームの開発
Development of Integrated Transport Simulation Scheme
for Impurity Control in Tokamak Plasmas
N. Kasuya 1,2, S. Mochinaga 2, A. Fukuyama 3, I. Shimomura 2, and M. Yagi 4
1
RIAM, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
2
IGSES, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
3
Kyoto University, Nishikyo-ku, Kyoto 615-8540, Japan
4
QST, Obuchi, Rokkasho, Aomori 039-3212, Japan
E-mail: kasuya@riam.kyushu-u.ac.jp
One of the most important problems in nuclear fusion reactor is control of impurities from thereactor wall, because high z materials are planned to be used for the diverter plates, which haslarge impact on the confinement of core plasmas. Some experiments show that the radio
frequency wave heating, as electron cyclotron heating (ECH) and ion cyclotron heating (ICH),
is effective for impurity control in the central region of tokamak [1]. Variation of the
background plasma profile is known to affect the neoclassical transport of impurities [2]. In
this research, integrated transport simulation code TASK [3] is developed to establish acontrol scenario of tungsten in the core plasmas.
Transport calculation module for impurities has been newly developed to include neoclassicaland turbulent transport processes and atomic processes [4]. The impurity transport module iscoupled with the transport module of main plasma components and wave heating module toevaluate response of both main and impurity components to ECH modulation (Fig. 1). Data
communication between the modules is implemented using data interface BPSD in TASK for
self-consistent calculation.
Transport analyses are carried out using the simulation scheme with plasma parameters on
middle size tokamak devices. Figure 2 shows time evolutions of the temperature of deuterium
ion and electron, density, power loss, and density of impurity ions. Neoclassical inward pinchin the region with large density gradient of the deuterium ions near the plasma edge induces
accumulation of impurities at the center of the plasma. The increase of the impurity density
leads increase of the line radiation losses as in Fig. 2(c), which causes decrease of thetemperature near the plasma center. Radiation collapse occurs in the larger impurity influxcase.
Response of neoclassical transport of impurities is focused on in the central region of
tokamaks, where the background pressure profile is not steep. It is demonstrated that radiation
collapse can be avoided with additional ECH. The case in Fig. 3 has the large impurity influx
enough to induce the radiation collapse, and the 2nd higher harmonic ECH at the center isapplied. Increase of the temperature and decrease of the density of the bulk ion by theadditional ECH reduce the impurity inward pinch to avoid the radiation collapse. In this way,
integrated transport simulation code TASK is developed to establish a control scenario oftungsten in the core plasmas. The simulation with fusion reaction is future task for next
generation tokamaks.

Fig. 1: Schematic of module configuration in TASK integrated simulation code used in thisresearch.

Fig. 2: Time evolutions of (a) temperature of deuterium ion and electron, (b) density, (c)power loss, (d) density of impurity ions. This is the case with ASDEX-U parameters.

Fig. 3: Time evolution of impurity ions when ECH is applied. (a) Temperature of deuterium
ion and electron, (b) density of deuterium ion and impurity ions, (c) input power, (d) loss
power are shown. The increase of the impurity density is halted with ECH injection from t =0.5 sec.

Acknowledgements
This work was partly supported by JSPS KAKENHI Grant Number (17H06089,
JP20K03905), by the collaboration program of NIFS (NIFS19KNST144) and RIAM of
Kyushu University.
References[1] R. Neu, et al., Plasma Phys. Control. Fusion 44 (2002) 811.
[2] C. Angioni, et al., Nucl. Fusion 54 (2014) 083028.
[3] https://bpsi.nucleng.kyoto-u.ac.jp/task/.
[4] I. Shimomura, et al., 17th BPSI meeting (Fukuoka 2019) 2-6.

Ion mass number dependence of resistive drift wave
turbulence
ISHIDA Masanobu1, KASUYA Naohiro1,2, IMAHASHI Yudai1, YAGI Masatoshi3
1

IGSES, Kyushu Univ., 2RIAM, Kyushu Univ., 3QST

1. Introduction
Since turbulence structures have great influence on plasma transport, understanding of theirformation mechanism is important for the control of fusion plasmas. Experiments in linear device
PANTA have shown that there exist several types of nonlinear structures generated by drift wave
instabilities, depending on parameters [1], such as streamers [2], zonal flows and solitary drift waves[3]. Numerical simulations also have shown structural bifurcation [4,5]. On the other hand, it is known
that the change of ion mass number leads turbulence suppression in fusion plasmas [6]. In this study,
we perform numerical simulations for the linear device using turbulence code Numerical Linear
Device (NLD) [7]. Difference of nonlinear saturation states of resistive drift wave instability isilluminated with change of discharge gas.

2. Simulation code
The target of this study is the linear device PANTA in Kyushu University, which produces acylindrical plasma with homogeneous magnetic field in the axial direction. Resistive drift wave
instability in PANTA plasmas is calculated using the three-field reduced fluid equations [7]. The set ofmodel equations is solved time evolutions of density N, electrostatic potential , electron velocity in themagnetic field line direction V for a plasma with a cylindrical shape (Fig. 1). Neutral particles exist
even in the center of this rather low temperature plasma, so the effect of neutrals is included by theion-neutral collision in terms. The boundary conditions are given as a periodic boundary in the axial
direction, and zero fluctuations at the plasma center and outer boundary in the radial direction. The mean
component is also calculated in the nonlinear simulation.
The simulation parameters are the followings; magnetic field B = 0.1T, electron temperature Te = 3eV,
plasma radius a = 7cm, device axial length  = 4m, electron-ion collision frequency νei/Ωic = 300,
electron-neutral collision frequency νen/Ωic = 10. The value of in is evaluated depending on the ion
species.

Fig. 1 Schematic of linear device plasma
3. Simulation results
First, a linear analysis of eigen modes is performed. The model equations are linearized, and the radialmode structure and the eigen values are calculated for each mode with azimuthal mode number m andaxial mode number n. Using the PANTA device parameters, the cases with discharge gas argon (Ar),neon (Ne), and helium (He) are analyzed. A fixed background density profile is used. Linear growthrate dependence for Ar, Ne and He with the same collision frequency in = 0.05 shows that there areunstable modes in the range of m < 10 and n < 3. In the He case, the linear growth rates in the wide(m, n) region increases (the region with growth rate > -0.01 is wider), but most of them are less thanzero. The maximum of the growth rate for He is also smaller than the cases for Ar and Ne. Actually,
the collision frequencies are different for different ion species. By considering the effective radius ofatoms, in of He is smaller than those of Ne and Ar. Figure 2 show the linear growth rate dependencewith consideration of the in value. The unstable modes exist in the wider range in the He case. Thisis due to the difference in the normalized distance. The distance is normalized by effective Larmor
radius 𝜌𝑀𝑇 ⁄𝑒𝐵 evaluated by using the electron temperature Te in this model, which depends

on the ion mass M. The longer normalized distance, even with the same plasma radius, enables modes
with larger mode number m and n to be unstable in the He case.

(a)

(b)

(c)

Fig. 2 m and n dependence of the linear growth rate for (a) Ar (in = 0.1), (b) Ne (in = 0.08) and (c)He (in = 0.01). The dotted line indicates the unstable boundary with zero growth rate.
Next, nonlinear calculations are performed for resistive drift wave turbulence. Figure 3 shows thetime evolution of mode energy in the case of Ar with in = 0.1. In the linear phase, modes with (m, n)= (2,1) and (3,1) become unstable in the linear phase, and are still dominant in the nonlinear phase.

Those higher harmonics are also excited, but the limited number of modes are dominant to form a
streamer. On the other hand, in the He case with in =0.01, the modes around (m, n) = (7,4) are excited(Fig. 4(a)), and their mode couplings give the nonlinear saturation (Fig. 4(b)). The energy of theunstable modes is transferred to the other modes, and a wide range of modes have the same level ofamplitudes, resulting a turbulent state (Fig. 4(c)). A zonal flow with (m, n) = (0,0) is also excited. The
linear analysis as in Fig. 2 predicts that the number of unstable modes is larger in the He case. The
difference in the linear mode spectra gives the difference in the nonlinear saturation states (zonal flowor streamer states) by changing the number of unstable modes.

Fig. 3 Time evolution of internal energy of each mode. The cases of Ar with in = 0.1 are shown. The(m, n) = (0, 0) component indicates the evolution of the background plasma density.

Fig. 4 Mode spectra in (a) linear, (b) developing from linear to nonlinear and (c) nonlinear phase of
the He case.
4. Conclusion
The ion mass dependence of the resistive drift wave instability, which is the driving source of
turbulence, is investigated for the linear device PANTA in the case of Argon, Neon and Helium
discharge gases. Nonlinear calculations are performed with the conditions obtained from the linear
analysis, and it is shown that the difference of the characteristic length (Larmor radius) causes thedifference of the spectra and the turbulent structure in the nonlinear phase. We will investigate thecharacteristics in mode spectra and turbulent states obtained in this analysis by comparing withexperiments.

Acknowledgement
This study is partially supported by the Grant-in-Aid for Scientific Research (JP17H06089,
JP20K03905) of JSPS, the collaboration program of NIFS (NIFS19KNST144) and of RIAM of
Kyushu University.

References[1] T, Kobayashi, et al., Plasma Fusion Res. 12, 1401019 (2017)[2] T. Yamada, et al., Nature Phys. 4, 721 (2008)[3] H. Arakawa, et al., Plasma Phys. Control. Fusion 53, 115009 (2011)[4] N, Kasuya, et al., Phys. Plasmas 15, 052302 (2008)[5] M, Sasaki, et al., Phys. Plasmas 22, 032315 (2015)[6] M, Nakata, el al., Phys. Rev. Lett., 118, 165002 (2017)[7] N, Kasuya, et al., J. Phys. Soc. Jpn. 76, 044501 (2007)

Eﬀects of Toroidal Field Ripple on Toroidal Torqueby Electron Cyclotron Heating in JT-60U Plasma
Hiroki ISHIGURO1 , Sadayoshi MURAKAMI1 , Yasuhiro YAMAMOTO1 , Mitsuru HONDA2
1. Department of Nuclear Engineering, Kyoto University
2. National Institutes for Quantum and Radiological Science and Technology

1Introduction
In recent years, spontaneous toroidal ﬂow due to electron cyclotron heating(ECH) has been observed in many tokamaks and helical devices such as
JT-60U [1], HSX, and LHD. ECH drives the radial electron current je byincreasing the radial motion of supra-thermal electrons. In the steady state,
the return current jr (=− je ) ﬂows to maintain the quasi-neutrality, andthe resulting torque (so-called J×B torque) drives the rotation of the bulk
plasma. When the supra-thermal electrons are generated by ECH, a toroidal
ﬂow is generated by the precession of the electrons. This toroidal ﬂow ofsupra-thermal electrons gives momentum to the bulk plasma by friction andthe torque is called the collisional torque. In the case of an axisymmetricﬁeld, the J×B torque and the collision torque cancel each other out [2].
In our previous studies, we have shown that ECH generates toroidal ﬂow in the plasma due to J×B
torque and collisional torque in helical plasmas [3]. Although the ideal tokamak is axisymmetric, themagnetic ﬁeld strength varies in toroidal direction called as

toroidal magnetic ripple

since the system

has a ﬁnite number of TF coils (JT-60U has 18 TF coils.). In this study, we evaluate the toroidal torquedue to ECH by varying the ripple strength using the GNET code. We assume the magnetic conﬁgurationof JT-60U and the toroidal ﬁeld ripple as
B = Baxisym. + B0,18 cos(18ϕ)

2where B0,18 = −δB0,0 (δ = 0 to 0.2 ％).

(1)

Simulation model
In this study, we use the GNET code [4]. GNET code solves the drift kinetic equation Eq. (2) forenergetic electrons in 5-D phase space (r, θ, ϕ, v∥ , v⊥ ) using the Monte Carlo technique.
∂δf
∂δf
∂δf
+ (vd + v∥ ) ·+ v̇ ·− C coll (δf ) = S ql (fMaxwell )
∂t
∂r
∂v

(2)

where δf is the displacement of the electron velocity distribution function from the Maxwell distribution
due to ECH. For the source term S ql , we use a quasi-linear diﬀusion term that takes into account the
spread in phase space. Since we assume second harmonic X-mode EC waves in this study, we use anapproximate expression for second harmonic X-mode for S ql .
(
)2 [
)2
) ]
(
ql (
3
D
k
v
v
γ
v
∥
∥
⊥
⊥
ij
2
S ql (fMaxwell ) = − 21−− v⊥
X exp(−X 2 )fMaxwell2−
vth vthvth∆c2
ωwhere X = [γ(1 − k∥ v∥ /ω) − 2Ωce /ω]/∆ and ∆ is a very small parameter (≪ 1).

1

(3)

3

ResultWe have simulated JT-60U with typical plasma parameters. (large radius R = 3.44m, small radiusa = 1.05m, magnetic ﬁeld in magnetic axis B = 3T). Fig.2 show typical trajectories of a trapped electron(initial electron energy is 10 keV). When there is no ripple (left. δ = 0.00 ％), the electron drifts along
the magnetic surface, but when there is a ripple (right. δ = 0.50 ％), the electron moves signiﬁcantly inthe radial direction. Fig.3 shows that as the ripple strength increases, the electron radial ﬂux increases.
Fig.4 shows the dependence of the toroidal torque on the ripple strength. The J×B torque increases as
the ripple strength increases. This is thought to be due to an increase in the number of particles trappedin the ripple and an increase in the electron radial ﬂux. On the other hand, the collisional torque isalmost unchanged. In the case of an axisymmetric magnetic ﬁeld conﬁguration (δ = 0.00 ％), the J×B
and collisional torques cancel and no net torque appears. The small amount of torque remaining isdue to the ﬁnite banana width eﬀect. On the other hand, for the non-axisymmetric ﬁeld conﬁguration(δ > 0.00 ％), a ﬁnite torque appears. The toroidal torque of JT-60U obtained with δ = 0.0 ∼ 0.20 ％ isabout 1/10 ∼ 1/100 of the typical torque obtained in a helical magnetic conﬁnement LHD.
Next, we compare the toroidal torque obtained from the simulation with the torque obtained fromthe theoretical model [5]. Considering that the particles trapped in the ripple drift directly in the radial
direction as shown in Fig.2(right), and assuming that the trapped supra-thermal electrons drift only forthe time they are trapped, the maximum value of the toroidal torque is given by[()]∫ ∆abs1
C2 n
ftr PECH
δ 3/2
ECH
Fmax ∼pr (l) = C1
1 − exp −.
dl∆abs 0⟨Ese ⟩nδ

(4)where pr (l) is the probability of existence of energetic electrons in the radial length l, ftr is the fractionof energetic electrons trapped in the ripple, PECH is the ECH absorption power, and ∆abs is the radialwidth of the ECH absorption distribution. In this case, δ = 0.0 ∼ 0.20 ％ around C1 ∼ 0.25, C2 ∼ 1.0.
The toroidal torques obtained by the simulation show a similar trend to the torques obtained by themodel.

[1] M. Yoshida, et al., Nucl. Fusion 49 (2009) 115028.
[2] M. N. Rosenbluth and F. L. Hinton, Nucl. Fusion 36 (1996) 55.
[3] Y. Yamamoto, et al., Plasma Fusion Res. 14 (2019) 3403105[4] S. Murakami, et al., Nucl. Fusion 40 (2000) 693.
[5] S. Murakami, et al., American Physical Society Division of Plasma Phyics (2019)2

Development of ECCD calculation module for
quasioptical ray tracing code
Kota Yanagihara
National Institute for Quantum and Radiological Science and Technology,
Naka, Ibaraki 311-0193 Japan
PARADE (PAraxial RAy DEscription) is newly developed quasioptical ray tracing code, which traces
propagation of arbitrary wave beams with refraction, diffraction, polarization, and dissipation in
inhomogeneous anisotropic media. In fusion sciences, ray tracing techniques have been used forprediction of electron cyclotron resonance heating and current drive (ECRH/CD). PARADE code canpredict ECRH power deposition profile from dissipation calculation with higher accuracy thanconventional ray tracing codes, but cannot ECCD for now. Thus, so as to predict the ECCD efficiency
profile in tokamak devices, development of ECCD calculation module optimized for PARADE has
been started. Adjoint technique, which can estimate the driven RF current with Spitzer function
instead of distribution function, has been widely used in many ECRH/CD codes, such as
EChamamatsu code. For sufficient result with reasonable calculation time, we also adopted the
adjoint technique for our ECCD module. In this paper, predictions performed by PARADE with
under-developing new ECCD module are preliminarily compared with the results of conventional
EChamamatsu code. Power deposition profiles and driven current profiles by PARADE with ECCD
module and EChamamatsu code in JT-60SA plasma (Fig.1) are shown in Figs.2. Toroidal injection
angle is fixed to 15°and poloidal injection angle is scanned. Width of each power deposition profiles
are similar between both codes. This mean that multi ray imitation of EChamamatsu code is good
approximation at least for defocused beam in large JT-60SA plasma with sufficiently low density andmagnetic field strength comparing to wave frequency, where diffraction is not crucial. Positions ofeach power deposition profiles are shifted to inner side for our ECCD code. We expect this might bedue to the influence of relativistic effect because PARADE code uses weakly relativistic dispersion
tensor for arbitrary wave vector while EChamamatsu code use weakly relativistic tensor for qasitransverse wave vector. The reason of this shift should be cleared as a future work. Driven current
profiles are in qualitative agreement but have some differences. Namely, current densities foroutside injections with 2nd heating and core injection with 3rd heating of our ECCD module arelower than those performed by EChamamatsu code. This might occur by difference between Spitzer
function solvers adopted in our ECCD module and EChamamatsu code. Our module analyticallysolves Spitzer function while EChamamatsu code numerically solves it. Influence of this difference
for accuracy and computational resources should be studied as a future work.
Baxis=0.00 (T), `=0.00
3
2.5

JT-60SA plasma
EC beams

2

-15

1.5

fce
rho(psi)
fuh=138GHz
frc=138GHz
fce1=138GHz
fce2=138GHz
fce3=138GHz

ly (m)

Z (m)

1
0.5
0

-0.5
Baxis=0.00 (T), `=0.00
-1

3

35
fce
rho(psi)
fuh=138GHz
frc=138GHz
fce1=138GHz
fce2=138GHz
fce3=138GHz

-1.5

2.5

-2

2

-2.5

1.5

-3

1

0

0.5

1

1.5

2

ly (m)

0.5

2.5 3 3.5
lR
(m)
x (m)

4

4.5

5

5.5

6

0

Fig. 1 Poloidal cross section of JT-60SA plasma flux surfaces, magnetic field strength contours, and-0.5

typical frequencies
for 138 GHz EC wave beams, whose poloidal injection angles are schematically-1indicated with-1.5
arrows.
-2-2.5

P (MW/m3)-3

Toroidal injection angle → 15
0 Poloidal
0.5 1 1.5injection
2 2.5 3 angle
3.5 4 →4.5scanning
5 5.5 6

0.3
0.25
0.2
0.15
0.1
0.05
0

lx (m)

25°

15°

5°

-5°

-15°

35°

0

0.1

0.2

0.3

0.4

Broad 3rd harmonic heating is appeared.
J (MA/m2)

Color : EChamamatsu
Black : PARADE-C

0.08
0.06 Due to 2nd. 25°
0.04
15°
0.02
35°
0
Due to 3rd.
-0.02
-0.04
0
0.1 0.2 0.3 0.4

0.5

0.6

0.7

0.8

P/m3
0.9

1

0.9

1

rho

5°

0.5

-15°

-5°

0.6

0.7

0.8

Figs. 2 Power deposition profiles and driven current profiles by PARADE with ECCD module (blacklines) and EChamamatsu code (colored lines) in JT-60SA plasma shown in Fig.1.

Data assimilation system based on integrated transport
simulation applying a reduced model of neutral beam injection
heating
Y. Morishita1, S. Murakami1, M. Yokoyama2,4, G. Ueno3,4,5
1Department of Nuclear Engineering, Kyoto University, Nishikyo, Kyoto, 615-8540, Japan
2 National Institute for Fusion Science, National Institutes of Natural Science, Toki 509-5292, Japan
3The Institute of Statistical Mathematics, Research Organization of Information and Systems, Tokyo 190-8562,

Japan
4The Graduate University for Advanced Studies, SOKENDAI, Toki 509-5292 and Tokyo 190-8562, Japan
5The Joint Support-Center for Data Science Research, Tokyo 190-0014, Japan

1. Introduction
Integrated simulation for fusion plasmas has various uncertainties in each of theemployed simulation models, specifically turbulent transport model. Because of this, thesimulation results also have uncertainties. To solve this problem, we are developing a data
assimilation system, ASTI[1], based on the integrated transport simulation code, TASK3D[2].
The purpose of this study is to develop a numerical system that can predict the
behavior of fusion plasmas with high accuracy employing the data assimilation. Dataassimilation can also improve simulation models involving turbulence models. As data
assimilation methods, We use the ensemble Kalman filter (EnKF) to predict the behavior of
fusion plasmas with high accuracy and the ensemble Kalman smoother (EnKS) to estimate
the models which can reproduce experimental data. To predict the behavior of fusion plasmas
in shorter time than real time, it is necessary to reduce the calculation cost of the NBI heating
power profile. In this paper, the results of data assimilation applying a reduced NBI heating
model are shown.

2. Data assimilation system
In TASK3D, one dimensional (1D) diffusive transport equation is solved. As aturbulent transport model for the LHD plasma, the gyro-Bohm model for electron and the
gyro-Bohm gradT model for ion are employed based on previous TASK3D simulations. In theprevious study, we employed HREYA and GNET-TD code to evaluate the NBI heating power.
However, this calculation takes several hours for a second in actual plasma. To shorten thecalculation time to evaluate the NBI heating power, we introduce a reduced model. At the firststep of the model, fast ion birth profiles by NBI heating for each beam are evaluated by theinterpolation of pre-calculation results of HFREYA and MCNBI code. This interpolation is
excited by the Gauss process regression, whose inputs are electron density and temperature
profiles. At the second step, the heating power profile are calculated using the simple
analytical solution of the Fokker-Planck equation. This calculation takes about a second for asecond in actual plasma. We adjust the model parameters (fast ion birth profile, criticalvelocity, and slowing down time) by the EnKF. Therefore, the state vector is composed of theelectron and ion temperature, density, numerical coefficients of turbulent transport models,
and the reduced NBI heating model parameters. Every state variable is defined on 60
computational grid points (radial direction, 1-D). The time series data of temperature anddensity profiles are assimilated into TASK3D.

3. Simulation results
We apply ASTI to a data set simulated by TASK3D using the NBI heating power
profiles calculated by HFREYA and GNET-TD base on the LHD shot 114053. Figure 1 shows
the assimilation results of Te and Ti by EnKF, that have been performed with 2000 ensemble
members for the cycle of assimilation, τDA = 80 msec. ASTI can accurately reproduce boththe electron and ion temperature in terms of those profiles and temporal variations through the
optimization of the employed turbulent transport model and the reduced NBI heating model.
Figure 2 shows the assimilation results of the NBI heating power. The solid lines
labeled “Prediction” are the results using the filtered estimate of the reduced NBI heating
model parameters, and the dashed lines labeled “TASK3D*” are the results using thesmoothed estimates. These results, especially using the smoothed estimates, agree well withthe results by HFREYA and GNET-TD (high cost model). It indicates that the reduced NBI
heating model supported by data assimilation can reproduce the time-space distribution of
NBI heating power calculated by HFREYA and GENT-TD at low calculation cost.

Fig. 1 : Prediction and filtered estimate of
temporal change of Te and Ti for ρ = 0.0 and ρ= 0.6. Here, ρ is the normalized minor radius.

Fig. 2 : Assimilation result of Pe and Pi for ρ =0.0 and ρ = 0.6.

4. Summary and future works
We are developing the data assimilation system, ASTI to predict and control the
behavior of fusion plasmas in real time. The fast prediction of LHD plasma temperature in ahigh accuracy has been realized by introducing the reduced NBI heating model into TASK3D
and optimizing the model parameters by ASTI. The calculation cost to evaluate the NBI
heating power has been reduced from several hours to about a second for a second in actual
plasma. To realize the real time prediction, we must speed up the predictive calculation bymaking the computational grid rougher and the time step longer, and introducing parallel
computing into TASK3D.

References[1] Y. Morishita, et al., Nuclear Fusion 60, 056001 (2020).
[1] S. Murakami, et al., Plasma Phys. Control. Fusion 57, 119601 (2015).

Evaluation of ballooning and kink instabilities
in PLATO tokamak
Y. Kuroda1, N. Kasuya 1,2, M. Sato 3, A. Fukuyama 4, M. Yagi 5, Y. Nagashima 1,2, A. Fujisawa 1,2
IGSES, Kyushu Univ. 1, RIAM，Kyushu Univ .2, NIFS 3, Kyoto Univ. 4, QST 5
1. Introduction
Tokamak type device PLATO is under construction in Kyushu university, which enables three dimensional
observations of plasma turbulence to improve understanding of plasma instabilities in torus devices [1]. The
flexible design of PLATO tokamak makes it possible to observe dynamics of 3-D turbulence structures.
Control of Magneto Hydro Dynamics (MHD) instability in tokamak devices is important for stable operationand improvement of confinement performance.
In this study, dynamics of MHD instabilities in PLATO are calculated. The case when ballooning and kinkmodes both exist is considered, and interaction between them is evaluated by changing the growth rate of theballooning mode.
2. Equilibrium in PLATO
The expected plasma parameters of PLATO tokamak are thefollowings; major radius R = 0.7m, minor radius a = 0.2m, plasma volume
VP = 0.9m3, electron temperature Te = 0.2keV, density n = 1.0×1019m3.
The plasma equilibria considering the experimental conditions areevaluated by using EQU module in integrated code TASK [2].
TASK/EQU calculates equilibria with the free boundary condition inexistence of external vertical coils [3]. Figure 1 shows the example of theobtained magnetic configuration. This calculation gives plasma current IP= 40kA, ellipticity κ = 1.6, triangularity δ = 0.3. Previous studies use
equilibrium by TASK/EQ, which is the module calculates the magnetic

Fig.1 Magnetic configuration

field in the last closed flux surface, for the initial magnetic fieldobtained by TASK/EQU.

distribution of MHD instability calculations [4]. In this study, aninterface using TASK/EQU equilibria is developed to introduce the effect of the external coil as anexperimental condition.
3. MHD simulation
MHD simulations are performed with equilibria obtained by TASK/EQU module as the initial condition.
Figure 2 shows the time evolution of the energy with different background pressure gradients. Larger central
beta value β0 gives larger linear growth rate of the ballooning mode, so the increase of the energy is enhanced.
The total pressure perturbation energy, which is integrated for the whole plasma volume, increases
monotonically in the linear phase, and show nonlinear saturation. Figure 3 shows the time evolution of thepressure profile. Here, the pressure is normalized by B02/µ0. The similar results are obtained for the ballooning
and kink mode excitation as in the previous research using TASK/EQ equilibria [4].

Fig.2 Time evolutions of the total

Fig. 3 Time evolution of the pressure profile in the poloidal

kinetic energy.

cross-section with 0 = 4.0%.

4. Characteristic mode structure
In this section the interaction between ballooning and kink modes is investigated by changing the plasmapressure. First, the case when the kink component is dominant for the whole evolution is presented. Here, the(m = 1, n = 1) component and (m = 14, n = 13) component correspond to the kink and dominant ballooning
mode, where m and n are the poloidal and toroidal mode numbers, respectively. Figure 4 shows the temporal
evolution of the mode structure with β0 = 0.5%. The (1 ,1) mode grows around ρ = 0.3 in the linear phase.
Figure 5 shows the pressure profile at that time. Internal kink mode [5] grows, and the destruction spreads
outward. After the spread, various modes are excited nonlinearly near the flux surface with safety factor q ~1 at t ~ 500. Then at t ~ 550 the (1, 1) mode grows again around ρ = 0.3. This repetitive excitation of the kinkmode is observed. The influence of the ballooning components is rather small in this case.
Next, the case when the ballooning component is dominant is presented. Here, the (m = 1, n = 1)
component and (m = 13, n = 12) component correspond to the kink and dominant ballooning mode,
respectively. Figure 6 shows the temporal evolution of the mode structure with β0 = 4.0%. In Fig. 6(a) the
kink mode structure can be seen around ρ = 0.3, as in the case with β0 = 0.5%. However, as in Fig. 7, whichshows the pressure profile at t = 200, the characteristic structure of the ballooning mode can be seen. The
linear growth rate of the ballooning mode is increase due to the larger pressure gradient. The (13, 12)component grows near the flux surface with q = 13/12 at t = 200 as in Fig. 6(b), and nonlinear couplings
among several modes induce the (1, 1) component at this region, connected to the kink mode excited at the
inner part of the plasma. Then, the mode amplitude becomes smaller after t = 270 as the pressure profile isflattened.

(b)Fig.4 Time evolutions of the mode structures with β0 = 0.5%. The (a)Fig.5 Pressure profile atballooning (m = 1, n = 1) and (b) kink (m = 14, n = 13) components aret = 400 with β0 = 0.5%.

shown.
0.6

(a)

(b)

0.022

0.02
0.4
0.018

0.016

0.2

0.014
0
0.012

0.01
-0.2
0.008

-0.4

0.006

0.004
-0.6
0.4

0.6

0.8

1

Fig.6 Time evolutions of the mode structures with β0 = 4.0%. The (a)Fig.7 Pressure profile atballooning (m = 1, n = 1) and (b) kink (m = 13, n = 12) components aret = 200 with β0 = 4.0%.

shown.

5. Conclusion and future plan
In this study, time evolutions in energy and perturbation amplitudes were evaluated by changing thepressure gradient when both ballooning and kink modes are excited. There is difference in the evolution of
the radial structure of the n = 1 component, and the magnitude of internal collapse is reduced. The observed
difference in the growth of the (m = 1, n = 1) component is due to the nonlinear interaction between theballooning and kink modes. In future we will investigate the unstable condition in PLATO by evaluating thedependence of MHD instabilities on the pressure gradient and magnetic shear, and compare it withexperiments.
References[1] A. Fujisawa, AIP Conf. Proc. 1993 (2018) 020011.
[2] http://bpsi.nucleng.kyoto-u.ac.jp/task/[3] K. Matsuoka, RIAM Reports, Kyushu Univ. 141 (2011) 51 (in Japanese).
[4] S. Tomimatsu, et al., Plasma Fusion Res. 15 (2020) 1403052.
[5] T. Tsunematsu, Fusion Res. 51 (1984) 101 (in Japanese).