Plasma and Fusion Research: Overview Articles

Volume 14, 3503059 (2019)

Recent Progress in the Numerical Simulation Reactor Research
Project∗)
Hideo SUGAMA1,2) and the Numerical Simulation Reactor Research Project Group1)1)2)National Institute for Fusion Science, National Institutes of Natural Sciences, Toki 509-5292, Japan
Department of Fusion Science, SOKENDAI (The Graduate University for Advanced Studies), Toki 509-5292, Japan(Received 27 December 2018 / Accepted 12 February 2019)
Fusion plasmas are complex systems which involve a variety of physical processes interacting with eachother across wide ranges of spatiotemporal scales. In the National Institute for Fusion Science (NIFS), we areutilizing the full capability of the supercomputer (Plasma Simulator) and propelling domestic and international
collaborations in order to conduct the Numerical Simulation Reactor Research Project (NSRP). Understanding
physical mechanisms of complex plasma phenomena for the systematization of fusion science, NSRP aims at
realization of the Numerical Helical Test Reactor, which is an integrated system of simulation codes to predict
behaviors of fusion plasmas over the whole machine range. In NSRP, eight task groups are organized to cover awide range of fusion simulation subjects: plasma fluid equilibrium stability, energetic-particle physics, integrated
transport simulation, neoclassical and turbulent transport simulation, peripheral plasma transport, plasma-wallinteraction, multi-hierarchy physics, and simulation science basis. Verification and validation researches are in
progress in these task groups collaborating with each other as well as with experimental and engineering groups.
Successful examples of validations of large-scale simulations of energetic particle driven instabilities and neoclassical and turbulent transport against experimental results from tokamaks and helical systems are highlighted. In
addition, recent achievements in advanced simulation studies on ion heating processes and plasma-wall interactions, as well as those in the application of Virtual-Reality (VR) technology to fusion engineering, are presented.
c 2019 The Japan Society of Plasma Science and Nuclear Fusion ResearchKeywords: numerical simulation, fusion plasma, helical system, tokamak, instability, transport, virtual reality
DOI: 10.1585/pfr.14.3503059

1. Introduction
Based on past research activities in large-scale numerical simulations, the Numerical Simulation Research
Project [1] was launched in 2010 to continue their tasks
and develop them in more systematic ways on the occasionof the re-organization of the National Institute for Fusion
Science (NIFS). The project was renamed as the Numerical
Simulation Research Reactor Project (NSRP) [2] in orderto accelerate the research activity towards the constructionof the Numerical Helical Test Reactor in 2014. Missions of
the NSRP are i) to systematize understandings of physicalmechanisms in fusion plasmas for making fusion science awell-established discipline and ii) to construct the Numerical Helical Test Reactor, which is an integrated systemof simulation codes to predict behaviors of fusion plasmas
over the whole machine range.
To fulfill the missions of the NSRP, the supercomputer
system Plasma Simulator has been made full use of. The
Plasma Simulator was replaced by Fujitsu PRIMEHPC
FX100 with the total peak performance of 2.62 Petaflops,
and the total main memory of 81TB as of 2015. The
author’s e-mail: sugama@nifs.ac.jp∗) This article is based on the plenary presentation at the 27th International Toki Conference (ITC27) & the 13th Asia Pacific Plasma Theory
Conference (APPTC2018).

Fig. 1 The Plasma Simulator, PRIMEHPC FX100.

snap shot of the Plasma Simulator, PRIMEHPC FX100, isshown in Fig. 1. The Plasma Simulator is widely utilized
not only for NSRP but also for numerical simulation studies done by many researchers in plasma physics and fusion
science. Figure 2 shows progress in the peak performanceof the Plasma Simulator and the number of submitted jobs
per month in past years.
The structure of NSRP is shown in Fig. 3. In 2017,
Prof. Ritoku Horiuchi, the Executive Director of NSRP,
was succeeded by Hideo Sugama. The Plasma Simulator Working Group makes management of the supercomputer system for eﬃciently performing simulation research. The steering committee consists of the Executive

3503059-1

c 2019 The Japan Society of Plasma
Science and Nuclear Fusion Research

Plasma and Fusion Research: Overview Articles

Volume 14, 3503059 (2019)

Fig. 2 Progress in the peak performance of the Plasma Simulator and the number of submitted jobs per month from2009 to 2017.

Fig. 3 The structure of the Numerical Simulation Research Reactor Project (NSRP).

Director, leaders of eight task groups, and directors of two
divisions, namely, the Fusion Theory and Simulation Research Division and the Fundamental Physics Simulation
Research Division in NIFS. All project members work insome of the task groups and participate in the NSRP Council to discuss research plans carried out through collaborations with domestic and foreign institutions as well as collaborations with the Large Helical Device (LHD) Project[3] and the Fusion Engineering Research Project [4] in
NIFS.
The eight task groups are organized to cover a wide
range of fusion simulation subjects: plasma fluid equilibrium stability [5], energetic-particle physics [6], integrated transport simulation [7], neoclassical and turbulent
transport simulation [8], peripheral plasma transport [9],plasma-wall interaction [9], multi-hierarchy physics [10],and simulation science basis [10]. Leaders and numbersof members of the eight task groups are shown in Fig. 4.
Each task group has from about 10 to 20 members. Many
collaborators from other projects in NIFS and from otheruniversities and institutions have joined these task groups
to enhance collaboration and cooperation studies. Most of
the simulation staﬀ members of each task group also participate in other task groups.
In the rest of this Overview, highlights of recent
achievements accomplished in NSRP are presented. In

Fig. 4 The eight task groups in the Numerical Simulation Research Reactor Project (NSRP). Leaders and numbers ofmembers of the eight task groups are shown.

Sec. 2, recent results from kinetic and magnetohydrodynamic (MHD) hybrid simulations of energetic particle
driven instabilities are shown. Next, neoclassical and turbulent transport studies based on drift kinetic and gyrokinetic simulations are explained in Sec. 3, where developments in turbulent transport modeling and an integrated
transport code are described, as well. Results from some ofthe other research activities in our project, such as particle
simulations of ion heating processes in spherical tokamaks,
multi-hybrid simulations of fuzzy structures on a tungsten
surface, and application of Virtual-Reality (VR) technology to fusion engineering, are reviewed in Sec. 4. Finally,summary and future plans are given in Sec. 5.

2. Simulation Studies on Energetic
Particle Driven Instabilities
Recent results from simulation studies on energetic
particle (EP) driven instabilities are presented in this section. Alfvén eigenmodes (AEs) driven by alpha particles
are one of the major concerns regarding the future burning plasmas because they can enhance alpha particle transport, reduce alpha particle heating eﬃciency, and deteriorate fusion reactor performance. For simulation studies
on such EP driven instabilities, the MEGA code was developed by Todo [11]. Hybrid simulations for EPs andmagnetohydrodynamics (MHD) are done by the MEGA
code, in which EPs such as fast ions and alpha particles
are treated by gyrokinetic Particle-in-cell (PIC) simulation
and a bulk plasma by MHD simulation. The coupling between EPs and MHD is taken into account through the EP
current in the MHD momentum equation. In addition, processes of neutral beam injection (NBI), EP collisions, andtheir losses are included in MEGA code simulations.
The MEGA code has been extensively applied to
study EP driven instabilities in tokamaks such as ITER[12], DIII-D [13], JT-60U [14], EAST [15], KSTAR [16],and AUG [17], as well as in stellarator/heliotron plasmas
such as LHD [18] and Heliotron J [19]. Validations of
MEGA code simulations have been performed on fast ion
profile flattening and electron temperature fluctuations in

3503059-2

Plasma and Fusion Research: Overview Articles

Volume 14, 3503059 (2019)

Fig. 5 Radial profiles of the fast ion pressure (left) and temperature fluctuations (right) obtained from the MEGA simulations and the DIII-D experiments (reproduced from
Ref. [13] with the permission of the IAEA publishing section).

Fig. 7 The EGAMs in LHD reproduced by the MEGA code simulation (reproduced from Ref. [24] with the permissionof the American Physical Society). (a) The poloidal velocity frequency spectrum including all of the frequency
components. (b) The poloidal velocity time evolution including only 50 kHz (red) and 100 kHz (green).

Fig. 6 The EP-driven TAE mode structures obtained from thesimulations (reproduced from Ref. [13] with the permission of the IAEA publishing section).

DIII-D [13], Abrupt Large Events (ALEs) in JT-60U [14],and AE bursts in LHD [18]. Also, fast ion profile stiﬀness due to the resonance overlap of multiple Alfvén eigenmodes is verified by the MEGA simulation [20]. Goodagreements between results from the MEGA simulations
and those from the DIII-D experiment regarding radial profiles of fast ion pressure and temperature fluctuations areshown in Fig. 5. The EP-driven TAE mode structures obtained from the simulations are given in Fig. 6.
For another example, the MEGA code is used to simulate synchronization of multiple AEs and fast ion energy flux profile evolution observed in the TFTR experiment [21–23]. The simulation shows that the time intervalof the bursts is about 3 ms and the maximum amplitude of
the radial MHD fluid velocity normalized by the Alfvén
velocity is given by vr /vA ∼ 3 × 10−3 , which are close
to the results from the TFTR experiment. The simulationis performed for the input power 10 MW by NBI whichis the same input power as in the TFTR experiment [23].
The fast ion energy flux was not experimentally found although the maximum fast ion energy flux shown by thesimulation reaches 60 MW, which is six times larger than
the input power. This large transport of fast ions causedby the multiple AEs implies a terrible deterioration of the
NB heating eﬃciency as observed in the TFTR experiment [23]. Physical mechanisms of the synchronizationof multiple AEs are elucidated by analyzing the time evolution of structures of the fast ion distribution function in
the space of particle energy and the canonical toroidal mo-mentum as follows. Before the AE burst, multiple AEs
grow to low amplitudes and flatten the fast ion distribution function, which takes a stepwise form. When the stepwise distribution reaches a critical distribution, the further
beam injection leads to broadening of the locally flattened
regions and their overlap. Then, the overlap of locally flattened regions brings about synchronized sudden growth of
AEs and global transport of fast ions. The MEGA simulation also confirms the profile resiliency which implies that
almost the same fast ion pressure profile and distributionfunction are obtained for the beam powers of 5 MW and
10 MW.
Energetic particle driven geodesic acoustic modes(EGAMs) are experimentally observed in tokamaks andhelical systems. Simulations are successfully performed
using the MEGA code to investigate the EGAMs observedin the LHD experiments [24]. Figure 7 shows that the frequency chirping of the primary mode and the sudden excitation of the half-frequency secondary mode observed in
LHD [25] are reproduced by the hybrid simulation using
the realistic physical condition and the three-dimensional
equilibrium. It is found from Figs. 7 (a) and (b) that thesecondary mode which has a half frequency (50 kHz) of the
primary mode frequency (100 kHz) is excited at t  1.8 ms
by energetic particles satisfying the linear and nonlinear
resonance conditions, respectively, for the primary andsecondary modes [24]. Both EGAMs have global spatial
profiles which are consistent with the experimental observations.
The bulk ion heating during the EGAM activity in
LHD is also demonstrated for the first time by the MEGA
simulation [26, 27]. The extended MEGA code which includes both kinetic EPs and kinetic thermal ions is devel-3503059-3

Plasma and Fusion Research: Overview Articles

Volume 14, 3503059 (2019)

oped to analyze EGAM physics such as energy channeling.
It is found from the MEGA simulation that, when the EPs
lose energy, the ions obtain energy through an energy channel established by the EGAM activity. The heating power
found by the simulation is 3.4 kW/m3 , which is close to thevalue evaluated from the experiment.
Another interesting application of the extended
MEGA code is stability analysis of LHD bulk plasmas. For
high beta LHD plasmas which are unstable against the resistive ballooning modes, linear and nonlinear simulations
using the extended MEGA code are performed. Then, theion kinetic eﬀects are found to reduce the linear growth rateof resistive ballooning modes [28]. Comparing nonlinear
MHD and kinetic MHD simulation results shows that thenonlinear MEGA simulation result is more consistent withhigh beta LHD experiments than the MHD simulation inthat the decrease of the central pressure for the former caseis smaller than that for the latter one.

3. Simulation Studies on Neoclassicaland Turbulent Transport
To understand and predict plasma confinement performance accurately, both neoclassical and turbulent transport
processes need to be studied especially for helical systems.
In this section, recent results from simulation studies on
neoclassical and turbulent transport in NSRP using various drift kinetic and gyrokinetic simulation codes are described.
The gyrokinetic Vlasov (GKV) code was developed
for Eulerian gyrokinetic simulation of turbulent transport
using a local flux-tube spatial domain [29]. At first, model
magnetic geometries of large-aspect-ratio toroidal systems
were employed in the GKV code although it can now treat
accurate magnetic configurations under the MHD equilibrium conditions corresponding to real experiments in tokamaks and helical systems. Validation of the GKV code
simulation is successfully performed for a JT-60U L-mode
plasma [30]. It is confirmed that the nonlinear GKV simulation can reproduce the experimental values of the ion
and electron heat transport fluxes by varying temperature
gradients within 30% from those observed in the experiment. In addition, the simulation shows that, in the coreregion, the ion temperature gradient (ITG) mode is a dominant instability while, in the edge region, the trapped electron mode (TEM) is dominant.
The GKV code is also applied to evaluate turbulent
transport fluxes in the high ion temperature (T i ) LHD discharge (#88343) [31, 32]. The nonlinear GKV simulation
using the flux tube located at the normalized minor radiusρ = 0.65 gives the turbulent electron energy flux which isin good agreement with the experimental observation. It
is found from the simulation at ρ = 0.65 that the turbulent particle flux is directed radially inward [32]. It is also
shown that when reducing the ion temperature gradient at
ρ = 0.65 from the experimental value by 20% the turbu-

Fig. 8 Collisionality dependence of the radial heat flux s=i,e q s
2evaluated by the mixing-length diﬀusivity γ/k⊥ for the
linearTEM (L-TEM) and ITG (L-ITG) modes in the LHD
inward-shifted hydrogen (H), deuterium (D), and tritium(T) plasmas (reproduced from Ref. [34] with the permission of the American Physical Society). For qualitative
comparison, nonlinear TEM simulation results for the
H and D cases [NL-TEM(H) and NL-TEM(D)] are alsoplotted.

lent ion heat flux given from the simulation agrees with theexperimentally observed value [32]. In other words, thepredicted temperature gradient deviates from experimental
observation about 20%.
In 2017, LHD deuterium experiments started. Investigation of isotope eﬀects on plasma transport is one of theimportant subjects of the deuterium experiments. Prior to
the LHD deuterium experiments, predictive studies withgyrokinetic simulations for isotope eﬀects also started.
Gyrokinetic TEM and ITG turbulence simulations in helical plasmas with hydrogen isotope ions and real-mass kinetic electrons are done using the GKV code [33, 34]. It isfound from the simulation results that combined eﬀects ofthe collisional TEM stabilization by the isotope ions andthe associated increase in the impacts of the steady zonal
flows at the near-marginal linear stability lead to the significant transport reduction. Figure 8 shows collisionality
dependence of the radial heat flux s=i,e q s evaluated by
the mixing-length diﬀusivity γ/k⊥2 for the linear TEM (LTEM) and ITG (L-ITG) modes in the LHD inward-shiftedhydrogen (H), deuterium (D), and tritium (T) plasmas.
Nonlinear TEM simulation results for the H and D cases[NL-TEM (H) and NL-TEM (D)] are also plotted in Fig. 8.
We see that collisional eﬀects on stabilization of TEM andreduction of TEM turbulent transport are stronger for D
plasmas than for H plasmas. In Fig. 8, nonlinear TEM
simulation results for the heat fluxes in the LHD H and
D plasmas are shown only for the single collisionality because the nonlinear GKV simulations for the LHD configuration need huge computational costs. Note that theheat fluxes obtained from nonlinear TEM simulations for
H and D plasmas in the cyclone-base-case (CBC) like tokamak with diﬀerent collisionalities are given by Fig. 4 of

3503059-4

Plasma and Fusion Research: Overview Articles

Volume 14, 3503059 (2019)

Fig. 9 Spatial structures of TEM potential fluctuations in hydrogen (H) and deuterium (D) plasmas in LHD and the
CBC-like tokamak (reproduced from Ref. [34] with thepermission of the American Physical Society). Turbulence simulation results for the H and D plasmas in LHD
are shown in (a) and (b), respectively, while those for the
tokamak H and D plasmas are in (c) and (d), respectively.

Fig. 10 Radial profiles of the ion heat diﬀusivities (χi ’s) evaluated from the experimental power balance for the H and
D plasmas (reproduced from Ref. [38] with the permission of the IOP publishing). The ion heat diﬀusivities
normalized by the gyro-Bohm diﬀusivity χGB(s) (s = H,D) are plotted as well.

Ref. [34] in which it is clearly shown that the TEM turbulent transport decreases with increasing the collisionality
as predicted by the estimate using the mixing-length diffusivity. Thus, the reduction of the TEM turbulent transport due to the collisional stabilization eﬀect can be expected for the LHD configuration as well although nonlinear TEM simulations for LHD with diﬀerent collisionalities still need to be done to get more definite conclusions.
Spatial structures of TEM potential fluctuations in H and D
plasmas in LHD and the CBC like tokamak are shown in
Fig. 9, where regulation of turbulence due to zonal flows ismore evident for deuterium plasmas in both LHD and tokamak cases. High-electron-temperature (T e ) LHD plasma
experiments relevant for TEM turbulence were done in order to examine the prediction of the GKV simulations as
described above [35]. In fact, better high-T e LHD plasma
confinement is experimentally confirmed for D rather thanfor H, which implies qualitative agreement with the simulation results although further detailed comparisons between simulation and experimental results remain as futuretasks.
Experiments of high-T i isotope plasmas are also donein LHD [36]. In the high-T i LHD experiments, formations of ion transport barriers (ITBs), where the ion heattransport is significantly reduced, are observed in the coreregions of both H and D plasmas although the heat transport and the density fluctuation amplitude are found to besmaller in the D plasma than in the H plasma. Radial
profiles of the ion heat diﬀusivities (χi ’s) evaluated fromthe experimental power balance for the H and D plasmas

are shown in Fig. 10, where the ion heat diﬀusivities normalized by the gyro-Bohm diﬀusivity χGB(s) (s = H, D)
are plotted, as well. Here, χGB(s) is defined by χGB(s) ≡23/2 ρ2ts vts /Rax , where vts ≡ (T s /m s )1/2 is the thermal velocity, ρts ≡ vts /Ω s the gyroradius, Ω s ≡ e s Bax /(m s c) the gyrofrequency, and Rax (Bax ) the major radius (the magnetic
field) at the magnetic axis in the vacuum field. It is clearlyseen from Fig. 10 that both χi and χi /χGB(s) are smaller inthe D plasma than in the H plasma.
Using the density and temperature profiles obtained
from the high-T i LHD experiments, GKV simulations arecarried out for the H and D plasmas [37, 38]. The mixinglength diﬀusivities for the H and D plasmas are evaluatedby χ ML ≡ ky (γ/ky2 )Δky , where the growth rate γ is givenby the linear GKV simulation as a function of the poloidal
wavenumber ky and Δky denotes the non-dimensional minimum grid size in ky . Figure 11 shows radial profiles of thenormalized mixing-length diﬀusivities, which are qualitatively consistent with the experimental results in Fig. 10
in that the normalized heat diﬀusivity is reduced in the Dplasma. Nonlinear GKV simulations are also performed
to obtain the turbulent ion heat diﬀusivities for the H and
D plasmas [37]. The simulations reproduce well the ratioof the ion heat diﬀusivity in the D plasma to that in the H
plasma evaluated from the experimental results although
the absolute values of the heat diﬀusivities are overestimated by the simulations. It is also found from the nonlinear GKV simulations that the ratio of the generated zonal
flow energy to the total turbulence energy at the normalized minor radius ρ ≡ r/a = 0.5 in the D plasma is about3503059-5

Plasma and Fusion Research: Overview Articles

Volume 14, 3503059 (2019)

Fig. 11 Radial profiles of the normalized mixing-length diﬀusivities for the H and D plasmas (reproduced from Ref. [38]with the permission of the IOP publishing).

30% larger than that in the H plasma.
The GKV code is also applied to investigate turbulent transport in quasi-axisymmetric systems. The Chinese First Quasi-axisymmetric Stellarator (CFQS) is a joint
project of the Southwest Jiaotong University (SWJTU) and
NIFS [39, 40]. GKV simulations of linear ITG modes and
ITG turbulent transport in the CFQS are performed [41]and their results are compared with those in the equivalent tokamak, which has a circular cross section withthe same safety factor q, magnetic shear parameter ŝ =(r/q)(dq/dr), and inverse aspect ratio r/R as those of the
CFQS. It is found from the linear GKV simulations that
the CFQS has wider poloidal wavenumber regions where
unstable ITG modes have higher growth rates than those in
the equivalent tokamak. However, the nonlinear GKV simulations show that the ITG turbulent heat flux in the CFQS
is comparable with or less than the heat flux in the equivalent tokamak, and that the stronger zonal flows are generated in the CFQS. Thus, the quasi-axisymmetric configuration is found to enhance the zonal flow generation and theresultant ITG turbulent transport regulation. The adiabaticelectron approximation is employed in the present nonlinear simulations for the ITG turbulent transport in the
CFQS. Eﬀects of kinetic electrons and finite beta on the
turbulence are to be examined in the future studies.
Neoclassical transport processes generally make significant contributions to plasma confinement in helicalsystems more than in tokamaks because radial drift motions of particles trapped in helical magnetic ripples bringabout high particle and heat transport fluxes. To investigate transport fluxes, viscosities, radial electric fields, andbootstrap currents produced by neoclassical transport processes in three-dimensional toroidal configurations, global
δ f Monte Carlo simulation code, FORTEC-3D [42, 43],
has been developed and applied to many devices. In non-Fig. 12 Time evolutions of GAM oscillations (top) and residual
zonal flow levels as functions of the radial wavenumber (bottom) in the LHD configuration obtained from theGT5D simulations and the analytical theory (reproduced
from Ref. [51] with the permission of AIP Publishing).

axisymmetric systems, radial electric fields are determinedusing the ambipolarity condition for the neoclassical particle fluxes. Thus, FORTEC-3D is used to calculate radial electric field profiles in several helical systems suchas LHD [44], TJ-II [45], and W7-X [46], and reasonableagreements between the simulations and the experimental results are confirmed. Also, FORTEC-3D is applied
to predict the neoclassical heat flux in the helical demo reactor, FFHR-d1 [47], as well as the neoclassical toroidal
viscosity (NTV) in JT-60SA with magnetic perturbations
imposed [48, 49].
Developments in global codes aiming at gyrokinetic
simulations of turbulent transport processes in helical systems have also started in NSRP. Extension of the global gyrokinetic simulation code (GT5D) to helical geometry is inprogress. GT5D is a global full-f gyrokinetic Vlasov simulation code which was developed by Idomura originally fortokamak configurations [50]. Verification studies of the extended version of GT5D are now being done for a typical
LHD configuration [51, 52]. Time evolutions of GAM oscillations and residual zonal flow levels as functions of theradial wavenumber in the LHD configuration are shown in
Fig. 12, where reasonable agreements between results fromthe GT5D simulations and the analytical theory [53,54] areverified. Note that GT5D is a full-f code which in principle can simulate both neoclassical and turbulent transport
processes. It is verified in Fig. 13 that neoclassical particle fluxes and ambipolar radial radial electric fields calcu-

3503059-6

Plasma and Fusion Research: Overview Articles

Volume 14, 3503059 (2019)is theoretically derived [57] as
φ(∞) =

Fig. 13 Radial profiles of neoclassical particle fluxes (a) andambipolar radial radial electric fields (b) calculated
from GT5D and FORTEC-3D codes (reproduced from
Ref. [51] with the permission of AIP Publishing).

lated from GT5D are in good agreement with those from
FORTEC-3D.
Extension of another global gyrokinetic simulation
code (XGC) [55] to helical geometry has been advanced
in collaboration between NSRP and Princeton Plasma
Physics Laboratory (PPPL) [56]. The X-point included
Gyrokinetic Code (XGC) is a global particle-in-cell code
which was originally developed by Chang et al. [55] toinvestigate neoclassical and turbulent transport processes
in tokamak peripheral plasmas including regions where
nested magnetic surfaces do not exist. Benchmark testsof the XGC against the GT5D code and the analytical theory are successfully done for the collisionless damping of
GAM oscillations in the LHD configuration. GT5D and
XGC global simulations of linear ITG modes and ITG turbulent transport remain as future tasks and comparisons
between results from these global simulations and the local GKV simulations are also planned.
Global gyrokinetic simulations are expected to clarify how turbulent particle and heat fluxes are influenced by
eﬀects of the global profiles of background densities, temperatures, and radial electric fields, which are diﬃcult for
local gyrokinetic simulations to accurately evaluate. Effects of the background radial electric field profile is oneof the major concerns from the viewpoint of mechanisms
for improving plasma confinement. We note here that, forthe helical system with the background radial electric field
Er , the formula of the residual zonal flow potential φ(∞) inthe long time limit after the collisionless damping process

φ(0),1 + G p + Gt + M −2
p (G ht + G h )(1 + T e /T i )(1)where φ(0) are the zonal flow potential in the long time
limit and that at the initial time, respectively, the dimensionless geometrical factors G p , Gt , Ght , and Gh are defined in Ref. [57], and M p ≡ |(vE /r)(Rq/vti )| represents
the poloidal Mach number defined with the backgroundE × B drift velocity vE ≡ cEr /B0 , the ion thermal velocity
vti ≡ (T i /mi )1/2 , the safety factor q, and the minor (major)
radius r (R). The validity of Eq. (1) is verified by the linear
GKV simulation using multi-flux-tubes distributed over agiven flux surface [58]. Since the poloidal Mach number
M p is proportional to the product of the background radial
electric field and the square root of the ion mass, the heavier ion mass gives the larger Mach number and the higher
residual zonal flow level according to the formula shownin Eq. (1). Thus, the favorable eﬀect on the turbulent transport regulation due to the hydrogen isotopes with heavier
mass is expected, and it is one of the interesting subjects
for GT5D and XGC global simulations to focus on in the
future.
Gyrokinetic turbulent transport simulations describedabove are regarded as being based on the rigorous firstprinciple model although they generally require huge computer resources. However, for the purpose of application
to a large number of transport analyses for experimental
results and reactor designs, it is desirable to perform simulations using simpler transport models with lower computational costs than the direct gyrokinetic turbulence simulations. In NSRP, turbulent transport models, which canquickly reproduce results from nonlinear gyrokinetic simulations of turbulent transport in helical systems, have beenstudied. Based on the GKV turbulence simulation results,
the model formula for the ion heat diﬀusivity χi in the ITG
turbulence is derived [59] asA1 ( k γ̃k /k̃y2 )αχmodel
i
=
,
(2)

χGB
A2 + τ̃ZF /( k γ̃k /k̃y2 )1/2
iwhere A1 ≡ 1.8 × 101 and A2 ≡ 5.2 × 10−1 are the dimensionless coeﬃcients, k̃y ≡ ky ρti the normalized poloidal
wavenumber, γ̃k ≡ γk /(vti /R) the normalized wavenumberdependent ITG mode growth rate, and τ̃ZF ≡ τZF /(R/vti )the normalized linear ZF response time [see Ref. [59] fordetails].
The χi model in Eq. (2) can well predict the nonlinear GKV simulation results, and it is applied to the integrated transport code (TASK3D) [60] which is found tosuccessfully simulate the T i profile in the high-T i LHD experiment [61]. Extension of the χi model is also done toinclude eﬀects of trapped electrons on the ITG turbulent
transport in the LHD plasmas [62]. In addition, quasilinear models for the ion energy flux, the electron energy flux,

3503059-7

Plasma and Fusion Research: Overview Articles

Volume 14, 3503059 (2019)

and the particle flux are newly constructed for the case of
the ITG turbulence including eﬀects of trapped electrons,
for which reasonable agreements between the quasilinear
flux models and the nonlinear GK simulations are verified [63]. However, further extension of these models andtheir validation for the case of the TEM turbulent transport as shown in Fig. 8 remain as future tasks. These newmodels will be applied to transport simulations predicting
profiles of T i , T e , and ne (the electron density) in the LHD
plasmas.
At the end of this section, we note recent progress inthe development of the integrated transport analysis suite,TASK3D [60]. As for the LHD experiment analysis version (TASK3D-a), it has been further extended by improving the NBI module to be applicable to H/D/He mixture
plasmas, realizing loose coupling to the local neoclassical transport code (DKES/PENTA) [64, 65] with the momentum conservation correction technique [66], and implementing the modules for analyzing energetic particles.
These extensions are targeted to the LHD deuterium experiments and facilitate database creation for researches
on key issues such as energy transport and hydrogen isotope eﬀects. Exploiting such accumulated database by
TASK3D-a, thermal transport modeling are made throughstatistical analyses and neural network approaches. Dataassimilation methods are also employed for the integrated
predictive modeling in collaboration with Kyoto University. Modules from TASK3D-a are continuously utilized
for core plasma designs for the the helical demo reactor,
FFHR [47, 67].

4. Studies on Basic Physical Processes
and Simulation Science Basis
Technologies
In NSRP, simulation studies on basic physical processes in plasmas and researches on simulation science basis technologies are actively progressed as highlighted inthis section. By means of particle simulations which mimic
merging plasmas in a spherical tokamak (ST), ion heating processes during magnetic reconnection in the presence of a guide toroidal magnetic field are studied [68–72].
Simulation results in Fig. 14 show that a circle- or an arcshape ion velocity distribution is formed in the downstream, which implies that ions are eﬀectively heated. In
addition, the simulation results reproduce tendencies observed in ST experiments, for example, the dependence of
the ion temperature on the toroidal magnetic field.
Simulations of complex processes caused by plasmawall interactions are also carried out in NSRP. It was
found in the experimental device NAGDIS-II at Nagoya
University that the low energy helium plasma irradiationon the tungsten (W) surface produces fuzz nanostructures[73, 74]. In order to investigate mechanisms of the tungsten fuzz formation composed of the multi-scale and multiphysics processes, the multi-hybrid simulation method has

Fig. 14 Ion velocity distributions formed in the downstream fromthe reconnection region.

Fig. 15 Application of Virtual-Reality (VR) to the fusion engineering. The viewer grasps the component by his virtual
hand (top). The viewer watches the movement of component by the robot system (bottom).

been developed [75, 76]. In the multi-hybrid simulation,
the irradiation process of the helium (He) plasma onto the
W surface is solved by the binary collision approximation(BCA), the diﬀusion process of the He atoms in the W material is treated by the random-walk model based on kinetic Monte-Carlo (KMC), and the deformation process ofthe W material due to the pressure from the helium bubble
is simulated by the molecular dynamics (MD). The BCAMD-KMC simulation elucidates the process in which thetungsten fuzz is grown by the He bubble agglomeration.
As an example of research activities in the science basis technologies, application of the Virtual-Reality (VR) to
the fusion engineering is advanced in NSRP [77, 78]. The
VR technology is expected to play important roles in thedesign, construction, and maintenance of the fusion reactor. Also, the robotics system will perform regularly the
maintenance, such as component replacement. Then, it is3503059-8

Plasma and Fusion Research: Overview Articles

Volume 14, 3503059 (2019)not easy to determine how to move and rotate the components on the 2D monitor by using the CAD software because the information regarding depth is lost. We can seek
the eﬃcient assembly and replacement procedures by investigating the component movements including the robot
motion in the VR space. Such examples of the VR application are shown in Fig. 15.

5. Summary
In this article, recent progress in the Numerical Simulation Reactor Research Project (NSRP) is presented.
NSRP is promoting large-scale numerical simulation researches through close collaborations with experimental
and engineering groups aiming at realization of Numerical Helical Reactor, which is an integrated system of simulation codes to predict behaviors of fusion plasmas over
the whole machine range. To construct reliable simulation codes which can make accurate predictions, verification and validation (V&V) studies are actively carried out
based on comparisons among improved theoretical models, simulation results, and experimental observations intokamaks and helical devices.
Highlighted topics are successful accomplishments
recently made by the following studies: (i) MEGA code
simulations of the EP-driven Alfvén eigenmodes, fast-ion
transport, EGAMs, bulk-ion heating, and kinetic eﬀects on
the RBM, (ii) FORTEC-3D code simulations of neoclassical transport, ambipolar radial electric fields, neoclassical toroidal viscosities due to resonant magnetic perturbations, (iii) GKV code simulations of the ITG and TEM turbulent transport, zonal flows, and isotope eﬀects, (iv) the
benchmark tests of the global gyrokinetic codes (GT5D
and XGC) regarding neoclassical transport and zonal flow
damping in the LHD configuration, (v) modeling of turbulent transport by the turbulent ion heat diﬀusivity andby the quasilinear fluxes of particles and energies for ionsand electrons, (vi) development of the integrated transport
code (TASK3D) (vii) particle simulations of the eﬀective
ion heating through merging plasmas in a spherical tokamak, (viii) multi-hybrid (BCA-MD-KMC) simulations offuzz nanostructures on a tungsten surface irradiated by ahelium plasma, and (ix) application of the virtual-reality(VR) technology to the fusion engineering. Not all activities in NSRP are covered in this article due to the space
limitation although part of recent NSRP results which are
not shown here are presented in Refs. [79–92].
As future plans, NSRP promotes further extensionsof the various simulation codes, developments of the integrated code system, and their V&V studies based on LHD
deuterium experiments and other device experiments totackle key issues such as energetic particle confinement,
kinetic eﬀects on MHD stability, isotope eﬀects on turbulence, impurity transport in core and peripheral plasmas, and plasma-wall interactions. NSRP will thus makesignificant contributions to plasma confinement improve-

ment, optimization of configuration for high plasma performance, and reactor design through further advancing
domestic and international collaborations.

Acknowledgments
The author (HS) deeply thanks the International Program Committee of the 27th International Toki Conferenceon Plasma and Fusion Research and the 13th Asia PacificPlasma Theory Conference for giving him the opportunityof the plenary presentation which has lead to this Overview
article. Support from all collaborators in the National Institute for Fusion Science (NIFS), other institutes, and universities to the Numerical Simulation Reactor Research
Project (NSRP) is highly appreciated. Numerical simulations were carried out with the use of the Plasma Simulator
and the LHD Numerical Analysis System at NIFS.
[1] R. Horiuchi, J. Plasma Fusion Res. 90, 2 (2014) (inJapanese).
[2] R. Horiuchi, J. Plasma Fusion Res. 92, 785 (2016) (inJapanese).
[3] Y. Takeiri et al., Nucl. Fusion 57, 102023 (2017).
[4] A. Sagara et al., Nucl. Fusion 57, 086046 (2017).
[5] K. Ichiguchi et al., J. Plasma Fusion Res. 92, 787 (2016)(in Japanese).
[6] Y. Todo et al., J. Plasma Fusion Res. 92, 806 (2016) (inJapanese).
[7] M. Yokoyama, J. Plasma Fusion Res. 92, 814 (2016) (inJapanese).
[8] R. Kanno et al., J. Plasma Fusion Res. 92, 794 (2016) (inJapanese).
[9] Y. Suzuki et al., J. Plasma Fusion Res. 92, 810 (2016) (inJapanese).
[10] S. Ishiguro et al., J. Plasma Fusion Res. 92, 821 (2016) (inJapanese).
[11] Y. Todo, Phys. Plasmas 13, 082503 (2006).
[12] Y. Todo et al., Plasma Fusion Res. 9, 3403068 (2014).
[13] Y. Todo et al., Nucl. Fusion 55, 073020 (2015).
[14] A. Bierwage et al., Nature Comm. 9, 3281 (2018).
[15] Y. Hu et al., Phys. Plasmas 23, 022505 (2016).
[16] T.L. Rhee et al., 11th Asian Pacific Plasma Theory Conference (Jeju, Korea, 2014), P26.
[17] J. Gonzalez-Martin et al., 45th European Physical Society
Conference on Plasma Physics (Prague, Czech Republic,
2018), I5.013.
[18] Y. Todo et al., Phys. Plasmas 24, 081203 (2017).
[19] P. Adulsiriswad et al., 60th Annual Meeting of the
APS Division of Plasma Physics (Portland, USA, 2018),
PP11.00068.
[20] Y. Todo et al., Nucl. Fusion 56, 112008 (2016).
[21] Y. Todo, 27th IAEA Fusion Energy Conference, (Ahmedabad, India, 2018), TH/1-2.
[22] Y. Todo, 27th International Toki Conference on Plasma andFusion Research and 13th Asia Pacific Plasma Theory Conference (Toki, Japan, 2018), P1-38.
[23] K.L. Wong et al., Phys. Rev. Lett. 66, 1874 (1991).
[24] H. Wang et al., Phys. Rev. Lett. 120, 175001 (2018).
[25] T. Ido et al., Phys. Rev. Lett. 116, 015002 (2016).
[26] H. Wang et al., 27th IAEA Fusion Energy Conference,(Ahmedabad, India, 2018), TH/P2-11.
[27] H. Wang et al., 27th International Toki Conference on

3503059-9

Plasma and Fusion Research: Overview Articles

Volume 14, 3503059 (2019)

Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), O-3.
[28] M. Sato et al., 27th IAEA Fusion Energy Conference,(Ahmedabad, India, 2018), TH/P5-25.
[29] T.-H. Watanabe et al., Nucl. Fusion 46, 24 (2006).
[30] M. Nakata et al., Nucl. Fusion 58, 074008 (2016).
[31] M. Nunami et al., Phys. Plasmas 19, 042504 (2012).
[32] A. Ishizawa et al., Nucl. Fusion 55, 043024 (2015).
[33] M. Nakata et al., Plasma Phys. Control. Fusion 61, 014016(2016).
[34] M. Nakata et al., Phys. Rev. Lett. 118, 165002 (2017).
[35] K. Tanaka et al., 27th IAEA Fusion Energy Conference,(Ahmedabad, India, 2018), EX/P3-6.
[36] K. Nagaoka et al., 27th IAEA Fusion Energy Conference,(Ahmedabad, India, 2018), EX/5-1.
[37] M. Nakata et al., 45th European Physical Society Conference on Plasma Physics (Prague, Czech Republic, 2018),
I2.105.
[38] M. Nakata et al., Plasma Phys. Control. Fusion 61, 014016(2019).
[39] H. Liu et al., Plasma Fusion Res. 13, 3405067 (2018).
[40] M. Isobe et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P1-41.
[41] M. Nakata et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), O-1.
[42] S. Satake et al., Plasma Fusion Res. 3, S1062 (2008).
[43] S. Satake et al., 2nd Asia-Pacific Conference on Plasma
Physics (Kanazawa, Japan, 2018), MFP-41.
[44] N. A. Pablant et al., Plasma Phys. Control. Fusion 58,
045004 (2016).
[45] F. Sánchez et al., Nucl. Fusion 55, 104014 (2015).
[46] N. A. Pablant et al., Phys. Plasmas 25, 022508 (2018).
[47] T. Goto et al., Nucl. Fusion 57, 066011 (2017).
[48] S. Satake et al., Phys. Rev. Lett. 107, 055001 (2011).
[49] M. Honda et al., Nucl. Fusion 58, 112012 (2011).
[50] Y. Idomura et al., Comput. Phys. Commun. 179, 391(2008).
[51] S. Matsuoka et al., Phys. Plasmas 25, 022510 (2018).
[52] S. Matsuoka et al., 2nd Asia-Pacific Conference on Plasma
Physics (Kanazawa, Japan, 2018), O-1.
[53] H. Sugama et al., Phys. Rev. Lett. 94, 115001 (2005).
[54] H. Sugama et al., Phys. Plasmas 13, 012501 (2006).
[55] C.S. Chang et al., Nucl. Fusion 57, 116023 (2017).
[56] T. Moritaka et al., 27th IAEA Fusion Energy Conference,(Ahmedabad, India, 2018), TH/P5-5.
[57] H. Sugama et al., Phys. Plasmas 16, 056101 (2009).
[58] T.-H. Watanabe et al., Nucl. Fusion 51, 123003 (2011).
[59] M. Nunami et al., Phys. Plasmas 20, 092307 (2013).
[60] M. Sato et al., Plasma Fusion Res. 3, S1063 (2008).
[61] S. Toda et al., J. Phys.: Conf. Ser. 561, 012020 (2014).
[62] S. Toda et al., Plasma Fusion Res. 12, 1303035 (2017).
[63] S. Toda et al., 2nd Asia-Pacific Conference on Plasma
Physics (Kanazawa, Japan, 2018), B-I46.
[64] P. Hirshman et al., Phys. Fluids 29, 2951 (1986).
[65] D.A. Spong, Phys. Plasmas 12, 056114 (2005).
[66] H. Sugama et al., Phys. Plasmas 9, 4637 (2002).
[67] T. Goto et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), I-03.
[68] S. Usami et al., Phys. Plasmas 23, 092101 (2017).
[69] S. Usami et al., Plasma Fusion Res. 13, 3401025 (2018).

[70] S. Usami et al., 27th IAEA Fusion Energy Conference,(Ahmedabad, India, 2018), TH/P4-9.
[71] S. Usami et al., 2nd Asia-Pacific Conference on Plasma
Physics (Kanazawa, Japan, 2018), F-I15.
[72] S. Usami et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-66.
[73] S. Takamura et al., Plasma Fusion Res. 1, 051 (2006).
[74] S. Kajita et al., J. Nucl. Mater. 418, 152 (2011).
[75] A.M. Ito et al., 2nd Asia-Pacific Conference on Plasma
Physics (Kanazawa, Japan, 2018), MF-I43.
[76] A.M. Ito et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-50.
[77] H. Ohtani et al., 2nd Asia-Pacific Conference on Plasma
Physics (Kanazawa, Japan, 2018), B-I3.
[78] H. Ohtani et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-75.
[79] R. Kanno et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P1-15.
[80] G. Kawamura et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P1-22.
[81] S. Toda et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P1-26.
[82] S. Satake et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P1-30.
[83] S. Ishiguro et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P1-67.
[84] H. Miura et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-01.
[85] K. Ichiguchi et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-03.
[86] N. Mizuguchi et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-03.
[87] M. Toida et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-12.
[88] H. Hasegawa et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-25.
[89] Y. Suzuki et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-26.
[90] K. Fujita et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-37.
[91] A. Takayama et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-51.
[92] H. Nakamura et al., 27th International Toki Conference on
Plasma and Fusion Research and 13th Asia Pacific Plasma
Theory Conference (Toki, Japan, 2018), P2-59.

3503059-10