第 1９回
核燃焼プラズマ
統合コード研究会
19th Burning Plasma Simulation
Initiative (BPSI) Meeting
日時：2022 年 1 月 6 日(木)－7 日(金)場所：オンライン および九州大学筑紫キャンパス応用力学研究所 2 階大会議室(Ver. 1, 2022 年 2 月 18 日発行)

合同会合第 19 回核燃焼プラズマ統合コード研究会
19th Burning Plasma Simulation Initiative (BPSI) Meetingおよび
核融合エネルギーフォーラムサブクラスター会合
シミュレーションクラスター 炉心プラズマモデリングサブクラスター
プラズマ物理クラスター 定常運転・制御サブクラスター(Ver.1.1)
日時：2022 年 1 月 6 日(木)－7 日(金)場所：オンライン および 九州大学筑紫キャンパス

応用力学研究所 2 階大会議室

1 月 6 日(木)13:00 – 13:10 事務連絡 Business announcement
企画セッション「不純物輸送」 Topical session “Impurity transport”(座長：矢木)
13:10 – 13:20 講演 1-1
矢木 (量研)Introduction
13:20 – 13:45 講演 1-2
相羽 (量研)Current status and issues on developing operation scenarios of fusion DEMO
reactor
13:45 – 14:10 講演 1-3
藤田 (名大)Analysis of impurity transport with TOTAL code in tokamak DEMO
14:10 – 14:30 講演 1-4
持永 (九大)Impurity transport simulation using integrated code TASK in tokamak plasma
14:30 – 14:40 講演 1-5
林 (量研)
Plan for developing IMPACT: overall plan of the integrated code and its part for
impurity transport
14:40 – 14:50 講演 1-6
本多 (京大)Plan for developing TASK/TX
14:50 – 15:00 議論
Discussion
15:00 – 15:15 休憩(座長：村上)15:15 – 15:40 講演 2-1

Xu (名大)1

Reversal of heavy impurity pinch velocity in trapped electron mode turbulenceby decreasing the magnetic shear
15:40 – 16:05 講演 2-2
滝塚 (阪大)Double leap-frog method for large-time-step particle simulation to keep Larmor
radius small
16:05 – 16:30 講演 2-3
本多 (京大)Development of an integrated model GOTRESS+ and its application toassessment of JT-60SA operation scenarios including the pedestal
16:30 – 16:55 講演 2-4
杉山 (量研)Evaluation of heating power required for L-H transition in JA DEMO
16:55 – 17:20 講演 2-5
福山 (京大)Progress of wave analyses in toroidal plasmas using the integrated TASK
17:20 – 17:35 休憩(座長：糟谷)17:35 – 18:00 講演 2-6
佐々木 (日大)Study on interaction between micro-turbulence and meso/macro-scale
structures based on wave-kinetic formalism
18:00 – 18:20 講演 2-7
森下 (京大)Development of control system for fusion plasmas using data assimilation
18:20 – 18:35 講演 2-8
安齋 (京大)Introduction of finite orbital width effects into TASK/FP and comparison withneoclassical transport theory
18:35

散会

1 月 7 日(金)9:00 – 11:00 ポスター講演 Poster
コアタイム core time 9:00 – 10:00 Group A
10:00 – 11:00 Group B(開室 start 8:45、閉室 end 13:00)
11:00 – 11:10 休憩(座長：糟谷)11:10 – 11:35 講演 3 横山 (核融合研)
Information sharing on kick-off meeting of the IAEA-TM on AI application(Session Leader：村上)
11:35 – 12:00 議論 Discussion2

12:00 – 12:10 事務連絡 Business announcement
12:10 – 13:10 昼休み(核融合エネルギーフォーラムサブクラスターとの合同会合)
13:10 - 13:20 林（量研）
サブクラスター関係連絡事項（ITER の統合モデリングエキスパートグループの日
本委員交代について）
13:20 - 13:50 宮戸（量研）
IFERC 計算機シミュレーションセンターの現状報告
13:50 - 14:20 横山（核融合研）
ITER に於ける統合モデリング活動の報告
14:20 – 14:50 若月（量研）
ITPA 統合運転シナリオグループ活動報告
14:50 – 15:10 星野（慶應大）炉心プラズマ＋炉工学炉材料モデリングサブクラスター合同会合の報告
15:10 – 15:30 林（量研）
今後の予定の議論
15:30

散会

3

Poster number
19th BPSI meeting (2022.1)
Group A
PA-1 瀬⼾春樹 (量研) H. Seto
A BOUT++ extension for full annular torus tokamak edge simulations
PA-2 相澤カランザ正⾂ (⽇⼤) M. C. Aizawa
Prediction of turbulence temporal evolution in PANTA by long-short term memory
networkPA-3 ⿊⽥侑 (九⼤) Y. Kuroda
Analysis of three-dimensional structure of MHD instabilities in tokamak
PA-4 登⽥慎⼀郎 (核融合研) S. Toda
Prediction of temperature profiles in helical plasmas by integrated code coupled withgyrokinetic transport models
PA-5 糟⾕直宏 (九⼤) N. Kasuya
Numerical diagnostic of plasma fluctuations by using 3-D simulation data in torus plasmas
Group B
PB-1 梅崎⼤介 (九⼤) D. Umezaki
Large-angle scattering transport effects in edge plasmaPB-2 成⽥絵美 (量研) E. Narita
Development of semi-empirical particle and heat transport model and improvement in its
turbulent saturation rulePB-3 髙橋凌⼀ (⽇⼤) R. Takahashi
Spatio-temporal analysis of plasma turbulence in PANTA by using singular value
decomposition
PB-4 野中奨 (九⼤) S. Nonaka
Effect of impurities on improved confinement in tokamak discharges
PB-5 糟⾕直宏 (九⼤) N. Kasuya
Drift wave turbulent structure depending on ion mass number in cylindrical plasmas

4Impurity Transport Simulation Using Integrated Code TASKin Tokamak Plasma
S. Mochinaga1), N. Kasuya1,2), A. Fukuyama3), S. Nonaka1) and M. Yagi4)
1)

1.

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

Introduction

One of the important issues for realizing a fusion reactor is the control of impurities in the core plasma. In devices
that confine high-temperature plasmas, tungsten W can flow into the core from the divertor plates or the first wall,and significantly reduces the plasma performance [1]. Therefore, tungsten control using radio frequency (RF) heating
was demonstrated in tokamak plasma experimental devices [2].
In this study, the impurity transport simulation scheme in the integrated code TASK [3] is developed to establish
the impurity control scenario using RF heating. In this work, three extensions are made in TASK. Firstly, the transportmodule is linked to the equilibrium code to include magnetic configuration. Secondly, the scheme is extended to deal
with the transport of multiple kinds of impurities such as tungsten, argon, neon, and carbon. Since the impurity
transport calculation in TASK has been developed at first focusing on neoclassical transport, the turbulent transportcoefficients have been set as given parameters in it. Therefore, thirdly, turbulence models of mixed Bohm/gyro-Bohm
model (mBgB) [4] and Multi-Mode-Model (MMM95) [5] are applied. Simulation results with the three extensions
is described here.

2.

Transport simulation models

The TASK code is one of the integrated codes developed in Japan, which has a modular structure including many
codes such as equilibrium, transport, and wave analysis code. The multiple modules exchange data through the datainterface BPSD [6], which enables self-consistent simulation.
In the transport calculations of plasma main component using the TASK/TR code, the following particle and heatdiffusive one-dimensional transport equations are solved for plasma species s [7].
1 𝜕 𝑉𝛤
1 ∂
(1)
𝑛𝑉
𝑆 ,
𝑉
𝜕𝜌
𝑉 ∂𝑡
1 ∂ 𝑉𝑄1 ∂ 3(2)𝑛 𝑇 𝑉 5/3
𝑃 ,5/3 ∂𝑡 2𝑉∂ρ𝑉
Here, ns is the plasma density, Ts is the plasma temperature, Ss is the particle source, Ps is the heat source, ρ is thenormalized radius, and V’ is the radial derivative of plasma volume. The particle flux Γs and heat flux Qs are given as
follows;𝜕𝑛〈|𝛻𝜌| 〉𝐷(3)
𝛤 〈|𝛻𝜌|〉𝑉 𝑛
,
𝜕𝜌
and
3
𝜕𝑇 3
⟨|𝛻𝜌|⟩ 𝑛 𝑇 𝑉⟨|𝛻𝜌|2 ⟩𝑛 𝜒
(4)
𝑄
𝛤𝑇 .
2
𝜕𝜌 2
Here, ˂|∇ρ|˃, ˂|∇ρ|2˃ are the metrics, Vs is the particle convective velocity, Ds is the particle diffusion coefficient,
VEs is the heat convective velocity, and χs is the heat diffusion coefficient. The diffusion coefficients are given by thesum of neoclassical and turbulent diffusion coefficients, and the convective velocities consist of the neoclassicalcomponent only. For the neoclassical transport model, NCLASS [8] is applied. For the transport model of turbulence,
mBgB or MMM95 models are applied. The mBgB model can be expressed as follows [9].
(5)𝜒
8.0 10-5 𝜒
3.5 10-2 𝜒 ,(6)
1.6 10-4 𝜒
1.75 10-2 𝜒 ,𝜒𝜒 𝜒𝐷,
0.3 0.7𝜌
.
(7)𝜒𝜒
The impurity particle diffusion coefficient is DZ = Di. The Bohm term χB and gyro-Bohm term χgB are given by𝑇 𝜌 1.0
𝑎 𝑑𝑝 ⁄𝑑𝑟 𝑇 𝜌 0.8(8)
𝜌 𝑐 𝑞2,
𝜒𝑇 𝜌 1.0𝑝𝑑𝑇 ⁄𝑑𝑟
(9)
𝜌2 𝑐
,
𝜒
𝑇

where, ρs is the gyro radius, cs is the sound speed and pe is the electron pressure. The MMM95 model is a theorybased model represented by a linear combination of coefficients by Ion Temperature Gradient (ITG) and Trapped
Electron Modes (TEM) of the Weiland model, and Resistive Ballooning (RB) mode and Kinetic Ballooning (KB)
mode of Guzdar-Drake model. The diffusion coefficients of MMM95 model are given as follows;
0.80𝜒 , /1.00𝜒 ,
1.00𝜒 , ,
𝜒(10)
𝜒0.80𝜒 , /0.65𝜒 ,
1.00𝜒 , ,(11)𝐷
0.80𝜒 , /0.65𝜒 ,
1.00𝜒 , ,(12)𝐷
0.80𝜒 , /1.00𝜒 ,
1.00𝜒 , .
(13)
In the impurity transport calculations using the TASK/TI code, the following particle transport equations including
atomic processes are solved for k-valent ionized impurity ions.
1 ∂
1 ∂
(14)
𝑉𝑛
𝛤
𝛾 𝑛 𝑛
𝛾 𝑛 𝑛
𝛼 𝑛 𝑛
,
𝑉 ∂𝑡
𝑉 ∂ρ
in the case of k = 1,
1 ∂
1 ∂
(15)
𝑉 𝑛1
𝛤
𝛾1 𝑛1 𝑛1 𝛼2 𝑛 𝑛2 𝛼1 𝑛1 𝑆 ,𝑉 ∂𝑡
𝑉 ∂ρ 1
where, γk is the ionization rate, αk is the recombination rate and S0 is the impurity source. The impurity temperatures
are assumed to be equivalent to bulk ion temperature. The γk, αk and radiation loss power of impurity are evaluatedwith OPEN-ADAS [10] database. In this simulation, the impurity source is set for the 1+ charged impurity ion,assuming that the ions of density n1 have the thermal velocity vth at the plasma edge. For impurity transport in TASK,
the time evolution is solved by linking TR and TI codes. The plasma density and temperature profiles data aretransferred from the TR code to the TI code, while the radiation loss power and effective charge number data aretransferred from the TI code to the TR code.

3.

Impurity transport simulation considering equilibrium

In this section, the effect of introducing an equilibrium in impurity transport is described. In the fixed boundary
equilibrium, the major radius R = 3.0 m and the minor radius a = 1.0 m are kept constant and only the ellipticities arechanged from κ = 1.4 - 1.8. As shown in Fig. 1, the volume derivative increases and the metric |∇ρ| decreases withincreasing ellipticity.
Assuming that the impurity transport coefficient profiles are constant and the background plasma is given to be
fixed as a parabolic profile with the density at the magnetic axis ne0 = 0.5×1019m-3 and the temperature at the magnetic
axis Te0 = Ti0= 1.5 keV, carbon influx from the plasma edge ΓC = 1.0×1015 m-2s-1 gives the dependence of the carbon
density profiles on the ellipticity in Fig. 2. The impurity density decreases with increasing ellipticity because thedecrease in |∇ρ| results in a decrease in the effective diffusion. This shows the effect of equilibrium in impurity
transport.

Fig. 1 Radial profiles of (a) the radial derivative of plasma volume, and (b, c) the metric values.

Fig. 2 Ellipticity dependence of carbon density profiles.

4.

Transport of multiple impurity species

In this section, simulations for the ITER hybrid scenario [11] with several impurity particle species are described.
The plasma parameters of ITER tokamak are R = 6.0 m, a = 2.0 m, κ = 1.7, toroidal magnetic field Bt = 5.3 T, Ip = 15
MA. 30MW of neutral beam (NB) injection and 20MW of ion cyclotron resonance heating (ICRH) are applied as theexternal heating. The main ion and electron profiles at the steady state for the case when only W, with a concentrationof nW / nD = 10-6, is contained in the plasma are shown in Fig. 3. The plasma density profiles are fixed as in Fig. 3(a).
Figure 4 shows the density profiles calculated for Ar only, W only, and both Ar and W cases. The Ar concentration isabout 10-3 in each case. The turbulent transport coefficients of impurities in these calculations are fixed to be constantin parabola oneas. There is almost no difference in the Ar density profiles, but the W density is smaller in the caseincluding Ar. This is because outward particle pinch of W is slightly larger due to the decrease in the electron
temperature caused by the radiation loss of Ar.

Fig. 3 Radial profiles of (a) plasma density, (b) plasma temperature, and (c) heating power.

Fig. 4 Radial profiles of (a) Ar and (b) W density. The dashed lines show the density calculated with single
impurity species, and the solid lines show the density with multiple impurity species.

5.

W transport simulation using turbulence models

In this section, the results of W transport simulations are shown by introducing mBgB and MMM95 models as aturbulent transport model, applied to L-mode plasmas in ASDEX-U tokamak discharges [12]. The plasma parameters
are R = 1.65 m, a = 0.5 m, κ = 1.6, toroidal magnetic field Bt = 2.0 T, Ip = 0.8 MA. As shown in Fig. 5, the plasma
density and temperature profiles are fixed, referring to the experimental data in Ref. 12. W influx is set to be ΓW =5.0×1015 m-2s-1.
Figure 6 shows a comparison of the density profiles at steady state using the turbulent diffusion coefficients ofthe mBgB model, the MMM95 model, and the neoclassical equivalent model that assumes the turbulent diffusion is
two times larger than the neoclassical diffusion obtained by NCLASS. The turbulent diffusion coefficients in thosecases are shown in Fig. 7. The neoclassical particle pinch component here is inward for ρ = 0.9-1.0, and outward up
to around ρ = 0.3. In the case of the mBgB model, the impurity profile is flattened because the value of the diffusion
coefficient is evaluated to be larger in a wider region than the other models. In the case of the neoclassical equivalentmodel, the value is smaller, so the impurity profile dependents on the pinch component. In the case of the MMM95
model, the Resistive Ballooning mode component is larger near the plasma surface, and the ITG/TEM mode
component is larger near the plasma center ρ < 0.6. Then, the impurity profile is a little bit hollow.

Fig. 5 Radial profiles of plasma (a) density and (b) temperature.

Fig. 6 Radial profiles of W density.

Fig. 7 Radial profiles of the turbulent diffusion coefficients using (a) MMM95 model, (b) mBgB model and (c) the
neoclassical equivalent model.

6.

Conclusion

Three extensions have been made to the impurity transport simulation scheme in the TASK code: linking transportand equilibrium, enabling the calculation of multiple impurity species, and applying the mBgB and MMM95 models
to the impurity turbulent transport model. In the future, the other turbulence model, GLF23 [13], will be applied andimpurity control simulation by ECH will be carried out.

References[1] R. Neu et al., Fusion Eng. Des. 65, 367 (2003).
[2] R. Neu et al., Plasma Phys. Control. Fusion 44, 811 (2002).
[3] M. Honda et al., J. Plasma Fusion Res. SERIES 6, 160 (2004).
[4] M. Erba et al., Plasma Phys. Control. Fusion 39, 261 (1997).
[5] Bateman, Glenn, et al., Phys. Plasmas 5, 1793 (1998).
[6] A. Fukuyama, J. Plasma Fusion Res. 85, 597 (2009).
[7] S. Mochinaga et al., J. Plasma Fusion Res. 16, 403093 (2021).
[8] W. Houlberg, et al., Phys. Plasmas 4, 3230 (1997).
[9] M. Erba, et al., Nuclear Fusion, 38, 1013 (1998).
[10] ADAS website, http://www.adas.ac.uk/[11] Na Y.-S. et al., J. Korean Phys. Soc. 35 3788 (2008).
[12] J. Schirmer et al. Nuclear Fusion, 46, 780 (2006).
[13] R. E. Waltz et al., Phys. Plasmas 4, 2482 (1997).

Double Leap-Frog Method for Large-Time-Step Particle Simulation
to Keep Larmor Radius Small
T. Takizuka 1, K. Ibano 1, S. Togo 2
1
2

Graduate School of Engineering, Osaka University, Suita 565-0871, Japan
Plasma Research Center, University of Tsukuba, Tsukuba 305-8577, Japan

Modification of the leap-frog (LF) method to keep the Larmor radius ρ small is presented forthe particle simulation of a plasma even when ΩΔt >> 1 (Ω: cyclotron frequency, Δt: time step).
The conventional LF method makes ρ very large for ΩΔt >> 1, and such a numerical condition
has been avoided in general. If the LF method is applicable to the case of ΩΔt >> 1, new particle
simulation codes can be more easily developed for a wide area of plasma physics. By repeating
the LF steps doubly and adopting the averaged velocity to advance the particle position, the ρ size
is kept real independently of the ΩΔt value. Proper nature on the energy conservation, magnetic
moment conservation and drift-velocity realization is safely inherited from the LF method.

1. Introduction
Numerical simulation is a powerful tool to understand underlying physics in experimental
observations for various fields of plasma science. It is also an indispensable means to predict plasma
natures in future experiments, such as fusion reactor etc. Global plasma simulations are especially
important for fusion plasma study. MHD/fluid modeling has mainly been applied to them, but it does
not sufficiently describe kinetic effects essential for low-collisionality.
Kinetic particle model is then valuable even for global plasma simulations, though it takes a highcomputational cost. In first-principle particle simulation codes including cyclotron gyration effects, theleap-frog (LF) method has been widely used to integrate equations of motion of an individual
numerical particle [1]. The LF method, however, gives a smaller cyclotron frequency ΩLF < Ω = qB/m,
and a larger Larmor radius ρLF > ρ = v⊥/Ω (q: electric charge, m: mass, B: magnetic field, v⊥: speed
perpendicular to B). Especially in the case of ΩΔt >> 1, ΩLF is limited to ~1/Δt << Ω and ρLF is enlarged
to ~v⊥Δt >> ρ due to the finite time step Δt. Of course particle simulations with gyration motions havegenerally been carried out under the condition of ΩΔt << 1. In PARASOL code for edge plasma
simulation (ΩeΔt >> 1 > ΩiΔt) [2, 3], the ion motion including gyration is solved by the LF method,while the motion of electron guiding center is solved by the predictor-corrector method.
We are studying the ion heat conduction in the mirror field by using PIXY code [4, 5], to find proper
formula of the ion heat flux in fluid simulations with anisotropic ion pressure (AIP) model [6, 7]. The
PIXY, which treats various kinds of particles to study the plasma-wall interaction, e.g. vapor shielding,is applying only the LF method. This is because the same subroutine can be simply used for motionsof multi species. During test simulations of SOL-divertor plasma (B intersects divertor plate obliquely,
Bx/Bt = 1/4), we surprisingly found that PIXY realized plasma profiles including presheath similar to
those simulated by PARASOL, even for ΩeΔt >> 1 (Fig. 1).

Fig. 1 Density profiles in a SOL-divertor plasma simulated by PIXY for large ΩeΔt values (a) ~ 6 (Bt = 0.1 T),(b) ~ 30 (Bt = 0.5 T) and (c) ~ 100 (Bt = 2 T). Magnetic presheath was well realized for all cases.

Why does the LF method work well even for very large ΩΔt value? If the LF method is applicable tothe large-time-step particle simulation of ΩΔt >> 1, new particle simulation codes can be more easily
developed and simulation studies (longer time scale with the larger time step) can be promoted for awide area of plasma physics. This is our motivation to develop a new algorithm based on the LF
method, which is improved to keep the Larmor radius small for large-time-step simulations.
2. Leap-frog method
The LF algorithm for charged-particle motion is written as follows,m (v+Δt/2 – v–Δt/2) / Δt = q {E(X0) + ½ (v+Δt/2 + v–Δt/2) × B(X0)}(1)(X+Δt – X0) / Δt = v+Δt/2(2)

where the velocity v and position X are set differently from each other with half-time-step separation.
The electric field E and magnetic field B in Eq. (1) are chosen at the particle position X0. One of notable
merits of the LF method is to assure the energy conservation; v+2 − v–2 = (2qΔt/m) E ⋅ ½ (v+ + v–).
2.1 Cyclotron frequency and Larmor radiusConsidering v⊥,+Δt/2 = v⊥ (sin θ+, cos θ+, 0) with θ+ = ΩLFΔt/2 + θ0 and X⊥,+Δt = ρLF (– cos θ, sin θ, 0)with θ = ΩLFΔt + θ0 in a magnetic field B = (0, 0, B), Eqs. (1) and (2) are reduced toΩLFΔt = 2 tan−1(ΩΔt/2)(3)

ρLF = ρ / cos(ΩLFΔt/2)(4)

The simulated cyclotron frequency ΩLF becomes smaller than real Ω, and the Larmor radius ρLF islarger than real ρ (see Fig.2). For ΩΔt < 1, the relative errors are ΩLF/Ω − 1 = −(ΩΔt)2/12 and ρLF/ρ − 1= (ΩΔt)2/8. Note that the relation “ρLF = v⊥/ΩLF” does not hold. When ΩΔt → ∞, ΩLF and ρLF converge
to π/Δt and v⊥Δt/2, respectively.
In order to correct the gyro-phase delay, the Boris algorithm has been adopted [8]. Although ΩBoris is justthe real Ω, the Larmor radius is affected as ρBoris = ρ |(ΩΔt/2) / sin(ΩΔt/2)|. When ΩΔt < 1, the relative
error, ρBoris/ρ − 1 = (ΩΔt)2/24, is smaller than that of LF method. On the other hand, when ΩΔt > 1, theradius is varied bizarrely as seen in Fig. 2 (b); ρBoris = v⊥Δt/2 is similar to ρLF at ΩΔt = (2l + 1)π, but it
becomes infinitely large at ΩΔt = 2lπ (l: integer).
The enlarged Larmor radius is a fatal demerit of the LF method if applied to ΩΔt >> 1. Strangely, theelectron Larmor radius, ρe,LF ∼ veΔt/2, becomes larger than that of deuterium ion, ρi,LF ∼ vi/Ωi, when 1> ΩiΔt > 2vi/ve ~ 1/30 (ve,i :thermal speed of electron and ion) as shown in Fig. 2 (c). Resultantly, theclassical diffusion perpendicular to B, D⊥ ≈ ρe2/τe (τe: electron collision time), becomes abnormallylarger by (ΩeΔt)2. Reduction of the wave coupling due to k⊥ρe,LF > 1 >> k⊥ρe is also the problem.

Fig. 2 Time-step dependence of (a) cyclotron frequency and (b) Larmor radius for LF method and Boris
algorithm. (c) Electron ρe,LF becomes larger than ion ρi,LF when ΩiΔt > 1/30 (deuterium ion).

2.2 Magnetic moment and drift
We examine to which extent the LF method correctly simulates the charged particle motion inelectric and magnetic fields when ΩΔt >> 1. As described above, ΩLF and ρLF converge to π/Δt << Ωand v⊥Δt/2 >> ρ , respectively. On the other hand, (i) the energy conservation is assured, and (ii) themagnetic moment, µ = mv⊥2/2B, is kept constant within the relative error ~ O(ρLF/LB) or O(1/ΩLFτB).
Therefore, (iii) the mirror force parallel to B is realized, FM = –µ∇||B. Here the characteristic length LB
is of the spatial variation of B, and the characteristic time τB is of its temporal variation. This favorable
property is based on a moment M = mρLF|½ (v⊥+ + v⊥–)| calculated from the r.h.s. of Eq. (1) being fully
independent of ΩΔt, M = mρLFv⊥ cos(ΩLFΔt/2) = mρLFv⊥. As for the drift perpendicular to B, (iv) the
E×B drift, VE×B = (E×B)/B2, is correctly simulated, and (v) the polarization drift, Vpolar = (dE/dt)/BΩ,
as well. Although ρLF >> ρ for ΩΔt >> 1, (vi) the curvature-∇B drift, V∇B = (2v||2 + v⊥2)(B × ∇B)/2B2Ω,
can be simulated without worry (v||: speed parallel to B).
3. Double leap-frog method
In order to keep the simulated Larmor radius small independently of the ΩΔt value, we propose a
new scheme based on the LF method [9]. The usual LF steps are repeated doubly as described below.
We call this scheme, therefore, “double leap-frog (DLF)” method.
1st LF step (Virtual position X ʹ′ and X # )( i ) v+Δt/2 – v–Δt/2 = (q/m){E ʹ′ + ½ (v+Δt/2 + v–Δt/2) × B ʹ′} Δt(5)#(ii) X Δt − X ʹ′0 = v+Δt/2 Δt(6)#

2nd LF step (Tentative velocity v +3Δt/2 )(iii) v #+3Δt/2 – v+Δt/2 = (q/m){E # + ½ (v #+3Δt/2 + v+Δt/2) × B #} Δt(7)(iv) XΔt − X0 = w+Δt/2 Δtaveraged velocity w+Δt/2 = (α/2) (v #+3Δt/2 – v–Δt/2) + (1 − α) v+Δt/2
with
α = 1/ 2(1 + cos(ΩLFΔt/2))

(8)
(9)
(10)

Resetting of virtual position( v ) X ʹ′Δt − XΔt = (α/2) (v+Δt/2 – v #+3Δt/2) Δt(11)

Gyration orbits of X (solid line) and X ʹ′ (dashed line) are schematically shown in Fig. 3 for ΩΔt =2.5 (ΩLFΔt = 1.79 and ρLF = 1.60ρ). The virtual position X ʹ′ (or X #) moves along an enlarged gyration
orbit with ρLF. To keep the real Larmor radius ρ, we use an averaged velocity w+Δt/2 to advance X in
Eq. (8). A simple reduction of the perpendicular movement, dX ⊥ = v⊥,+Δt/2 Δt × cos(Ω LFΔt/2),
cannot be applied because the E×B drift is also reduced. The gyration frequency ΩLF is not corrected
as it is, and the phase relation is the same as that of LF method; X⊥ ~ exp(iΩLFt) and v⊥ ~ i exp(iΩLFt).

Fig. 3 Gyration orbit of X with correct
Larmor radius ρ [2nd step] for ΩΔt = 2.5.
Velocity v, averaged velocity w, virtual
orbit X ʹ′ (or X #) with larger ρLF = 1.60ρ
[1st step] and resetting of X ʹ′(Δt) are
schematically shown.

Electric and magnetic fields in Eqs. (5) and (7) are chosen at the virtual position; E ʹ′ - B ʹ′ at X ʹ′0 andE # - B # at X #Δt. By this setting, proper nature on the magnetic moment conservation and drift-velocity
realization is safely inherited from the LF method. The electrostatic component Es = –∇φ (φ: potential)
could be chosen at the real position X0. Since XΔt is still not determined before Eq. (7), how to
self-consistently treat E # and B # is the future problem.
As for the particle flux Γ = nU (n: density, U: flow velocity), the component parallel to B iscalculated at the real position, while the perpendicular flux is calculated at the virtual position to
include diamagnetic flow; Γ (X) = Σj {S(X − Xj) ½ (v||,+Δt/2 + v||,–Δt/2) + S(X − X ʹ′j) ½ (v⊥,+Δt/2 + v⊥,–Δt/2)},
where S is a shape function in a Particle-in-Cell simulation (j: particle tag). Note that the diamagnetic
flow cannot directly be obtained in the guiding-center system.
The initial value of virtual position X ʹ′(t = 0) is determined from X(t = 0) and v(t = −Δt/2) by utilizing
the resetting Eq. (11), X ʹ′(0) − X(0) = (α/2) {v(−Δt/2) − v #(Δt/2)} Δt, with the tentative-velocity Eq. (7),v #(Δt/2) − v(−Δt/2) = (q/m) [E(X(0)) + ½ {v #(Δt/2) + v(t−Δt/2)} × B(X(0))] Δt.
Tentative variables X # and v # do not continue to the next time-step calculation. Collisional change
in v+Δt/2 then can be given after the above DLF process by using a collision model, such as theMonte-Carlo binary collision model [10].
4. Summary and discussion
Double leap-frog (DLF) method to keep the Larmor radius small is presented for the particle
simulation of a plasma even when ΩΔt >> 1 [9]. By repeating the LF steps doubly and adopting the
averaged velocity to advance the particle position, the Larmor radius is kept real independently of theΩΔt value. Proper nature on the energy conservation, magnetic moment conservation and driftvelocity realization is safely inherited from the LF method.
The DLF method for long temporal-scale particle simulation becomes more powerful when coupledwith the ingenious model [11,12] for large spatial-scale simulation. The DLF algorithm will be tested
during our development of the PIXY code (ΩeΔt >> 1 > ΩiΔt) in the near future. The possibility ofusing this new modeling in kinetic simulations instead of using the gyrokinetic approach (ΩiΔt >> 1)[13] will be examined in future studies. A weak point of the GK modeling is the absence of thepolarization drift movement. This weak point can be solved by the DLF method.
References[1] J.P. Verboncoeur, Plasma Phys. Control. Fusion 47, A231 (2005).
[2] T. Takizuka, Plasma Phys. Control. Fusion 59, 034008 (2017).
[3] T. Takizuka, Plasma Sci. Technol. 13, 316 (2011).
[4] K. Ibano, S. Togo, T. Lang, et al. Contrib. Plasma Phys. 56, 705 (2016).
[5] K. Ibano, Y. Kikuchi, S. Togo, Y. Ueda, T. Takizuka, Nucl. Fusion 59, 076001 (2019).
[6] S. Togo, T. Takizuka, M. Nakamura, et al., J. Comput. Phys. 310, 109 (2016).
[7] S. Togo, T. Takizuka, D. Reiser, et al., Nucl. Fusion 59, 076041 (2019).
[8] J.P. Boris, Proc. 4th Conf. on Numerical Simulation of Plasmas, p.3 (1970).
[9] T. Takizuka, K. Ibano, S. Togo, Plasma Fusion Res. 16, 1203100 (2021).
[10] T. Takizuka, H. Abe, J. Comput. Phys. 25, 205 (1977).
[11] T. Takizuka, K. Ibano, M. Yagi, Plasma Fusion Res. 13, 1203088 (2018).
[12] T. Takizuka, K. Ibano, M. Yagi, Plasma Fusion Res. 14, 1203091 (2019).
[13] X. Garbet, Y. Idomura, L. Villard, T.H. Watanabe, Nucl. Fusion 50, 043002 (2010).

Development of an integrated model GOTRESS+ and its
application to assessment of JT-60SA operation scenarios
including the pedestal
M. Honda, N. Aiba, *H. Seto, E. Narita and N. Hayashi
Kyoto University, Nishikyo-ku, Kyoto 615-8530 Japan
National Institutes for Quantum Science and Technology, Naka, Ibaraki 311-0193 Japan*National Institutes for Quantum Science and Technology, Rokkasho, Aomori 039-3212 Japan

The integrated model GOTRESS+ with EPED1
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 GOTRESS, the equilibrium and current profile solver ACCOME and theneutral beam (NB) heating code OFMC for predictions of temperature profiles consistent
with an equilibrium and heating profiles. For the development of operation scenarios in
JT-60SA, ITER and a DEMO reactor, it is indispensable to predict the pedestal, which
accounts for a large proportion of plasma stored energy in an H-mode plasma. Recently
the EPED1 model [3] has been incorporated into GOTRESS+ for predictions consistent
between the core and the pedestal plasmas. The EPED1 model is based on the hypothesisthat the pedestal is determined by two constraints: the peeling-ballooning stability, andthe empirical scaling of the pedestal height observed in DIII-D, the latter of which relieson 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
flowchart of a GOTRESS+ execution can be seen in Fig. 1. The workflow is regulated by
the Python script with the aid of a job scheduler installed in the machine where
GOTRESS+ runs. GOTRESS+ is now extended to predict the plasma profiles entirely
inside the plasma, provided that the density profiles are given.

Application of GOTRESS+ to JT-60SA plasmas
With the CDBM transport model used, GOTRESS+ has been applied to thedevelopment of two operation scenarios in JT-60SA: an ITER-like inductive operation
scenario (#4-1) and a high 𝛽! fully non-inductively current driven operation scenario(#5-1), where 𝛽! is the normalized beta. We here demonstrate the latter one: 𝐵" =1.72T , 𝐼# = 2.3MA , 𝑃!$ = 16MW and 𝑃%& = 7MW . The setting of the plasma
parameters is targeted at developing the DEMO-relevant plasma regimes in a steady state
without the inductive current. In general, it is difficult to simultaneously achieve them ina self-consistent manner, as all the current must be sustained by the current drive of theauxiliary heating systems and the bootstrap current, and the MHD stability must also be
met. The boundary of GOTRESS is set at 𝜌 = 0.8 and outside of it the profiles are
determined by EPED1.
As shown in Fig. 2, the GOTRESS+ simulation revealed that the consistent solution
could be obtained, which satisfied all of the target parameters such as 𝛽! =4.33(> 4.3), the energy confinement improvement 𝐻'((*,,) of 1.61(> 1.3) and thebootstrap current fraction 𝑓$. ~0.68(~0.68)[5]. Of course, the plasma was MHD stable.
The double barrier structure of the given density profile makes it possible to attain thesolution. The negative magnetic shear in conjunction with the alpha stabilization
suppresses the turbulent transport in the core region and forms the ITBs, which bringabout the steep pressure gradients and the resultant high 𝑓$. .
Currently the implementation of alpha heating due to fusion reactions in
GOTRESS+ is in progress for the purpose of predicting DEMO plasmas.

References[1] M. Honda, Comput. Phys. Commun. 231, 94 (2018).
[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).
[5] M. Honda et al., Nucl. Fusion 61, 116029 (2021).

Fig. 1: Flowchart of a GOTRESS+ execution with EPED1.

Fig. 2: Final state of (a) the electron and ion temperature profiles, (b) the electron and ion heatsources and (c) the electron density and the safety factor profiles for #5-1.

Evaluation of heating power required for L–H transition in JA DEMO
S. Sugiyama, N. Aiba, N. Asakura, N. Hayashi, Y. Sakamoto,and Joint Special Design Team for Fusion DEMO
National Institutes for Quantum Science and Technology

1 Introduction
In the DEMO reactor, it is desirable to reduce the volume of the ports for the heating and current
drive components from the viewpoint of tritium breeding. It is necessary to clarify the role of theheating and current drive in each phase of a plasma operation and to evaluate the required heating
power. In JA DEMO [1], which is a design concept of a tokamak DEMO reactor, the total heating
power of 80–150 MW can be prepared from the study of tritium breeding. The H-mode operation
is essential to obtain high-performance burning plasmas; the L–H transition is a stage requiring highheating power and should be completed within the heating power injected in the steady burn phase.
The heating power required for the L–H transition can be evaluated by comparing the loss powerfrom the plasma with the threshold power. Since the density evolution during the transition and theback transition conditions are not fully understood, assumptions on these are necessary to evaluate
the required heating power. The higher heating power is required when the density rise time is shorter
because of the increase of the change rate of the plasma stored energy, dW/dt [2]. In DIII-D and
ASDEX Upgrade, the H-mode can be maintained even when the loss power drops to half of thethreshold power, while in JET, no or weak hysteresis has been observed. It has been pointed out thatdW/dt should not be subtracted from the loss power after the L–H transition [3]; however, it is safer
to consider the eﬀect of dW/dt and assume no hysteresis for evaluating the heating power required
for the L–H transition of the DEMO reactor.
To reduce the heat load on the divertor, the loss power should be controlled below the allowablevalue by impurity injection. It is necessary to consider the scenario such that the L–H transition andthe divertor heat load are compatible by adjusting heating power and the impurity amount.
We have investigated the dependence of the heating power required for the L–H transition on thedensity rise time and the control of the loss power by using Ar injection in JA DEMO.

2 Analysis model
An L–H transition model was introduced into the TOPICS code [4]. The threshold power, PTH , isobtained from Martin’s ﬁtting expression [5] and multiplied by 0.8 to consider the ion mass dependence.
The L–H transition is started when the net loss power across the separatrix, Psep , exceeds PTH . Here,
Psep is deﬁned as the absorbed heating power minus the radiation power and dW/dt. The temperature
proﬁle is solved with the prescribed density proﬁle evolution and the CDBM model. The H-mode
temperature proﬁle outside the normalized minor radius, ρ, of 0.85 is ﬁxed using the pedestal proﬁlemodel [6]. The pedestal proﬁle during the L–H transition is determined so that the L- and H-mode
proﬁles are smoothly connected taking a prescribed pedestal formation time, ∆tped . The impurity ion
density proﬁle in the core region is solved by giving the time evolution of the density at the separatrix.
The density proﬁle was assumed to be [4(1 − ρ2 )0.6 + 1] × 1019 m−3 in the L-mode and the ﬂat proﬁlewith the Greenwald density fraction of 0.85. Argon was considered as the impurity. The anomalous
diﬀusion coeﬃcient of 2.2 m2 /s and the pinch velocity of −0.28ρ m/s were assumed for argon ions.

3 Results and discussion
100

Figure 1 shows the time evolution of Psep for ∆tped = 5 s. Af-90

ter the plasma current reaches the ﬂat top, NB is injected from

Dt

ped

P

LH

TH

NBI

80

66 s and the L–H transition occurs at 70.5 s. The loss power

P

70

P [MW]

drops by the increase of dW/dt after the transition started; intP

this case, the back transition occurs if no hysteresis exists.

sep

P

abs

60

Pa

NBI

abs

50

d

W/dt

40

Figure 2 shows the ∆tped dependence of the minimum value

30

P

20

of Psep /PTH during the transition. The higher heating power

rad

10

allows the shorter transition time. If the density can be con-065

trolled so that the transition takes 18 seconds, the transition is

66

67

68

69

70

Fig. 1

72

73

74

75

76

Time evolution of Psep .

achieved without hysteresis with the NBI power of 80 MW.
Figure 3 shows the dependence of Psep /PTH during the tran-

1.00

sition on Ar density fraction at the separatrix, fAr . The loss

P

NBI

0.95

= 80 MW120 MW
0.90

(

sep

P /P

)

TH minpower decreases as the radiation power increases by impurity
injection. When the impurity amount is too much, Psep drops

71Time [s]below PTH . The loss power cannot be suppressed by increasing
the impurity amount in the case of the high NBI power. By

0.85

selecting an appropriate combination of impurity amount and0.80

NBI power, Psep /PTH can be controlled to be almost constant.

5

6

7

8

9

10

11

12

Dt

13

ped

4 Conclusion
We have investigated the heating power required for the L–H

Fig. 2

14

15

16

17

18

19

20

[s]

Dependence of the mini-

mum value of Psep /PTH on ∆tped .
1.5transition of JA DEMO, introducing the L–H transition modelP

NBI= 120 MW
f

Ar

1.4

into the TOPICS code. The required heating power varies de-= 0.5%

1.3

pending on the density rise time and the back transition confraction, it is possible to control Psep and achieve both L–H

TH
sep

f

Ar

0.9

P

tained the prospect that the L–H transition is possible within

NBI

f

Ar

= 80 MW= 0.5%

0.8

the heating power supposed in the current design activities of

0.7
0

2

JA DEMO. Our future work is to evaluate the speciﬁc heating
describe the density and temperature evolutions.

= 0.3%

1.1

1.0transition and reduction of the divertor heat load. We have ob-

power required for the L–H transition, improving the model to= 0.1%Ar

P /P

dition. By adjusting the heating power and impurity density

f
1.2

4

6

8

t

Fig. 3
on fAr .

[1] K. Tobita, et al., Fusion Eng. Des. 136 (2018) 1024.
[2] C. Kessel, et al., Nucl. Fusion 55 (2015) 063038.
[3] O. Mitarai, et al., Fusion Eng. Des. 55 (2001) 477.
[4] N. Hayashi, et al., Phys. Plasmas 17 (2010) 056112.
[5] Y. Martin, et al., J. Phys.: Conf. Ser. 123 (2008) 012033.
[6] J. Garcia, et al., Nucl. Fusion 54 (2014) 093010.

10

-t

LH

12

14

16

18

20

[s]

Dependence of Psep /PTH

Progress of wave analyses in toroidal plasmas
using the integrated code suite TASK
Atsushi Fukuyama
Professor Emeritus, Kyoto University, Kyoto, Japan
There are various applications of wave in plasmas for control of plasma by waves including
plasma production, heating, current drive and rotation generation, and for diagnostics of plasma
by waves including Bremsstrahlung and synchrotron radiation, wave cutoﬀ and interference,scattering of waves by ﬂuctuation, and excitation by energetic particles. Various modules forwave modeling and analysis have been implemented in the integrated code TASK:
Modules related to wave analysis in TASKDP
dielectric tensor (cold, Maxwellian, arbitrary f , relativistic)
WR geometrical optics (ray tracing, beam tracing)
WM full wave optics (Fourier expansion, diﬀerential)
WF full wave optics (ﬁnite element method, integro-diﬀerential)
FPvelocity distribution function (kinetic transport, kinetic wave)
1

DP

Recent progress in wave dispersion analysis using DP has been made in the implementation of
perpendicular cold beam contribution in the dielectric tensor. It has been applied to the analysis
ion cyclotron emission by NBI in LHD plasmas, and the dispersion relation and the growth rate
shown in Fig. 1 are consistent with PIC simulation by Dr. Toida.
2

WR

Recent ECCD analysis for ITER plasma using the full-relativistic numerical estimate of theanti-Hermite term of dielectric tensor has suggested that the down-shifted second-harmonic EC
resonance strongly absorb the EC waves in high-temperature central plasma. Therefore lowﬁeld-side horizontal injection is not eﬃcient for ECCD in the central plasma (see Fig. 2), andtop-launch injection is promising for ECCD in the central and inner region of reactor plasma.
BB = 5.953
RR = 8.481
RA = 2.574
RKAP = 1.816

Q0 = 1.000
QA = 3.000
RF = 1.700E+05
NPHII= 0.1736

PA(1) = 0.0005
PZ(1) = -1.0000
PN(1) = 0.9000
PNS(1) = 0.0300

PTPR(1)= 35.0000
PTPP(1)= 35.0000
PTS(1) = 1.0000
MODELP1=206.
11

6PA(2) = 2.0000
PZ(2) = 1.0000
PN(2) = 0.4000
PNS(2)= 0.0133

PTPR(2)= 35.0000
PTPP(2)= 35.0000
PTS(2) = 1.0000
MODELP2= 0.

ABS POWER

10
9
8
7
4

6
5
4
3
2

2

1
0
0.0
0.2
POWER FLUX

0.4

0.6

0.8

1.0

0.4

0.6

0.8

1.0

1.0
0
0.8

0.6
-2
0.4

0.2

-4
6

7

8

9

10

11

0.0
0.0

0.2

Fig. 2: Strong second-harmonic
absorption
in high T e (0) = 35 keV
Fig. 1: Dispersion relation obtained by
TASK/DP for the ion cyclotron emission plasma prevents eﬃcient ECCD due torelativistic downshift of ω = 2ωcedriven by NBI in LHD
resonance
1

3

WM

Improvement of mode-expansion full wave code TASK/WM for the analysis of ICRF heating in
LHD plasma is under way. More poloidal mode numbers are required for accurate descriptionof wave structure in helical plasmas.
4

WF

There are several versions of FEM full wave modules in TASK:
Modules for FEM full wave analysis in TASK
W1
1D hot O-X-B mode conversion of EC waves in tokamaks
WF2
2D cold Helicon wave excitation in cylindrical plasmas
WF2D
2D cold O-X mode conversion EC waves in LHD plasmas
LH, Slow wave of IC
WF3D
3D cold IC waves in tandem mirror conﬁguration
In order to describe the kinetic eﬀects of plasmas in the full wave analysis using the ﬁnite element method (FEM) requires the integro-diﬀerential analysis which is free from wave number.
We have formulated the integro-diﬀerential dielectric tensor for unmagnetized plasmas (laserplasma interaction), magnetized longitudinal (magnetic beach heating) and magnetic transverse(ﬁnite Larmor radius eﬀects).
Accurate analysis of O-X-B mode conversion of EC wave requires kinetic full wave analysis
because both evanescent region and Bernstein waves have to be taken into account. Figure 3
shows the wave structure of the EC wave in a slab plasma simulating tokamak plasma. The EC
wave is excited as an O-mode (E ﬁeld parallel to the external magnetic ﬁeld in z direction) by anwave guide antenna located near the lower corner on the low-ﬁeld low-density right-hand-side
wall. The O-mode with optimum injection angle is mode-converted to the X-mode near the
O-mode cutoﬀ. The X-mode is reﬂected to the shorter wave-length X-mode in the high-density
region, and it is reﬂected again at the upper-hybrid-resonance layer and converted to the veryshort wave-length Bernstein mode which is absorbed at the fundamental cyclotron resonanceor by collision damping. Though the O-X-B mode conversion eﬃciency is similar to the cold
O-X mode conversion eﬃciency, the mode structure the Bernstein wave was ﬁrst obtained by
the integro-diﬀerential full wave analysis.
O-modeImEz (x, z)

X-modeReEy (x, z)

Bernstein-mode
ReE x (x, z)

Fig. 3: Kinetic full wave analysis of O-X-B mode conversion of EX waves in a slab plasma2

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, Japa
5The Joint Support-Center for Data Science Research, Tokyo 190-0014, Japan

1. Introductio
To develop a numerical system that can control the behavior of fusion plasmas withhigh accuracy, we are developing a data assimilation system, ASTI[1,2], based on the
integrated transport simulation code, TASK3D[3]. Data assimilation can also improve
simulation models involving turbulence models. As data assimilation methods, We use theensemble Kalman lter (EnKF) to predict the behavior of fusion plasmas and the ensemble
Kalman smoother (EnKS) to estimate the models which can reproduce experimental data. We
have developed a new control method for fusion plasmas using data assimilation. In thispaper, the results of numerical experiments to control the LHD plasma using the proposed
method are shown

2. Data assimilation syste
ASTI is a data assimilation system which employs the integrated transport simulation
code, TASK3D as the system model. ASTI can assimilate the experimental time series data
into the integrated simulation to predict and analyze the behavior of fusion plasmas. We haveextended ASTI to a control system by implementing the control method.
The proposed control algorithm consists of two components. The first is a filter forcontrol, which estimates the control input u* to achieve the target state z and computes thepredicted distribution under the control. This filter finds the optimal control u* by
assimilating z to the predicted distribution and calculates the predicted distribution under thecontrol by assimilating u*. This filter is assumed to be executed at Δtz intervals.
The second is a filter for observation, which optimizes the state vector using theobservation data obtained sequentially. In a situation where the observations are obtained at
Δty intervals (Δty > Δtz), it is required to assimilate the observations separated in time,
because the time at which the observation are obtained and the time of the latest prediction innumerical space differ by Δty. This can be achieved by assimilating the observation at a timet, yt, to the joint distribution of the two time points, p(xt , xt+Δt y ).

3. Simulation results

n

 

m

n

fi

.

n

A numerical experiment is conducted to investigate the control performance bycontrolling the virtual LHD plasma (TASK3D). We consider controlling the ECH input power
to make the electron temperature follow the target state time series as shown in Figure 1. The
ECH input power is adjusted by ATSI in the range 0.5 - 5 MW in 0.5 MW increments(Δtz = 0.1 s). We take the electron temperature at the plasma center (ρ = 0) as the target state  	Development of control system for fusion plasmas using data
assimilatio

z and assume that the radial profiles of electron temperature and ion temperature can beobserved in the cycle of Δty = 0.5 s. The electron density at the center is assumed to be
1.0 × 1019 m−3 . The state vector consists of electron temperature, ion temperature, ECH
input power, and numerical factors for thermal turbulent diffusivity of electron and ion. This
numerical experiment is performed by 400 ensemble members. Different models of electron
thermal turbulent diffusivity in ASTI and virtual plasma are assumed. Gyro-Bohm model andthe model that averages gyro-Bohm model and Alcator model are assumed in ASTI and thevirtual plasma, respectively.
Figure 2 shows the simulation results of the electron temperature controlled by ASTI.
It can be seen that the electron temperature under the control (green curve) agrees well for thetarget state time series (red curve). We have confirmed that the numerical factor for thermal
turbulent diffusivity is optimized to bridge the differences between the models in ASTI andthe virtual plasma.

Fig. 1 Target state time series of electron temperature at the plasma center.

Time [s]Fig. 2 Simulation results of the electron temperature controlled by ASTI.

4. Summary
We are developing the data assimilation system, ASTI to analyze and control the
behavior of fusion plasmas in real time. We have developed the control algorithm using data
assimilation and implemented it in ASTI. The numerical experiment has been conducted todemonstrate the performance of the control method. ASTI is expected to be a robust system
for analysis, prediction, and control of fusion plasma that connects an actual plasma and aplasma in numerical space.

References[1] Y. Morishita, et al., Computer Physics Communications 274, 108287 (2022).
[2] Y. Morishita, et al., Nuclear Fusion 60, 056001 (2020).
[3] S. Murakami, et al., Plasma Phys. Control. Fusion 57, 119601 (2015).

Implementation of Finite Orbit Width Bounce-averaged Fokker-Planck Simulation CodeTASK/FP
Akiyoshi Anzai, Keigo Ota, Atsushi Fukuyama, Mitsuru Honda, Sadayoshi Murakami
Department of Nuclear Engineering, Kyoto University, Nishikyo, Kyoto, 615-8140, Japan
1

In terms of computational grids, the set I = (p, θ0 , ρ0 ) is used in
TASK/FP[2]. Here, p is the momentum, θ0 is the pitch angle at aminimum magnetic ﬁeld, and ρ0 is the minor radius at a minimumﬁeld. Using these COM, the set (ϵ, µ, Pξ ) is derived from I = (p, θ0 , ρ0 ).
Another merit using this COM space is that the relation between thesevariables and orbits are obvious. For each ρ0 point, the values of particle
momentum p and pitch angle θ0 determine one orbit.
In TASK/FP, the diﬀusion tensor D uu and the advection vector F u
at a local space u = (p, θ, ρ, ϕ) are averaged over particle orbits, whichhave zero orbital width. Here, p is the momentum, θ is the pitch angle
at the minimum ﬁeld, ρ is the minor radius, and ϕ is the poloidal angle.
The tensor and vector driven by Coulomb collision can be expressed as
Dpp Dpθ 0
Fpuuu
D =  Dθp Dθθ 0  ,F =  Fθ  .
(5)
0
0
0
0

Introduction

In a fusion DEMO reactor, it is necessary to accurately simulate theplasma behaviors in order to construct an eﬃcient start-up scenario ofburning plasma. Currently, fusion reactor simulation codes are being developed in many countries, and research is being conducted to improve
the accuracy of the simulation scheme and to reduce the computation
time. One way to simulate plasmas is to solve the magnetohydrodynamic equation, but this is not suﬃcient to describe the behavior ofenergetic particles. It is necessary to solve the kinetic equation describing the behavior of the distribution function of plasma species including
energetic particles (EPs).
In an axisymmetric magnetic conﬁnement system such as tokamak,bounce-averaged Fokker-Planck equation (FP) where the distributionfunction is averaged over the particle gyro motion, the bounce motion in a magnetic mirror, and the axisymmetric direction, has been
introduced[1]. The bounce-averaged FP code TASK/FP was developed
to analyze the behavior of the velocity distribution function in toroidal
plasmas, and has been employed to analyze the interaction of plasmas
with various waves and toloidal electric ﬁelds in the presence of various
collision models[2]． However, the conventional TASK/FP analysis uses
an approximation that neglects the orbit width of the particles. As aresult, the behavior of EPs with large orbit width cannot be describedaccurately.

2

Zero orbital width approximation causes that radial variations of thediﬀusion tensor and advection vector due to collisions, are not taken
into account. The radial variable is constant in the original FP code.
Therefore we cannnot evaluate the radial transport due to collision andthe radial broadening of the collisional energy transfer.
orbit in (r,z)

Original TASK/FP code2

In the original TASK/FP code, we use the bounce-average FP equation with zero orbital width and calculate time evolutions of the constants of motion(COM) in the midplane. The bounce-average FP
equation is the kinetic equation which is averaged over for three timescales, ﬁrstly gyro-phase, secondly bounce-phase, and thirdly toroidal
angle[1][3]. On the assumption that gyro-phase and bounce-phase areshorter than the time scale of other interest plasma variations, these average operations are applied. As a result, we can reduce 3 dimensions in
the kinetic equation, and obtain the following Fokker-Planck equationin the COM space.
[
]
∂
∂
∂
JI f (I, t) =JI Dij− F i f (I, t) + JI H(I)
(1)
∂t
∂Ii
∂Ij
⟨ i
⟩
∂I j
∂I
Dij =
· D uu ·
(2)
∂u
u
⟩
⟨ i
∂I
Fi =
· Fu
(3)
∂u

0

-2

5

6

7

8

(b) The ﬁnite-width orbit for atrapped particle(a) Zero width orbit

Fig. 1: The Particle orbits in the tokamak

In the following section, we introduce the FOW variables into
TASK/FP and consider some conditions to solve these problems.

3

Here f is distribution function of particles, J is jacobian, I = (I1 , I2 , I3 )is three components COM vector, Dij is the diﬀusion coeﬃcient, F i isthe advection coeﬃcient and H is the source term. D uu is the diﬀusion
tensor and F u is the advection vector in a local space u. The notation⟨⟩ denotes the bounce-average. The summation convention implied.
Hereafter, the bounce-average operation is to be taken over poloidal
motion and ⟨⟩ is
∫
1 τp
dt
(4)
⟨⟩ =
τp 0

Introduction of FOW variables

To introduce FOW eﬀects into TASK/FP, we changed variables(p, θ0 , ρ0 ) to the new COM space I = (p, θm , ρm ). Here, p is the momentum, θm is the pitch angle at the maximum poloidal magnetic ﬂux ψm ,
ρm is the minor radius at ψm of an orbit. This is the diﬀerence between
our study and the previous work by Petrov[4]. In the new COM space,
the transform matrix is given by.

 ∂p
∂p
∂p
∂p
∂I
 m
=  ∂θ
∂p∂u∂ρmwhere τp is the bounce time on the drift orbit.
As COM variable, we choose the set I = (ϵ, µ, Pξ ) in which ϵ is theenergy, µ is the magnetic moment, and Pξ is the canonical toloidal angular momentum. However, these variables are inconvenient to calculate
the time evolution of distribution functions because these variables do
not contain isolated conﬁguration space coordinates, and also the momentum p is mixed in Pξ . This causes the diﬃculty in computational
grids.

∂p

∂θ
∂θm
∂θ
∂ρm
∂θ

∂ρ
∂θm
∂ρ
∂ρm
∂ρ




(6)

To obtain ∂I/∂u and the jacobian J , we calculate these mathematicalquantities (ϵ, µ, Pξ ) for u = (p, θ, ρ, ϕ),which are written as√
ϵ = m2 c4 + |p|2 c2 − mc2
(7)
µ=

1

1 p2 sin2 θm
1 p2 sin2 θ=
2 mB(ρ, ϕ)
2 mBm

(8)

Pξ =

F (ρ, ϕ)Fmp cos θ − qψ(ρ) =p cos θm − qψm .
B(ρ, ϕ)Bm(9)where m is the mass of a particle, c is the speed of light, ψ(ρ) is themagnetic ﬂux at ρ, and q is the charge. F = RBϕ denotes poloidal
currents where R is the toloidal radius, and Bϕ is the toroidal ﬁeld.
These conserved quantities are wrote down.
p=p

(10)

sin2 θ
sin2 θm=
B(ρ, ϕ)
Bm

(11)

F (ρ)
Fm
p cos θ − qψ(ρ) =p cos θm − qψm
B(ρ, ϕ)Bm(12)

Using these relations, the transform matrix can be obtained. Hereafter,
the prime indicates the deivative with respect to ∂/∂ρm .
In the kinetic equation, like the Boltzmann equation, the variables ofthe distribution function are in the phase space (x, p). On the other
hand, for computational grids, we use the new COM I = (p, θm , ρm ),which conserves (ϵ, µ, Pξ ). Therefore, the jacobian should be expressedin those three variable spaces.
We transform the variables (x, p) to (z, z̃), where z = (ϵ, µ, Pξ ), z̃ =(ϕ, ξ, Θ) and Θ is the gyro phase.Then, we average over z̃ and transformz to the new COM space (p, θm , ρm ). Finally, the jacobian for thetransformation from (x, p) to (p, θm , ρm ) is expressed as
JI =44π 2 τp
p3 sin θm
√
· 2
p 2·
2
m |q| m Bm 1 + ( mc)[{p(F´m Bm − Fm B´m ) cos2 θm
2
Bm
}
]
1
− F´m Bm sin2 θm − q ψ´m cos θm .
2

Fig. 3: Orbit region at ψm = 0.63

surface includes the point where vZ = 0 and are fragile for collisional
perturbation. Figure3 shows the cross-section at ψm = 0.63.

5

Conditions on the orbit boundaries

As discussed in the previous section, it is necessary to set the properboundary condition on pinch orbits in order to conserve particle numbers. The COM variables may change discontinuously when the particle
cross the pinch orbit surface. Figure4 shows the change of the COM
variables in phase space and Figure5 show the change of the particle
orbit in real space. Figure4 and ﬁgure5 show the orbit variation due

(13)

ψm1

ψm2

Orbit classiﬁcations

When we use the zero orbital width approximation, we neglect theeﬀects of collisional variations of drift orbits. However, with the FOW
eﬀects, it is necessary to classify orbits and to evaluate the variationof them since the diﬀusion and advection coeﬃcients are obtained byintegrating over the bounce motion.
There are several previous papers on the orbit classiﬁcation[5][6]. According to Rome and Peng (1979)[5], the orbit region is classiﬁed into
4 regions with ψm by unique COM space surfaces. Figure2 shows the

P
B

A

C

Fig. 4: Particle ﬂux and orbital boundaries

to the collision. In ﬁg4, the each point A, B, C, P corresponds to theorbit A (co-passing), B (trapped), C (counter-passing), and P (pinch)
respectively. Small collision perturbation on the orbits may cause thetransition between orbit regions represented by arrow. It should be
noted the particle ﬂux from the trapped region (red) has to be equal to
the sum of the ﬂux to the co-passing region (green) and the one to thecounter-passing region (blue).

6

Comparison of particle diﬀusion estimates

We compare the diﬀusion coeﬃcient of the FOW FP code with thediﬀusion coeﬃcient of the neoclassical theory. In order to evaluate only
the diﬀusion coeﬃcient, we set the advection therm as zero and thesource term are neglected. Figure6 shows the radial diﬀusion coeﬃcients of the electron and the deuterium ion as a function of minor
radius. The radial diﬀusion coeﬃcients are evaluated by the ratio of
the particle ﬂux and the density gradient. For FOW cases, the radialdiﬀusion is caused by orbit change due to collisional velocity change.
For neoclassical estimates, the neoclassical radial diﬀusion is averaged
over the bounce orbit. The FOW results are the same order with theneoclassical estimates, and shows larger diﬀusion near the magnetic axis
and smaller near the peripheral. The origin of the discrepancy requires
further analysis of the FOW eﬀects.

Fig. 2: Orbit classiﬁcation with orbit surfaces

7

Summary and future works

The FOW eﬀect has successfully been introduced in TASK/FP. The
schemes to properly classify orbits and set boundary conditions havebeen implemented. Radial diﬀusion coeﬃcients of FOW became larger
than the conventional neoclassical estimates near the magnetic axis.

surfaces in the 3D COM space represented by (v/c, cos θm , ψm ). The
green surface separates the passing and forbidden regions, and the blue
one the passing and trapped regions. The stagnation orbit on the green2ਤ 5.4: v/c − θm ฏ໘Ͱͷ‫ي‬ಓྖҬ(a) ρm = 0.55

ਤ 5.4: v/c − θm ฏ໘Ͱͷ‫ي‬ಓྖҬ(b) ρm = 0.23

ਤ 5.4: v/c − θm ฏ໘Ͱͷ‫ي‬ಓྖҬ

TrappedPinch

ナ領域における新古典輸送論との比較
バナナ領域における新古典輸送論との比較

Counter-passing
計測によるLHDプラズマ中の高速イオン速度分布関数再構成手法の開発
A Co-passing

閉じ込めは重布)を実験的した, FIDA計

存するため,ー)が期待でorbit in (r,z)

鈴木 航介1, 村上 定義1, 大舘 暁2, 藤原 大2, 山口 裕之2
1京都大学大学院工学研究科, 2核融合科学研究所

LHDにおけるFIDA計測• 実験時期• Shot number 146995• 重水素NBIの重水素運転期間中• 磁場配位• Rmax=3.60m
Bt=2.75T(反転)gamma=1.254
Bq=100%• 加熱ビーム• NB1: 174keV 0.8MW• NB2: 150keV 1.1MW• NB3: 175keV 1.5MW• NB4: 57keV 5.9M

2

FIDA計測によるLHDプラズマ中の高速イオン速度分布関数再構成手法の開発
鈴木 航介1, 村上 定義1, 大舘 暁2, 藤原 大2, 山口 裕之2
1京都大学大学院工学研究科, 2核融合科学研究所• 重み関数

各位置・周波数ごとで速度空間のFIDA信号への寄与を表す.
LHDにおけるFIDA計測
ピッチとエネルギーを軸に持つ.
• 磁場閉じ込め型核融合の実現に向けて高速イオンの閉じ込めは重• 実験時期
要な課題のひとつである.
• Shot number 146995• プラズマ中の高速イオン情報(速度分布関数・空間分布)を実験的• 重水素NBIの重水素運転期間中に取得する必要がある.
• 磁場配位FIDA計• 近年, 高速イオンの荷電交換反応で発生する光を利用した,再構成• Rmax=3.60m測が注目されている.
-2
Bt=2.75T(反転)• FIDA信号は発生元高速イオンの位置・速度情報に依存するため,
gamma=1.254• ある位置 r を測定するとき, 速度分布関数・空間分布 f (E, p, r) とFIDA計測からの速度分布関数の再構成(トモグラフィー)が期待でBq=100%周波数のFIDA計測信号 s (r, λ) は重み関数 w (E, p, r, λ) を通じて記
きる.
• 加熱ビーム
述できる.
NB1: 174keV 0.8MW
ヘ•• トカマク型では速度分布関数の再構成が比較的成功しているが,1 ∞7 p, r, λ)f8 (E, p, r)dpdE
リカル型では三次元性により再構成が困難となっている.
• NB2: 150keV 1.1MW
s (r,5λ) = 6w (E,∫−1 ∫0• NB3: 175keV 1.5MW
LHD実験• 本研究ではGNETコードおよびFIDASIMコードを用いて,• NB4: 57keV 5.9M
• 上式を離散化して以下の線形方程式を得る.
時のFIDA計測情報から高速イオンの速度分布関数を再構成する手

ϵ = r/R背景00• 重み関数R0：大半径 a：小半径

各位置・周波数ごとで速度空間のFIDA信号への寄与を表す.
ピッチとエネルギーを軸に持つ.

ϵ = r/R0• 重水素において逆アスペクト比が大きく• 重水素において逆アスペクト比が大きく再構成なるほど，拡散係数の値が新古典輸送の理論値から外れるのは，大アスペクト近IDA信号似が破れるためと考えられる．
ਤpoints
5.5: ఺
P ʹ͓͚Δϐϯν‫ي‬ಓͱྡ઀͢Δ‫ي‬ಓ• MHD不安定性のない実験のために接線NBIはイオン源を１
Fig. 5: Orbit shapes at eachਤ 5.5:
఺ P ʹ͓͚Δϐϯν‫ي‬ಓͱྡ઀͢Δ‫ي‬ಓ
S = (s , s , . . . , s )

しているが, ヘいる.
て, LHD実験
再構成する手つだけ使用することで入射パワーを制限している.

39

1

2

FIDA信号
M

因果F = ( f1, f2, . . . , fN )なるほど，拡散係数の値が新古典輸送のw (E, p, r, λ)f (E, p, r)dpdE
理論値から外れるのは，大アスペクト近∫ ∫• 上式を離散化して以下の線形方程式を得る.
似が破れるためと考えられる．
S = (s , s , . . . , s )周波数のFIDA計測信号 s (r, λ) は重み関数 w (E, p, r, λ) を通じて記述できる.
s (r, λ) =• MHD不安定性のない実験のために接線NBIはイオン源を１つだけ使用することで入射パワーを制限している.

速度分布関数の再構成• プラズマ中心において，有限軌道幅FPと
推定

と呼ぶ.

高速イオン情報の計測

39

GNETコードの計算による速度分布関数を元に, FIDASIMコードはFIDA信号を計算できる.

子ビームとの荷電交換によって光(656.3nm)を放出する.

• この光は発信源の粒子が高速で運動しているため波長が偏移(ドップラーシフト)する.

• この波長偏移量を計測することで高エネルギー粒子の速度情報を(a) Radial diﬀusion coeﬃcients ofelectron(b) Radial diﬀusion coeﬃcients of• 線形モデルを仮定すると,+ 最小二乗法によりフィッティング直線の
0
0
+
+ HNB
→ Hfast+ HNB
deuteriumion Hfast
係数が決まる.
• 一般的に, 観測は積分操作を伴うため観測量Sは物理量Fよりも低次元となる.
• 低次元情報から高次元情報を復元しようとするため, 解の不定性が生じる.

Fig. 6: Diﬀusion coeﬃcient of FOW and neoclassical theory

There is a possibility that the dependence is due to the potato orbit
eﬀects. In the following study, we implement the NBI heating term
into the FOW TASK/FP and evaluate the FOW eﬀects.The eﬀort to
reduce computation time by parallel processing in the presence of FOW
is under way.

References[1] M.G.McCoy A.A.Mirin J.Killeen, G.D.Kerbel. Computational Methods for Kinetic Models of Magnetically Conﬁned Plasmas. Springer-Verlag, 1986.
[2] H.Nuga and A.Fukuyama. Progress in nuclear science and technology, 2:78–84,
2011.
[3] K.Kupfer. Proc.IAEA Technical Committee Meeting on Fast Wave Current
Drive in Reactor Scale Tokamaks, page 149, 1991.
[4] Y.V.Petrov and R.W.Harvey. Plasma Phys. Control. Fusion, 58:115001, 2016.
[5] J.A.Rome and Y.M.Peng. Nuclear Fusion, 19:1193–1205, 1979.
[6] M.Okamoto S.Satake, H.Sugama and M.Wakatani. Plasma Phys. Control. Fusion, 77:573–585, 2001.

3

1

2

M

F = ( f1, f2, . . . , fN )• プラズマ中心において，有限軌道幅FPとと呼ぶ.

比べ新古典輸送理論の径方向拡散係数が線形回帰問題としての再構成
小さくなるのは，通常の新古典輸送理論
線形回帰問題に帰着される.
ではプラズマ中心に存在する特殊な軌道による輸送の効果が含まれていないことが原因の一つと考えられる．

• S = W F を F について解くとはWを入力としてSを出力とする• FIDAは中性粒子ビームを利用した高エネルギー粒子計測方法.
• 荷電粒子としてプラズマ中に存在する高エネルギー粒子が中性粒得ることができる.
s (r,
λ)

∞

• 上式を速度分布関数Fについて解くことを再構成(トモグラフィー)

FIDA信号計算の順過程•
比べ新古典輸送理論の径方向拡散係数が線形回帰問題としての再構成
Fast Ion D Alpha (FIDA)W F を F について解くとはWを入力としてSを出力とする• S = 小さくなるのは，通常の新古典輸送理論
線形回帰問題に帰着される.
39
ではプラズマ中心に存在する特殊な軌道による輸送の効果が含まれていないことが原因の一つと考えられる．

1

−1 0

S = WF• 単一計測視点からの再構成であり, 非常に高難度である.

• 上式を速度分布関数Fについて解くことを再構成(トモグラフィー)• GNETコードの計算による速度分布関数を元に, FIDASIMコードはFIDA信号を計算できる.

速度分布関数S = WFFIDA信号計算の順過程の速度情報を0+
ast + HNB法を開発する.

速度分布関数の再構成

計測方法.
子が中性粒
る.
が偏移(ドッ• ある位置 r を測定するとき, 速度分布関数・空間分布 f (E, p, r) と

ਤ 5.5: ఺ P ʹ͓͚Δϐϯν‫ي‬ಓͱྡ઀͢Δ‫ي‬ಓ• 単一計測視点からの再構成であり, 非常に高難度である.

R0：大半径 a：小半径s (r, λ)• 線形モデルを仮定すると, 最小二乗法によりフィッティング直線の係数が決まる.

• 一般的に, 観測は積分操作を伴うため観測量Sは物理量Fよりも低次元となる.

• 低次元情報から高次元情報を復元しようとするため, 解の不定性が生じる.

A BOUT++ extension for full annular tokamak edge turbulence simulation
SETO Haruki1), DUDSON Benjamin Daniel2,3), XU XueQiao2), YAGI Masatoshi1)
1) QST, 2) LLNL, 3) U. York
It is one of key issues for ITER and DEMO to understand the dynamics of edge localized modes (ELMs)[1] in metal wall devices. This is because intermittent
heat fluxes triggered by ELMs should be avoided or
mitigated to levels low enough to remain within heatload constraints on plasma facing components [2].
In H-mode discharges with moderate to high D2
gas puff rate in JET ITER-like wall (JET-ILW),
which is required for W control for stable H-mode
operations, the increase of neutral gas pressure degrades pedestal temperature, and the operation point
is in ballooning mode branch where the high-n pressure-gradient driven instabilities determine the operation point like type-III ELMs [3]. Here n is the toroidal mode number. This is qualitatively different
from type-I ELMs driven by moderate-n peeling-ballooning mode.
On the other hand, it is reported that a highly reproducible stationary small ELM regime is observed
in EAST experiments, which is suitable for future
steady-state fusion devices [4]. The small ELM is akind of grassy ELM, but its operation point is in the
low-n peeling mode branch rather than the high-nballooning branch like grassy ELMs in JT-60U [5].
The nonlinear simulations show that the key mechanism for the grassy ELM in EAST is the expansionof the peeling boundary due to radially localized
steepening of pedestal pressure gradient by a radiallylocalized collapse. Further verification study is required in both experiments and simulations to elucidate the mechanism of grassy ELM in EAST.
It is therefore important for understanding physics
of ELMs in metal wall devices including ITER and
DEMO to develop a nonlinear reduced MHD codesolving off-resonant (n=0), low-n (n~1), moderate-n(n>1) and high-n (n>>1) modes simultaneously. In
this work, an edge MHD/turbulence code for solving
interplay among off-resonant, low-n, moderate-n andhigh-n modes using BOUT++ framework [6] is developed.
BOUT++ employs both field-aligned coordinates[7] and flux-surface coordinates for tokamak edge
simulations in combination with the radial derivative
method and the shifted metric method to simulateballooning mode instabilities with reasonable computational cost and high accuracy. It has provided
qualitative understandings on ELMs triggered bymiddle-n ~ high-n ballooning modes [8]. Taking

fluctuation-driven n=0 and very low-n flow andmagnetic field however has been an open issue dueto the numerical problem in the flute-ordered Poisson solver required for BOUT++ coordinates.
These limitations have been partially removed by
the flute-ordering-free Poisson solver designed for
n=0 mode [9] which enables BOUT++ to simulate aninterplay between fluctuation driven n=0 net flow
and MHD turbulence during ELM crash triggered bymiddle-n ~ high-n ballooning modes. In our previous
works [10,11], it is revealed that the zonal net flow
generated during the pedestal collapse can trigger
subsequent turbulence bursts accompanied with cyclic oscillations between pressure gradient, zonal
flow and turbulence and enhance turbulent energy
transport after the pedestal collapse. This framework
however has a problem in handling very low-nmodes and simulation domains are still limited to annular wedge torus domain.
To simulate tokamak edge MHD/turbulence infull annular torus domain, a novel hybrid Poissonsolver consisting of a flute-ordering-free Poissonsolver designed for n=0 and low-n modes and theconventional flute-ordered Poisson solver for moderate-n and high-n mode is developed. The numerical detail of the hybrid Poisson solver is left for future publications.
The hybrid Poisson solver is tested with a simulation of ELM crash triggered by resistive ballooning
mode (RBM) and drift wave (DW) instabilities in afull annular torus domain with shifted circular crosssection. Here the simulation setup including the governing equations and the equilibrium but the grid resolution and the toroidal domain length is same as isreported in Ref [11]. For the computational domain,
the radial grid number is 512 for 0.6 ≤ 𝑥 ≤ 1.2, theparallel grid number is 128 for 0 ≤ 𝑦 ≤ 2𝜋 and thebinormal grid number is 256 for 0 ≤ 𝑧 ≤ 2𝜋 respectively. For the hybrid Poisson solver, the fluteordering-free Poisson solver is applied from n=0 to
n=4 modes, and the flute-ordered Poisson solver isapplied from n=5 to n=80 modes and n>80 modes
are numerically removed respectively.
Figure 1 shows the dispersion relation of initiallyunstable modes where n=32 is the most unstablemode in this simulation, and Fig. 2 also shows thetime evolution of perpendicular kinetic energy spectra with respect to toroidal mode numberrespectively. In this simulation, the large energy loss
by ELM driven by RBM and DW occurs at t~250tA.
The time evolution of perpendicular kinetic energy spectra shows that the initially unstable modes
get saturated at t~150tA and the very low-n modes(n=1~3) are then driven in prior to the zonal mode(n=0) via nonlinear couplings among initially unstable modes as is shown in Fig.2(b) and Fig. 2(c). Afterthe saturation of very low-n modes and zonal mode,
the large pedestal collapse triggered by the initiallyunstable high-n modes (n~32) and the moderate-n
modes (n~20) in Fig. 2(d), and the zonal mode becomes dominant after the pedestal collapse in Fig.
2(e). This multi-phase collapse is not observed in ourprevious works [10, 11] using 1/5th annular wedge
torus domain due to lack of very low-n modes.
In summary, the hybrid Poisson solver has beendeveloped for full annular tokamak edge simulationin BOUT++ framework. It has been demonstrated
that the usage of full annular torus domain canchange the dynamics of ELM driven by high-n RBM
and DW modes. The detailed analyses such as energy transfer among toroidal modes, the spatial
structure of low-n modes, are left as future works.

Fig. 1: dispersion relation of initially unstable RBM+DW instability, where themost instable mode is n=32

Acknowledgements:
This work was partially supported by JSPS KAKENHI Grant No. 20K14448 and performed under
EU-Japan BA collaboration, under JIFT collaboration between QST and LLNL and under the auspicesof the U.S. Department of Energy by LLNL underContract No. DE-AC52-07NA27344. Computations
were carried out on JFRS-1 supercomputer at International Fusion Energy Research Centre Computational Science Cenre (IFERC-CSC) and on HPE
SGI8600 supercomputer at QST and JAEA.
References:[1] H. Zohm, Plasma Phys. Control. Fusion 28,(1996) 1213
[2] J.P. Gunn et al. Nucl. Fusion 57, (2017) 046025[3] C.F. Maggi et al., Nucl. Fusion 55, (2015) 113031[4] G.S. Xu et al., Phys. Rev. Lett. 112, (2019)
255001[5] N. Oyama et al., Nucl. Fusion 50, (2010) 064014[6] B.D. Dudson et al., Comput. Phys. Commun. 180(2009) 1467
[7] M.A. Beer et al., Phys. Plasmas 2, (1995) 2687
[8] X.Q. Xu et al., Phys. Rev. Lett. 105 (2010)
175005[9] B.D. Dudson and J. Leddy, Plasma Phys.
Control. Fusion 59 (2017) 054010[10] H. Seto et al., Phys. Plasmas 26 (2019) 052507[11] H. Seto et al., Contrib. Plasmas Phys. 60 (2020)
e201900158

Fig. 2: (a) time evolution of perpendicular
kinetic energy spectra and time slices ofnormalized perpendicular kinetic energy
spectra at four time slices, where (b) and (c)are in nonlinear phase before the pedestal
collapse, (d) is during the pedestal collapseand (e) is after the pedestal collapse respectively.

Prediction of turbulence temporal evolution in PANTA by Long-short term memory network
College of Industrial Technology, Nihon University
M.Aizawa Caranza, M.Sasaki, H.Minagawa, Y.Nakazawa, Y.Liu, Y.jazima, K.Hara

１. Introduction

The prediction of the plasma behavior including turbulence is an urgent issue for the dynamic control of the magnetic
confinement plasma. In this study, we perform time series prediction using LSTM for multi-scale turbulent time series
observed by Langmuir probes in the linearly magnetized plasma experimental device PANTA, and describe theprediction characteristics when the training conditions of the networks are varied.
２. Prediction Method・Evaluation Method

LSTM (Long-short term memory) is a kind of recurrent neural network[1], which is suitable for time series forecastingof multi-scale fluctuations because it can learn while retaining past information. The LSTM has the network structure
shown in Fig. 1, with multiple units called LSTM Blocks in the middle layer, cells that store past information in theunits, and three gates (input gate, forget gate, and output gate) that control the flow of information. By adjusting thedegree of flow of this gate, storing it in a cell, and inputting the information from the output and the intermediate layer
unit into the next intermediate layer, it is possible to learn based on past information.
In this study, Prediction of plasma turbulence using the waveform of ion saturation current in turbulent plasmameasured by PANTA. This turbulence data is, Short-period drift-wave turbulence is dominant, and long-period abrupt
transport is driven by its nonlinear processes [2]. The prediction of long-period abrupt transport is an important
technique for the control and operation of fusion plasmas.
We used 80% of the first half of the turbulent time series as training data and 20% of the second half for prediction (seeFig. 2). For the evaluation of the learning efficiency, the cross spectrum evaluated from the predicted and observed
values was used to examine the prediction accuracy.
ht-1

ht

メモリ

✖

学習結果を
更新・記憶

Xt-1

➕
⼊
⼒
ゲ
$
ト

忘
却
ゲ
$
ト

At-1
σ

ht+1

σ

✖

tanh

出
⼒
ゲ
$
ト

tanh

At+1

✖

σ

Xt

Fig.1 Structure of the LSTM to be used

Xt+1

Fig.2 Range of training and prediction

３.Result

In the case of normal prediction using LSTM, the dominant drift wave frequency band of this turbulence data is about
10 kHz, and the dominant frequency components could be predicted, but it was difficult to improve the prediction
accuracy for low-frequency components of about 5 kHz. Since the prediction of the low-frequency component, whichcorresponds to nonlinear abrupt transport, is intrinsically important, the high-frequency component and the averaged

component were removed from the prediction. By narrowing down the target frequency, we succeeded in improving theprediction accuracy of the low-frequency component, as shown in Fig. 3.

Fig.3 Prediction results when high frequency components are excluded.
（Top: Comparison of predicted and observed values, Bottom: Spectrum and cross spectrum of both）

４.SummaryTime series prediction of plasma turbulence as a typical multiscale fluctuation is carried out using LSTM network andthe prediction characteristics are clarified. It is found that the prediction of dominant coherent fluctuations is possible,
and that the performance of the prediction of low-frequency fluctuations can be improved by removing the highfrequency and averaged components.
REFERENCE[1] Sepp Hochreiter, Jürgen Schmidhuber (1997) ”LONG SHORT-TERM MEMORY”, Neural Computation 9 (8):1735–1780.
[2] F. Kin, Phys. Plasmas,26, 042306 (2019)

Analysis of three-dimensional structure of MHD instabilities in tokamak
Y. Kuroda1, N. Kasuya 1,2, M. Sato 3, M. Sasaki 4
IGSES, Kyushu Univ. 1, RIAM, Kyushu Univ.2, NIFS3, CIT, Nihon Univ.4
1. Introduction
In tokamak devices, Magneto-Hydro-Dynamics (MHD) instabilities is induced by large pressure gradientand plasma current, which cause confinement degradation and plasma collapse. Therefore, it is essential toevaluate their excitation and saturation mechanisms. Numerical simulations give three-dimensional (3-D)
structures of the plasma instabilities for comparison with experimental measurements [1].
In a previous study [2], simulations of MHD instability have been performed for the tokamak device
PLATO [3]. From the simulation results, ballooning instability and kink instability both can be excited andthe characteristics of the instabilities with different safety factor profiles have been evaluated.
In this study, under the condition when several characteristic instabilities are excited, we extract theircharacteristic structures by Fourier mode decomposition, and dynamic mode decomposition (DMD), which
has recently been used for dynamic structure extraction.
2. Simulation code
In this study, MHD simulations are carried out for tokamak device PLATO by using MIPS code [4]. Since
the simulation requires initial equilibrium data, this routine uses TASK/EQU [5], which allows equilibrium
calculation under free boundary conditions including the fields induced by the external coil current in theexperimental condition. The equilibrium obtained in the calculation is introduced into the MHD simulation
code MIPS to perform a three-dimensional simulation of the MHD instability. The expected plasma
parameters of PLATO tokamak are the followings; major radius 𝑅
plasma volume 𝑉

0.9 𝑚 , electron temperature 𝑇

0.7 𝑚, minor radius 𝑎

0.2 𝑘𝑒𝑉, density 𝑛

1.0

10

0.2 𝑚 ,𝑚 .

3. Analysis method
In a previous study, the mode spectrum analysis was performed using Fourier expansion. Fourier mode
expansion decomposes a waveform into periodic sine and cosine waves with different amplitudes,
frequencies, and phases. However, it is difficult to evaluate the characteristic mode structure in Fouriertransform when there exists frequency temporal variation or spatial localization. Therefore, mode spectrum
analysis is carried out using DMD. DMD focuses on the time evolution of time series data for performing
mode decomposition, Characteristic modes and their corresponding frequencies can be obtained.
4. MHD simulation result
The MHD simulations are carried out under the condition that both ballooning and kink modes are excited.
The used plasma parameters are as follows; central safety factor 𝑞
viscosity 𝜈

10 , resistibility 𝜂

0.8, central beta value 𝛽

10 , thermal diffusion coefficient 𝜒

1.5,

10 . Figure 1 shows time

evolution of the pressure profile. In Fig.1(b), a ballooning mode structure whose perturbation component
exists in the low field side of torus plasma can be seen with finite poloidal mode numbers. The ballooning

mode consists of multiple Fourier modes with poloidal mode number m = 14 - 25 as shown in Fig. 2. In
Fig.1(c), a kink mode structure exists near the plasma center. In this way, both of ballooning and kink modes
are excited in this simulation.

(a)

(b)

(c)

(d)

Fig.1 Time evolution of the pressure profile in thepoloidal cross-section (t = (a) 0, (b) 350, (c) 460, (d)500)

Fig.2 Radial profiles of the mode spectrum(t = 350)

5. Fourier mode expansion result
Next, dynamics of the modes is analyzed with the time evolution of Fourier modes. Figure 3 shows the radial
profiles of Fourier modes when the (a) ballooning and (b) kink mode are excited. The kink mode component
is excited near the center of the plasma and propagates to the outer region. In Fig.3(b), the kink mode reaches
the position where the ballooning mode exists. Then the kink mode is weakened, and propagates to the outerregion with reduction of the amplitude. When the ballooning mode is not excited, the kink mode directly
propagates to the outer region without amplitude reduction (Fig. 3(c)). The mitigation of the mode
propagation occurs near the rational surface q = 1 in this case, but there is no qualitative change with thevariation of the q profile. Therefore, the magnetic configuration is not the cause of this phenomenon.

Fig.3 Radial profiles of the Fourier components at t = (a) 350 (only ballooning mode is dominant) and (b) 490(competition between the ballooning and kink modes) with q0 = 0.8. The case only with a kink mode bychanging the q profile with q0 = 0.5 is also shown in (c).

Mode coupling between the kink and ballooning modes is the candidate. To confirm this, time evolutionsof the mode components are evaluated. When the kink mode reaches at the edge of the ballooning mode
region ρ ~ 0.5, the amplitude of the ballooning mode begins to be reduced, compared with the case without
the ballooning mode. This is related to flattening of the background plasma profile due to the ballooning
mode (quasi-linear effect). On the other hand, the time evolution of the ballooning mode amplitude has bump,when the kink component reaches that position ρ ~ 0.7. This might be the appearance of the energy transferfrom the kink to the ballooning mode. It is suggested that the interaction between the components affects the
propagation.
6. Dynamic mode decomposition result
Next, DMD is applied to analyze the mode structure [6]. By using DMD, it is possible to extract the whole
image of the characteristic structure, which is expressed as the addition of many modes in Fourier mode.
Figure 4 shows the example of DMD mode in the time evolution of the MHD instabilities. In Fig.4, the highwavenumber component and the kink component were extracted. Since the growth rates of the high wave
number component and the kink mode component are comparable, mixed structures of them were extracted
as the DMD modes. Figure 5 shows reconstruction of the time series of the pressure profile from thedecomposition modes. Figure 5(a)and (b) is the reconstruction from a few number of modes and the original
data, respectively. From these results, it is possible to reproduce the characteristics with a limited number of
DMD modes.

Fig.4 The example of DMD mode (localized mode, ballooning mode, mixture of ballooning and kinkmode, mean profile)Fig.5 (a) Reconstruction of time series of pressure profile from a few number of modes, and original data.

7. Conclusion and future plan
In this study, 3-D simulation of MHD instabilities were carried out for the PLATO device. Fourier mode
expansion analysis and dynamic mode decomposition for the perturbation mode structure was carried out inthe time evolution of competition between ballooning and kink instabilities. Interpretation of mode structures
extracted by dynamic mode decomposition and quantitative evaluation of interaction between modes will be
the future task.
References[1] C. Holland, et al, Phys. Plasmas 16 (2009) 052301.
[2] S. Tomimatsu, et al, Plasma Fusion Res. 15 (2020) 1403052.
[3] A. Fujisawa, AIP Conf. Proc. 1993 (2018) 020011.
[4] Y. Todo, et al, Plasma Fusion Res. 5 (2010) S2062.
[5] http://bpsi.nucleng.kyoto-u.ac.jp/task/[6] M. Sasaki, et al, Plasma Phys. Control. Fusion 61 (2019) 112001.

Prediction of temperature proﬁles in helical plasmas by integrated
code coupled with gyrokinetic transport models
S. Toda, M. Nunami, and H. Sugama
National Institute for Fusion Science, Oroshi-cho 322-6, Toki, Gifu, 509-5292, JapanTo realize fusion energy, prediction of turbulent transport in toroidal plasmas is one of the mostcritical issues to be solved. Study of zonal ﬂow eﬀects on plasma conﬁnement improvement isnecessary for accurate prediction of turbulent transport and plasma proﬁles. A large number ofgyrokinetic turbulence transport simulations have been performed for tokamak plasmas and forhelical plasmas [1, 2, 3]. In tokamak plasmas, transport simulation, which is coupled to gyrokinetic analysis results, is performed. Much larger computer resources are consumed to perform
gyrokinetic simulation in helical plasmas than in tokamak plasmas, because the former requires alarge number of mesh points along ﬁeld lines to resolve helical ripple structures. Since combining
nonlinear gyrokinetic simulation with integrated simulation in helical plasmas is not yet easy, there
is a strong need for predictive models which can quickly reproduce nonlinear simulation results.
Heat diﬀusivity models for electron and ion heat transport, and quasilinear ﬂux models for particle
and heat transport have been proposed in a kinetic electron condition for ion temperature gradient (ITG) turbulence in the Large Helical Device (LHD) [2]. The nonlinear simulation results arequickly evaluated by reduced models, because these models are basically functions of the quantity
related with a linear growth rate and the characteristic quantity for linear response of zonal ﬂows
(the zonal ﬂow decay time [4]). The heat diﬀusivity models have been installed into transport
code for simulating evolutions of the plasma proﬁles in the LHD when additional modeling by normalized ion temperature scale length is used [5, 6]. Since the density gradient which is necessary
to construct an eﬀective diﬀusivity model is not accurately evaluated for the hollow and ﬂattened
density proﬁles which are typical in LHD experimental results, particle diﬀusivity is hard to bemodeled. Furthermore, dynamical transport simulation with a quasilinear ﬂux model installed has
not been performed so far, because the quasilinear ﬂux model depends on a cross-phase between
ﬂuctuating potential and temperature (or density) and this cross-phase term is diﬃcult to be modeled by temperature or density gradients. The quasilinear ﬂux formulation enables us to model
particle and heat transport for ﬂattened or hollow density and temperature proﬁles. The nonlinear simulation results can be predicted by the reduced models for the ion heat transport [3]. The
simulation which includes multi-time scales is a challenging task, because the diﬀerence between
energy conﬁnement time (around 1s) and a typical gyrokinetic time scale (around 1µs) is extremely
large. To obtain the temperature proﬁles, the transport simulation with the gyrokinetic transportmodels installed is performed by the integrated code in this article, where the heat diﬀusivity andquasilinear ﬂux models are evaluated by the linear gyrokinetic simulation by GKV code [7] at eachtime step in the evolution of the integrated simulation. The zonal ﬂow decay time is ﬁxed in thetransport simulation. The linear simulations are performed for the various temperature proﬁles toevaluate the reduced transport models at each time step in the integrated simulation.

The reduced turbulent transport model is evaluated by the linear gyrokinetic simulation using
local ﬂux tube code, GKV. The integrated simulation is performed by the TASK3D code. In the
integrated simulation by TASK3D, the heat transport consists of turbulent and neoclassical parts.
The neoclassical diﬀusion coeﬃcient is derived from the DGN/LHD database. The turbulent
transport is evaluated by the reduced models. The linear growth rates are evaluated at tg =50R/vti in the linear gyrokinetic simulation at an equally spaced twenty radial points in the region
0.0 < ρ < 0.8. Here, ρ = r/a, tg is the time period treated in the linear gyrokinetic simulation, k̃xis the normalized radial wavenumber, r is the radial coordinate and a is the minor radius. The ﬂux
tubes of the GKV code are set at each radial point. The value of 50R/vti is typically ≃ 10−5 s in theplasma parameter region studied here. The integrated simulation by the TASK3D is performedat the eighty radial points in the region 0.0 < ρ < 0.8, where the time step in the integrated
simulation, dt is set to 10−3 s. The quantities L related with the mixing length estimate at the
sixty radial points, where the gyrokinetic simulation is not done, are interpolated from those which
are calculated by the GKV at the equally spaced twenty radial points. The value of the linear
growth rate is determined once one cycle of the integrated simulation is performed. The time of
50R/vti (≃ 10−5 s) for the gyrokinetic simulation is one hundredth of dt(= 10−3 s) in the integrated
simulation using the TASK3D. The results of the temperature proﬁle for one cycle of the integrated
simulation are used as initial conditions for the next cycle. The cycle of the integrated simulationis repeated twenty times until the temperature proﬁles become stationary in the parameter region
studied here.
The transport simulation is performed by the integrated code using the reduced gyrokinetic models. Especially, the quasilinear ﬂux model is used in helical plasmas. The turbulent transport isevaluated by reduced models using the linear gyrokinetic simulation results for the kinetic electron
response at each time step in the integrated simulation. By performing the integrated simulation
with the gyrokinetic models installed, multi-time scale simulation can be performed. The multitime scale simulation is done for the high-Ti plasmas in the standard ﬁeld conﬁguration and the
low-Ti plasmas in the inward shifted ﬁeld conﬁguration of the LHD. The stationary temperature
proﬁles predicted by the transport simulation, using the reduced models which reproduce the nonlinear simulation results, are predicted within a diﬀerence of 30% at most from the experimental
results in the LHD, where the ITG mode is unstable. The simulation using the quasilinear ﬂux
model, including eﬀects of zonal ﬂows and kinetic electrons, is performed for the ﬁrst time to predict the temperature proﬁles in the LHD. This is achieved by combining the two codes, TASK3D
and GKV, developed for helical plasmas. How to construct the models [1, 2] by linear simulation
results is applicable to other plasmas in helical devices where modes other than the ITG mode areunstable, and to plasmas in tokamak devices by capturing the dependence of turbulent transport
on the linear properties in the reduced models. By these reduced models, the integrated simulation
with the gyrokinetic model installed will be performed to predict plasma proﬁles and conﬁnement
properties. This is left for future studies.
This work is partly supported by the NIFS Collaboration Research Programs, NIFS20KNST168

(Plasma Simulator), the JSPS KAKENHI Grant No. 21K03514 and the Collaborative Research
Program of Research Institute for Applied Mechanics, Kyushu University, 2021FP-06.

[1] M. Nunami et al., Phys. Plasmas 20 092307 (2013)[2] S. Toda et al., Phys. Plasmas 26 012510 (2019)[3] S. Toda et al., J. Plasma Phys. 86 815860304 (2020)[4] H. Sugama and T. -H. Watanabe, Phys. Plasmas 13 012501 (2006)[5] S. Toda et al., J. Phys. Conf. Ser. 561 012020 (2014)[6] S. Toda et al., Plasma and Fusion Research 14 3403061 (2019)[7] T. -H. Watanabe and H. Sugama, Nucl. Fusion 46 24 (2006)

トーラスプラズマにおける 3 次元データを用いた揺動計測模擬
Numerical diagnostics of plasma fluctuations
by using 3-D simulation data in torus plasmas
N. Kasuya1,2, R. Yoshihara2, T. Ido1,2, Y. Kuroda2, M. Yagi3, M. Nunami4, Y. Nagashima1,2,
A. Fujisawa1,2
1
RIAM, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
2
IGSES, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
3
QST, Obuchi, Rokkasho, Aomori 039-3212, Japan
4
NIFS, Toki, Gifu 509-5292, Japan
E-mail: kasuya@riam.kyushu-u.ac.jp
Plasma turbulence largely affects particle and heat transport in torus plasmas, so it isimportant to clarify its 3-D feature by detailed measurements in experiments. Synthetic
diagnostics uses simulation data to represent experimental measurements and to makequantitative evaluations [1]. Simulations can show global features of plasma fluctuations, andquantitative comparison between experiments and simulations is possible by simulating
measurement signals numerically using simulation data. Simulation models are different eachother, depending on each simulation target, so comparison between models is also importantto capture various characteristics in magnetized plasmas. We are now developing integrated
numerical diagnostic platform, which can be applied to multiple simulation data that havedifferent data formats [2]. Three kinds of perturbation data are analyzed with a numericalheavy ion beam probe (HIBP) routine. The construction of the developed routines andcalculation results using them are described in this presentation.
Heavy ions with 1+ charge is injected into a core plasma and are ionized to 2+ charge ionsby collision. By measuring the energy of the ionized ions, HIBP observes the potential,
density and magnetic field at a local point in core plasma [3]. PLATO tokamak device [4] is
now constructed to begin the detailed plasma turbulence observation [4]. Multiple HIBPs will
be installed in PLATO for simultaneous measurement at different poloidal and toroidal
positions. For numerical HIBP, 3-D potential, density and temperature profiles are necessary.
Mean on the magnetic flux surface and perturbation quantities with 3-D variations makepossible to carry out a simulation of fluctuation measurement. Constant temperature isassumed for simplicity in this time, so the spatial variation of the collisional cross section is
not considered. Here 3-D profiles obtained from MIPS [5]、R4F [6]、FORTEC-3D [7] areanalyzed. All of them are the calculation results in torus plasmas. In general, particularmodels do not always include all the necessary quantities, so insufficient ones aresupplemented with any assumptions. In addition, spatial coordinates may be different, so datatransfer to the coordinates for the numerical measurement is necessary. Data mapping to thecommon format is often used, but here we introduce the object-oriented method in accordancewith each data format to gain the maximum information. MIPS is the 3-D MHD code, whichprovides 3-D data in torus meshes with same intervals. A MHD small collapse driven bypressure gradient is analyzed by using MIPS. R4F is the 4-field reduced MHD code, which
provides 1-D differential data in the radial direction and 2-D spectral expansion data in thepoloidal and toroidal directions. Global MHD mode and micro turbulence are analyzed byusing R4F. FORTEC-3D is the drift-kinetic code, which provides profile data in the flux
surface coordinates for the analyses. Inhomogeneous potential on the flux surfaces is analyzed
from neoclassical transport in helical plasmas by using FORTEC-3D.
[1] C. Holland, et al., Phys. Plasmas 16 (2009) 052301.
[2] N. Kasuya, et al., 28th Int. Toki Conf. (Toki, Gifu, 2019) P2-01.
[3] T. Ido, et al., Rev. Sci. Instrum. 77 (2006) 10F523.
[4] A. Fujisawa, AIP Conf. Proc. 1993 (2018) 020011.
[5] S. Tomimatsu, et al., Plasma Fusion Res. 15 (2020) 1403052.
[6] M. Yagi, et al., Plasma Fusion Res. 2 (2007) 025.
[7] S. Satake, et al., Plasma Fusion Res. 3 (2008) S1062.

Development of a semi-empirical particle and heat transportmodel and improvement in its turbulent saturation rule
E. Narita1 , M. Honda2 , M. Nakata3,4 , M. Yoshida1 and N. Hayashi11National Institutes for Quantum Science and Technology, Naka, Ibaraki 311-0193,
Japan
2
Graduate School of Engineering, Kyoto University, Nishikyo, Kyoto 615-8530, Japan
3
National Institute for Fusion Science, Toki, Gifu 509-5292, Japan
4
The Graduate University for Advanced Studies, Toki, Gifu 509-5292, Japan

1

Introduction

The gyrokinetic-based turbulent transport models, such as TGLF [1,2], are essential to predict
density and temperature proﬁles, but introducing detailed descriptions of turbulence physics
tends to increase the computational cost. To accelerate the proﬁle predictions, a neural-network(NN) based approach has been undertaken (see e.g., [3–6]). Our study is also developing a
NN-based turbulent transport model DeKANIS [7,8]. A turbulent saturation rule employed in
DeKANIS was based on experimental particle ﬂuxes estimated for JT-60U H-mode plasmas, andit was apt to overestimate temperatures. To reduce the overestimation, a diﬀerent saturation
rule is built including experimental heat ﬂuxes.

2

Transport model DeKANIS and its original saturation rule

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ﬀusive) and oﬀ-diagonal (non-diﬀusive) 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 as the ﬁrst oﬀ-diagonal terms of Γ̄e and Q̄e satisfy the Onsager
symmetry, the following relationship is obtained: χ̄e = D̄(CT + 1.5)/CN . 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,χ̄i,eﬀχ̄e,eﬀ LRT nnei TTei . Here, χ̄eﬀ is the
DeKANIS predicts the ion heat ﬂux, which is given as Q̄i = χ̄e,eﬀ
inormalized eﬀective heat diﬀusivity, which is deﬁned when the heat ﬂux is expressed with only
the diagonal term. Once χ̄e,eﬀ is obtained with Q̄e , Q̄i can be constructed through χ̄i,eﬀ /χ̄e,eﬀ .
DeKANIS estimates CT,P,N,HP , D̄, χ̄e and χ̄e,eﬀ /χ̄i,eﬀ with a NN model. The training dataset
has been developed based on 23 JT-60U H-mode plasmas. Turbulence of these plasmas isdominantly driven by the ion temperature gradient (ITG) mode or the ITG/trapped electron
mode. In the dataset, CT,P,N,HP and χ̄e,eﬀ /χ̄i,eﬀ are given by linear gyrokinetic calculations with
GKW [9]. On the other hand, D̄ has been determined by the two diﬀerent saturation rules. The
ﬁrst rule, the semi-empirical saturation rule, calculates D̄ to match it to the the experimental

Table 1. The exponents α and β and proportionality constant c of the scaling formulas.
Base
Γ̄exp
e,turbQ̄expe,turbα1.10081.1990β−2.8024−2.4732c [10−4 ]1.69967.8584turbulent particle ﬂux Γ̄exp
e,turb . The second rule will be shown in later. The output parametersof the NN model are CT,P,N,HP , D̄ and χ̄e,eﬀ /χ̄i,eﬀ , and with the coeﬃcients, χ̄e is computed not
to break the Onsager symmetry.
When DeKANIS is applied to predicting particle and heat ﬂuxes of 14
JT-60U H-mode plasmas that are notused for the training dataset, Q̄e andQ̄i tend to be underestimated in ρ >0.3 by comparison with Γ̄e as shownwith the mean oﬀset (MOFF) in ﬁgFigure 1. The mean oﬀsets of (a) Γ̄e and Q̄e,i and (b)ure 1(a), where ρ is the normalized
ne and Te,i predicted with the semi-empirical saturaminor radius. Here, MOFF of antion rule for 14 JT-60U test cases at each radial point.
arbitrary ξ is deﬁned as MOFF =∑(1/N ) (ξDeKANIS /ξexp − 1). Because of the underestimate, Te and Ti are overestimated when
the density and temperature proﬁles are predicted with the integrated code TOPICS [11], asshown in ﬁgure 1(b). The overestimated temperatures accompanied by steeper temperature
gradients cause an increase in the particle ﬂux through the oﬀ-diagonal terms. Since the valueof the diagonal term of the particle ﬂux needs to decrease to oﬀset the increase in that of theoﬀ-diagonal ones, the density proﬁle tends to be ﬂat, yielding an underestimate of the density.

3

Mixing-length-like saturation rules

The second saturation rule is based on the mixing-length theory, assuming that D̄ is expressedby the quasilinear diﬀusivity γ̄/k̄θ2 characterizing instabilities and the residual zonal ﬂow level
(
)0.5
(
)α
LZF [12]: D̄ = c γ̄/k̄θ2 LβZF . Here, LZF is deﬁned as LZF ≡ KRH γ̄/k̄θ2, where KRH = 1/(1+2
1/2
1.6q /ϵ ) with the safety factor q and the inverse aspect ratio ϵ [13]. When the mixing-lengthlike saturation rule is used, yet another NN model is applied, which forecasts the maximum
linear growth rate γ̄max and the wavenumber kθ with γ̄max based on the GKW calculations. D̄ iscomputed by using γ̄max and kθ given by the NN model with the scaling formula. The exponents
α and β and proportionality constant c of the original scaling formula were derived from a manual
regression analysis for D̄ evaluated based on Γ̄exp
e,turb in the same manner as the semi-empiricalsaturation rule. In this paper, the regression analysis is performed for D̄ estimated to match itto the experimental electron turbulent heat ﬂux Q̄exp
e,turb . Also, α, β and c are optimized with agenetic algorithm. To check the eﬀect of the change in the basic experimental value, the scaling
formula optimized for Γ̄exp
e,turb is also created.
exp
The obtained two scaling formulas based on Γ̄exp
e,turb and Q̄e,turb have α, β and c shown intable 1. Since the two scaling formulas share comparable α and β, they show similar eﬀects of
instabilities and zonal ﬂows. On the other hand, c increases if the basic experimental value isFigure 2. Radial proﬁles of (a) ne , (b) Te and (c) Ti predicted with the semi-empirical saturation
rule and the two mixing-length-like saturation rules for a JT-60U plasma outside the training
dataset. The mean oﬀsets of (d) ne , (e) Te and (f) Ti predicted with the three saturation rules
for 14 JT-60U test cases at each radial point.
exp
Q̄exp
e,turb . When these formulas are used in TOPICS, due to the larger c, the Q̄e,turb -based mixinglength-like saturation rule predicts lower temperatures as shown in ﬁgure 2(b) and (c) with thelarger heat ﬂuxes in comparison with the semi-empirical saturation rule and the Γ̄exp
e,turb -based
mixing-length-like one. As for the density proﬁle, the three saturation rules gives similar results(ﬁgure 2(c)). The tendencies are also observed in the MOFF proﬁles for the 14 JT-60U H-mode
plasmas outside of the training and regression datasets as shown in ﬁgure 2(d)–(f). One can ﬁndthat the temperature overestimation is reduced by introducing Q̄exp
e,turb in the saturation rule.
Further improvement in the saturation rule is, however, required, because the current scaling
formulas show weak dependence of D̄ on γ̄/k̄θ2 as α ∼ −0.5β. It contradicts the well-known
mixing-length rule. Furthermore, we need to investigate the nature of the inconsistency of c
between the two scaling formulas, since the inconsistency indicates that Q̄e /Γ̄e estimated by thequasilinear theory does not match the experimental value.

References[1] G.M. Staebler et al., Phys. Plasmas 12 102508 (2005).
[2] G.M. Staebler et al., Nucl. Fusion 61 116007 (2021).
[3] M. Honda and E. Narita, Phys. Plasmas 26 102307 (2019).
[4] M. Honda and E. Narita, Plasma Fusion Res. 16 2403002 (2021).
[5] O. Meneghini et al., Nucl. Fusion 61 026006 (2021).
[6] K.L. van de Plassche et al., Phys. Plasmas 27 022310 (2020).
[7] E. Narita et al., Nucl. Fusion 59 106018 (2019).
[8] E. Narita et al., Nucl. Fusion 61 116041 (2021).
[9] A.G. Peeters et al., Comput. Phys. Commun. 180 2650 (2009).
[10] E. Narita et al., Plasma Phys. Controlled Fusion 60 025027 (2018).
[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] M.N. Rosenbluth and F.K. Hinton, Phys. Rev. Lett. 80 724 (1998).

Spatio-temporal analysis of plasma turbulence in PANTA byusing singular value decomposition
R. Takahashi 1, M. Sasaki1,2,3
College of Industrial Technology, Nihon University

１

Research Institute for Applied Mechanics, Kyushu University, 816-8580, Japan

２

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

３

1 Introduction
In turbulent transport in magnetic confined plasmas, turbulence is excited bytemperature and density gradients, and the excited turbulence causes transport andaffects the temperature and density gradients. In addition, turbulence forms structures
such as flow and rotation, and the formation of these structures suppresses the turbulence,
which is a dynamic competitive state that needs to be understood. The conventional
analysis method, Fourier transform, uses trigonometric functions as basis functions,which is effective for linear data, but difficult to characterize time-varying data due to its
large degree of freedom. In this study, singular value decomposition was used to analyze
magnetic confinement plasma data, which is suitable for extracting features with small
degrees of freedom because the optimal basis can be obtained.

2 ModelSingular value decomposition of 64 channels of spatio-temporal plasma turbulence data
during 0.2[s] acquired by the linearly magnetized plasma system PANTA has been
analyzed using a computer. The resulting data are further analyzed.

3 Results
In Fig. 1, the vertical axis shows the spatial structure and the horizontal axis shows thetime. Since the spatial data is invariant in 32ch, the vertical axis shows the combination ofModes decomposed by SVD. The last column of each magnetic field is the singular value
Σ of the characteristic mode multiplied by the right singular vector Vt. In the case of
600G, the zonal flow seen in the Mode with small contribution seems to affect the ModeFig.1Results of SVD

with large contribution. On the other hand, at a magnetic field of 1200G, zonal flow-likephenomena are observed, but there is little effect on the Mode with high contribution. The
reason for this is that in the 600G case, there are low-periodic areas of high intensity, whilein the 1200G case, there are areas of low intensity. It was thought that because of the presenceof high periodic intensity areas, the next change would occur before it could affect the highcontribution Mode, and thus the effect would not be felt.
Fig. 2 shows the standard deviation (top)and frequency (bottom) of the zonal flowlike structures, which increase with higher
magnetic fields, and may be related to the
hypothesis.

Fig.2 Frequency of zonal flow _ like

Effect of impurities on improved confinement in tokamak discharges
Sho Nonaka1, Naohiro Kasuya1,2, Atsushi Fukuyama3,
Shota Mochinaga1, Masatoshi Yagi4
IGSES, Kyushu Univ.1, RIAM, Kyushu Univ.2, Kyoto Univ.3, QST.4
1 Introduction
Presently, International Thermonuclear Experimental Reactor (ITER) is under construction. For
realizing the fusion reactor, achievement of high performance plasma discharge is necessary. High
confinement discharges have been obtained, which are called high βp mode [1], reversed shear mode[2], and so on, with formation of internal transport barrier (ITB). On the other hand, high Z impurities
degrade plasma performance, as tungsten planned to be used for the divertor plate of ITER. In thisresearch, the effect of tungsten impurity is investigated in the reversed shear mode by using integrated
transport simulation code TASK [3]. Dependence on impurity influx is evaluated for the heating power
threshold and internal transport barrier performance.

2 Analysis
2.1 L-mode discharge
For the simulation of L-mode discharge, additional RF heating is applied at t = 1.0s, which is assumed
to be equally absorbed by electrons and ions. Tungsten impurity comes from the surface of the plasma
with influx in,W to give a peaked profile of line radiation, which affects to form a hollow temperature
profile in t = 1.5 - 3.0s as in Fig.1. On the other hand, the effect of the impurity is small in the steadystate, due to reduction of impurity accumulation by RF heating. The density profile of the impurity
depends on used turbulent diffusion models. Increase of the temperature which caused by RF heating
leads to increase the turbulent diffusion coefficient (Fig.2).

Fig.1 Temperature profile in the case of ΓW,in＝0.0m-2s-1 and ΓW,in＝2.0×1018 m-2s-1

Fig.2 Radial profile of the turbulent diffusion coefficient.

2.2 Reversed shear mode discharge
For the simulation of reversed shear mode, additional RF heating is applied as the same way as the
L-mode discharge, and the plasma current is increased from t = 1.0s to 2.0s in addition. In thissimulation, the obtained density profile of the impurity is rather flat as in Fig.3, due to large turbulent
diffusion coefficient as in Fig.4.
Dependence of the temperature on RF heating power is investigated. When the impurity influx is
ΓW,in＝1.0×1018 m-2s-1, the electron temperature is largely changed between PRF = 2MW and 3MW(Fig.5). With PRF = 2MW, the temperature at the center decreased after t = 2.0s and the profile becomes
hollow, while the plasma is sustained with PRF 3MW. The difference between two discharges isunderstood by comparison between the RF heating and line radiation powers. In the case when the RF
heating is greater than the line radiation, large temperature plasma can be sustained.

Fig.3 Density profile of the impurity

Fig.4 Radial profile of turbulent diffusion coefficientFig.5 Dependence on the heating power of electron temperature at the plasma center.

2.3 Comparison of L-mode and reversed shear
The L-mode described in Sec. 2.1 is compared with the reversed shear mode on the condition when
the RF heating is 10MW, the impurity flux is 2.0×1018 m-2s-1, the plasma current is 3MA. The effectof the impurity is small both in the L-mode and the reversed shear mode with rather large input power.
In the case of the reversed shear mode, the turbulent diffusion coefficient is greater than that in the Lmode due to the larger temperature and its gradient inside of the ITB. The line radiation of the reversed
shear mode is smaller than that of the L-mode, because the impurity density is reduced by the larger
turbulent diffusion. The turbulent diffusion coefficients in this simulation comes from the used
turbulent model; mBgB model. Another model gives the other results [4].

3 Conclusion
The effect of tungsten impurity is investigated in the reversed shear mode by using integrated
transport simulation code TASK. In the L-mode discharge, peaking of line radiation by impurities
affects the temperature profile in the transient phase. On the other hand, the effect of impurities issmall in the steady state, because impurity density and its radiation are reduced by RF heating. In thereversed shear mode, dependence of the temperature on RF heating power is investigated with
impurity influx. The case when the RF heating power is smaller than the line radiation by impurities,
a high temperature plasma is not sustained. The line radiation in the reversed shear is smaller than thatin the L-mode, because impurity density is reduced by lager turbulent diffusion due to ITB formation.
The simulation result depends on turbulence model.

References[1] Y. Koide, et al., Phys. Rev. Lett. 72, 3662 (1994)[2] T. Fujita, et al., Phys. Rev. Lett. 78, 2377 (1997)[3] M. Honda, et al., Nucl. Fusion 46, 580 (2006)[4] Erba, et al., Nuclear Fusion 38, 1013 (1998)

円筒プラズマにおけるイオン質量数とドリフト波乱流構造
Drift wave turbulent structure depending on ion mass number
in cylindrical plasmas
N. Kasuya1,2, Y. Imahashi2, H. Todoroki2, M. Ishida2 and M. Yagi3
1
RIAM, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
2
IGSES, Kyushu University, Kasuga, Fukuoka 816-8580, Japan
3
QST, Obuchi, Rokkasho, Aomori 039-3212, Japan
E-mail: kasuya@riam.kyushu-u.ac.jp
Variation of ion species in discharges affects performance of plasma confinement. Basicexperimental devices with cylindrical configuration are used to investigate fundamental
mechanism of plasma turbulence, such as formation mechanism of turbulence structures [1]and nonlinear interaction between plasma instabilities [2]. Effects of ion species have been
also studied in the basic plasmas [3]. In simulations, variations of characteristic mode
numbers and formed structures in plasma turbulence are predicted between different dischargegases [4]. Nonlinear simulations for PANTA linear device [5] using turbulence code NLD [6]have shown that a large number of modes with larger axial mode numbers become unstable inthe case of smaller mass ions, and their nonlinear couplings gives difference in formed
turbulent structures. Figure 1 shows comparison of snapshot of turbulent states in argon andhelium plasmas. There are preferable azimuthal and axial mode numbers corresponding to theion mass number.
For the calculation of resistive drift wave instability in cylindrical plasmas withhomogeneous magnetic fields in the axial direction, the following set of three-field reduced
fluid equations is used [8]. The time and distance are normalized by the ion cyclotron
frequency Ωci and the effective Larmor radiusevaluated by using the electron
temperature Te, respectively. Neutral particles exist even in the center of this rather low
temperature plasma, so the effect of neutrals is included by the ion-neutral collision in terms.
In addition, only charge number Z = +1 ions are considered with this range of the temperature.
Dependence of the linear growthrate is obtained numerically as in Fig. 2. The simulation
parameters are the followings; magnetic field B = 0.1 T, electron temperature Te = 3 eV, plasma
radius a = 7 cm, device axial length  = 4 m, electron-ion collision frequency νei/Ωci = 300,
electron-neutral collision frequency νen/Ωci = 10, viscosities W = V = 10-4, N = 10-2 and Ln =a / 5. There are several (m, n) modes, which have positive growthrate (unstable modes) with acertain in and A, and among them the mode with maximum growthrate is picked up, which
has mode number (mp, np). Here, m and n are azimuthal and axial mode numbers, respectively.
Larger in makes modes stable, and heavier ion makes mp and np smaller. (mp, np) = (15, 11)mode is the most unstable with helium plasma, though (3, 1) mode is with argon plasma,
when in = 0.04.
Linearization
eigen-equation

andsimplificationwith0  0 r  V0  V0 r  0

give

an

T  1i   P  in  *   0 ,
(1)

T 
T 


where P  M i k z 2  e is related to axial wave number kz, *  k Ln 1 is the diamagnetic

 2  i in  iP

frequency, M i  Amu me is the ion mass ratio, A is the ion mass number, mu is the atomic
mass constant, Ln 1   dN 0 dr is inverse of the density gradient length, and T  k 2  ikr Ln 1 .
Here viscosities are put to be zero W = V = N = 0 for simplicity. The solution of Eq. (1) is
2

 
T 1
 T 1
 T 1  in  2i *   P 2 i in  iP
  in 2  2 P 

(2)
T
T 
 T  .
 T2In these parameters, in and P have the same order of magnitude, so the expression is little bitcomplex. To evaluate the most unstable kz, an analytical expression with in << P << 1 isobtained as

T 1
(3)
 ~ iP
 iP * .
2T
TSolution (3) has the maximum growthratek 2
1 k Ln 1
1
,
when
(4)
Im  Pk L 1 .
2  n
2
2
4 1  k2 1  k Then, k2 = 1 + kr2 gives the maximum of Eq. (8) with
1 1  2k r 2
Ln 1
1
,
when
P

Ln 1 .
Im  
32
32
2
2
8 1  k 
8 1  k 

(5)

r

r

For the radial direction, smaller kr makes the mode unstable (kr  0), and the value isdetermined by the boundary condition, when the global mode structure is considered. Fromthese expression, dependence of a spectrum peak on the ion mass Mi is given to bem p  A 1 2 from k ~ 1(6)
andn p  A 3 4 from

M i kz 2

1

1~ Ln ,(7)e8
considering the normalization with s. In this way, preferable k and kz are obtained, whichmakes the plasma most unstable. Perpendicular and parallel flow balance and delay of the
response by collision is the cause of the resistive drift wave instability. In the balance,competition of the convective derivative term and the perpendicular drift (especiallypolarization drift) term is the key to determine the wave number spectrum.

Fig.1 Numerical simulations of potential perturbations in the saturated states of plasma
turbulence in cylindrical plasmas. Comparison between argon and helium discharges
are shown.

Fig.2 (a) Dependence of the maximum linear growthrate on ion mass number A and collision
frequency in. The dashed line shows the unstable boundary where the growthrate iszero. The maxima are given with different (b) azimuthal mp and (c) axial np, dependingon the parameters.
Acknowledgements
Authors acknowledge discussion with Dr. M. Sasaki. 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. Yamada et al., Nature Phys. 4, 721 (2008).
[2] H. Arakawa et al., Plasma Phys. Control. Fusion 53, 115009 (2011).
[3] T. Kobayashi et al., Phys. Plasmas 27, 062309 (2020).
[4] S. Tomimatsu et al., Plasma Fusion Res. 15, 1403052 (2020).
[5] S. Inagaki et al., Sci. Reps. 6, 22189 (2016).
[6] N. Kasuya et al., J. Phys. Soc. Jpn. 76, 044501 (2007).