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

第 11 回核燃焼プラズマ統合コード研究会
11th Burning Plasma Simulation Initiative (BPSI) Meeting(Ver.2)
日時：2013 年 12 月 19 日(木)－20 日(金)場所：九州大学筑紫キャンパス 応用力学研究所 2 階大会議室(20 min talk+5 min discuss or 15min talk+5min discuss)
12 月 19 日(木)9:00 – 9:10 はじめに
(座長：福山)9:10 – 9:35 講演 1-1

福山 (京大)若狭 (RIST)

Heat Transport Simulation in LHD Plasmas Using the TASK3D code
9:35 – 10:00 講演 1-2
本多 (原子力機構)
JT-60U における新古典トロイダル粘性の影響
10:00 – 10:25 講演 1-3
宮戸 (原子力機構)
トカマクプラズマにおけるＩＴＧ乱流とコールドパルスのシミュレーション
10:25 – 10:40 休憩(座長：宮戸)10:40 – 11:05 講演 1-4

矢木 (原子力機構)

Simulation study of nonlocal transport from edge to core in tokamak
plasmas11:05 – 11:30 講演 1-5
村上 (京大)LHD における高 Ti プラズマの統合シミュレーション
11:30 – 11:55 講演 1-6
登田 (核融合研)
ヘリカルプラズマでのイオン温度勾配乱流における熱拡散係数のモデリング11:55 – 12:15 講演 1-7
西岡 (京大)ヘリオトロン J の NBI プラズマにおける新古典イオンフローのモーメント法を用いた数値解析
12:15 – 13:15 昼休み(座長：糟谷)13:15 – 13:40 講演 2-1

林 (原子力機構)
1

ペレットによる ELM 制御の統合シミュレーション
13:40 – 14:05 講演 2-2
松山 (原子力機構)
ペレット入射プラズモイドのドリフト機構のモデリングと実験との比較
14:05 – 14:30 講演 2-3
河野 (九工大)
イオンサイクロトロン加熱におけるプラズマ波と RF シースの相互作用
14:30 – 14:50 講演 2-4
羽田 (京大)One-Dimensional Transport Analysis of ECRH Assisted Plasma Start-up
in JT-60SA
14:50 – 15:05 休憩(座長：林)
15:05 – 15:30 講演 2-5

滝塚 (阪大)Enhanced neo-classical resistivity due to the m/n = 1/1 MHD-mode
deformation for tokamak hybrid scenario
15:30 – 15:55 講演 2-6
佐藤 (核融合研)
陰的時間積分法を用いた MHD コードの開発
15:55 – 16:15 講演 2-7
糟谷 (九大)乱流計測シミュレータを用いた非拡散輸送過程の数値診断
16:15 – 16:35 講演 2-8
佐々木 (九大)直線型装置における乱流構造形成シミュレーション
16:35 – 16:50 休憩(座長：滝塚)16:50 – 17:10 講演 2-9

小菅 (九大)炉心－周辺結合領域における輸送モデリング17:10 – 17:30 講演 2-10
大澤 (九大)タングステン空孔中の水素の特異構造とその材料への応用
17:30 – 17:50 講演 2-11
東 (京大)
PARASOL コードによるトカマク周辺プラズマの粒子シミュレーション
17:50

散会

19:30 – 22:00 懇親会 (十徳や筑紫口店にて)

2

12 月 20 日(金)9:00 – 9:05 事務連絡(座長：矢木)
9:05 – 9:30 講演 3-1
福山 (京大)ITER 統合モデリングの進展
9:30 – 9:55 講演 3-2
横山 (核融合研)
統合輸送解析スイート TASK3D（a,p）開発と LHD 実験適用の進展
9:55 – 10:15 講演 3-3
金森 (京大)TASK/TX によるトカマク粒子輸送シミュレーション
10:15 – 10:20 休憩(座長：横山)10:30 – 10:55 講演 3-4
藤堂 (核融合研)
DIII-D におけるアルフベン固有モードと高速イオン輸送のシミュレーション
10:55 – 11:20 講演 3-5
Lesur (九大)Scheme-independent error in entropy conservation in turbulent kinetic
simulations
11:20 – 11:45 講演 3-6
小関 (原子力機構)
ITER 遠隔実験センター
11:45 – 12:00 まとめ矢木
12:00 – 13:00 昼休み
13:00 – 17:00 核融合エネルギーフォーラムサブクラスター会合
17:00

散会

3

Heat Transport Simulation in LHD Plasmas
Using TASK3D code
A. Wakasa1, A. Fukuyama2, S. Murakami2, N. Takeda2, and, M. Yokoyama3
1

Research Organization for Information Science and Technology, Tokai, 319-1106, Japan2

Department of Nuclear Engineering, Kyoto University, Sakyo-ku, Kyoto 606-8501, Japan
3

National Institute for Fusion Science,322-6 Oroshi-cho, Toki 509-5292, Japan

1. Introduction
TASK3D[1] is an integrated transport simulation code for the helical plasmas. It is under development
in collaboration with Kyoto University and NIFS and has been used to analyze a variety of LHD
experiments so far. Using TASK3D, we have performed self-consistent calculation of the heattransport and heating profile for parameters in experiments (Experimental analysis) and predictive
simulations assuming a variety of NBI heating conditions (Predictive analysis). In this study, in orderto improve the accuracy of the turbulent transport model in TASK3D, we have made comparison andthe validation of gyro-Bohm based transport models with LHD experimental results (14th, Sep. 2011).
2. Integrated simulation code, TASK3DTASK3D has a modular structure and each module describes different physics phenomena. Figure 1
shows the flow diagrams of the 1D heat transport simulation using TASK3D. Here, TR is the 1D
diffusive transport module, and solves particle transport equation, heat transport equation, andmagnetic field equation. The neoclassical transport coefficients are calculated by the neoclassicaltransport database, DGN/LHD[2], and the radial electric field is calculated by the ER module,
according to the ambipolar condition. The NBI power deposition profiles are calculated by FIT3D
module in response to changes
inthespatialtemperature

profilesandofdensitycalculated by the TR module.
Thus,

we

can

perform

self-consistent simulations ofheat transport in NBI plasmas.
3. Predictive simulation withsimple gyro-Bohm transport
model
Fig. 1. The flow diagrams of 1D heat transport simulation using TASK3D

We assume that the thermal diffusion coefficients are given by the sum of a turbulent term TB and aneoclassical term NC, = TB + NC. We also assume eTB=iTB and adopted a simple gyro-Bohm
TB
model,  gyroBohm
 CgyroBohm T eB  i a , compatible with International Stellarator Scaling 95，

ISS95[3]. Here, CgyroBohm is a constant factor, which is adjusted so that the simulation results well
reproduce experimental observations. In this study, CgyroBohm is determined as CgyroBohm = 21.2 tominimize the relative mean square error between experimental and simulation results of the reference
shots. Using this value of CgyroBohm, we have performed predictive heat transport simulation in LHD as
shown Fig. 3(a). In order to validate the simulation results, the LHD experiments in the similar NBI
heating condition and density profiles to the TASK3D simulation were performed as shown Fig. 3(b).
Te obtained by the predictive simulation is in good(a)4

obtained by the predictive simulation is about 25% lower
than the experimental results in the core region. In recent
high Ti LHD experiments, the Ti profiles obtained by
TASK3D with this simple gyro-Bohm model are alsorather low in the plasma core region.
5. Summary
We have been developing the integrated simulation code
TASK3D for non-axisymmetric plasmas. We have

Temperature [keV]

agreement with the experimental results, while Ti

Te

3
2

Ti

1
0
0

0.2

0.4

0.6

0.8

1

r/a

chosen the constant factor for a simple gyro-Bohm(b)transport model, Cgyro-Bohm, using the recent LHD
experimental results. Then, we have performed a
predictive simulation and found that the predicated Te
is in good agreement with the experimental results,
while the predicted Ti is less peaked compared withthe experimentally observed profile.
Acknowledgments

This work is supported in partby Grant-in-Aid for Scientific Research (S) (20226017)
from JSPS and by the NIFS Collaboration Research

Fig. 3. (a) Rsults of Predictive simulation by
TASK3D and (b) Experimental results underclose condition with TASK3D simulation.

program.
References
[1]
[2]
[3]

A. Wakasa et al., 23rd IAEA-FEC (2010).
A. Wakasa et. al., Contrib. Plasma Phys. 50, No. 6-7, 582 585 (2010)
U. Stroth et al., Nucl. Fusion 36, 1063 (1996).

Effects of the neoclassical toroidal viscosity in JT-60U
M. Honda1, S. Satake2, Y. Suzuki2, G. Matsunaga1, K. Shinohara1, M. Yoshida1, H. Urano1,
S. Ide1 and N. Hayashi1

1Naka Fusion Institute, Japan Atomic Energy Agency
2
National Institute for Fusion Science

The integrated code, the TOPICS suite, has recently gained the capability to

self-consistently calculate the radial electric field Er and the evolution of poloidal and toroidal

flows [1]. With the aid of the orbit-following Monte Carlo code OFMC, toroidal rotation
driven by neutral beam injection (NBI) can be calculated in plasmas with the toroidal field(TF) ripple and thus the ripple effects on fast ions regarding a torque has already been

investigated [1]. Recently the neoclassical toroidal viscosity (NTV) in tokamaks acting on
thermal ions has attracted attention in that the NTV damps toroidal rotation towards a certain

level of the rotation velocity, called the offset rotation. The NTV is a kind of the

non-axisymmetric 3D effect in tokamaks due to the toroidally-asymmetric perturbed magnetic

field that stems from not only the discrete TF coils but also their misalignment, MHD

activities and the resonant magnetic perturbation. Given that in reality tokamaks are notaxisymmetric, the NTV is ubiquitous and is therefore an indispensable piece of toroidal
rotation predictions as a torque source especially localized in the edge region. In contrast,

transport simulations with the NTV have seldom been performed. This is due partly to its

complicated 3D nature and a lack of tools to numerically estimate it. A drift-kinetic δf Monte

Carlo code, FORTEC-3D, which has originally been developed for neoclassical transport
analyses in heliotrons and stellarators, is able to compute the NTV [2,3], including effects of

Er in tokamaks [4]. A numerical tool like FORTEC-3D has many advantages in that it does
not essentially require any assumptions that are usually required in the course of thederivation of the analytic NTV formula and can adopt a real equilibrium geometry. The

magnetic equilibrium including the plasma response by the perturbed magnetic field in
vacuum can be computed by the 3D equilibrium code, VMEC [5].

The cooperative framework of TOPICS, VMEC and FORTEC-3D has beendeveloped for the NTV analysis [6], as in detail shown in Fig.1. In this work, it is applied toJT-60U experimental analyses. Since JT-60U had somewhat larger TF ripple amplitude due
partly to the circular TF coils, it can be speculated that the NTV potentially plays a
non-negligible role in the formation of the toroidal rotation profile. We therefore compute theNTV profile in JT-60U L-mode discharges w/ and w/o the ferritic steel tiles (FSTs) and study
its impact on toroidal rotation by TOPICS simulations. Note that the boundary is set at ρ=0.95.
As shown in Fig. 2(b), the NTV torque indeed lies in the edge region and its amplitude iscomparable to or greater than the j×B torque due to the NBI. It is found in Fig. 2(a) that

including the NTV largely improves the reproducibility of predicted rotation profiles, a factwhich indicates that the NTV is imperative as one of the major torque sources.

Figure 1. The cooperative framework of TOPICS, VMEC and FORTEC-3D.

Figure 2. (a) Carbon toroidal rotation profiles of the experimental observation (points with error bars anddotted lines), the predictions w/ and w/o the NTV (solid and dashed lines, respectively).
(b) Torque profiles of the collisional and the j×B transfers from the NBI and the NTV[1] M. Honda et al., Nucl. Fusion 53 (2013) 073050.

[2] S. Satake et al., Plasma Phys. Control. Fusion 53 (2011) 054018.
[3] S. Satake et al., Phys. Rev. Lett. 107 (2011) 055001.
[4] S. Satake et al., Nucl. Fusion 53 (2013) 113033.

[5] S. P. Hirshman and O. Betancourt, J. Comput. Phys. 96 (1991) 99.
[6] M. Honda et al., submitted to Nucl. Fusion.

Simulation study of ion temperature gradient driven turbulence with cold
pulse in a tokamak plasma
N. Miyato and M. Yagi
Japan Atomic Energy Agency, Rokkasho, Aomori
Nonlocal plasma responses which cannot be explained by the conventional local transportmodels are widely observed in toroidal plasmas [1]. For example, electron temperature in thecore region rapidly increases by edge cooling before its effect reaches the core region
diffusively. Recently, a nonlocal response of plasma to edge density sources is found in global
fluid simulations which are based on the 4-field reduced magnetohydrodynamics (RMHD)
model [2, 3]. In the nonlocal plasma response/transport found in the RMHD simulations both
0 Fourier modes are responsible,

nonlinear and toroidal couplings between axisymmetric,in particular

1, 0 pressure fluctuations play a key role, where m and n are poloidal

and toroidal mode numbers, respectively. Although the resistive ballooning modes are unstable
at the edge, they do not play an important role in the nonlocal transport.
In the above RMHD simulations no turbulence exists in the core region where the nonlocal
response appears. However, ion temperature gradient (ITG) driven turbulence is important fortransport of tokamak core plasmas. The

1, 0 pressure perturbations are the ingredient ofoscillatory zonal flows (ZFs) or geodesic acoustic modes (GAMs) also [4]. The GAM
frequency observed in global ITG simulations deviates from its local frequency and has some
identical value over a certain radial distance [5]. Therefore, the nonlocal transport observed in
the RMHD simulations may be affected by the ITG turbulence. We have modified the global
Landau-fluid code used for the ITG turbulence simulations to include the toroidally elongatedsource used in the RMHD simulation [2]
2

exp
where

is the source amplitude,,

cos cos2∆

2

sin sin

1

0.8, 0 and Δ=0.1 represent the source locationin the poloidal plane and the half width of source, respectively. As a first step we perform
electrostatic ITG simulations based on the 3-field Landau-fluid model in which the source isimplemented as a sink or cold pulse in an ion temperature equation [5]. We use the parameters
used in [5]:0.35⁄

1⁄200 ,⁄

0.65 1

⁄

and

⁄

4,
1.05

2 ⁄

1,

0.8

0.2exp2 ⁄,. Fourier modes included in thefollowing calculations are ones having resonant surfaces between 0.2⁄

0.8 withcold pulse

w/o cold pulse⁄⁄

Fig.1: Temporal evolutions of ion temperature fluctuation energy.
∆

4,

100 and nonresonant ones with

0 and | |20. Number of radial grids

is 512. The code is well parallelized using Message Passing Interface (MPI) and OpenMP.
Figure 1 shows temporal evolutions of ion temperature fluctuation energy. As shown in theleft panel, the 1, 0 component increases slowly and becomes larger than the turbulencefluctuations with

0, when no cold pulse is imposed.

The 1, 0 component seems to be due to the GAMs driven
by the turbulence. The cold pulse puts ion temperaturefluctuation energy into

, 0 components directly as

shown in the right panel in which the cold pulse is imposedin the hatched region ( 400 ≦ ≦ 600 ). The

0turbulent fluctuations also increase during the cold pulse.
This is due to profile steepening of ion temperature by thecold pulse. Excitation of the GAM oscillation of ZFs isobserved around ⁄ ~0.6

Fig.2: Zonal flow behavior
0.8 from ~450 a little after the start of the cold pulse. It is

supposed that direct excitation of GAMs by the cold pulse is not so strong, rather they are
excited by the turbulence strengthened.
The computations were carried out on the HELIOS system at the IFERC-CSC, Rokkasho,
Aomori, Japan.
[1] J. D. Callen and M. W. Kissick, Plasma Phys. Control. Fusion 39 (1997) B173.
[2] M. Yagi, et al., “Simulation study of nonlocal transport from edge to core in tokamak plasmas,” to be
published in Contrib. Plasma Phys.
[3] M. Yagi, et al., “Simulation study on nonlocal transport for peripheral density source,” 23rd
International Toki Conference (Toki, Japan, 21 Nov. 2013) P1-43.
[4] N. Miyato, Y. Kishimoto and J. Q. Li, Phys. Plasmas 11 (2004) 5557.
[5] N. Miyato, Y. Kishimoto and J. Q. Li, Plasma Phys. Control. Fusion 48, A335 (2006).

Simulation study of nonlocal transport from edge to core in tokamak plasmas
M. Yagi, M. Miyato, A. Matsuyama and T. Takizuka
Japan Atomic Energy Agency, Rokkasho, Aomori
Graduate School of Engineering, Osaka University, Suita, Osaka
The transient transport events are observed in the magnetic confinement devices[1,2]. The
time scale of the transient transport is much faster than the diffusive time scale, which isdifficult to explain based on the local transport model. Theoretically, the non-local transportmodel was proposed based on the one-dimensional integral heat flux model[3]. In this work,
the non-local transport for peripheral density source is investigated by using 4-field reduced
MHD model[4]. The spherical source is introduced in the density evolution equation# r 2 + rs2 − 2rrs cosθ + ζ 2 / ε 2 &S = SAMP exp % −
(
2Δ 2
$
'

(1)

where SAMP = 10 −2 , rs = 0.8 , Δ = 0.1 , ε = 1 / 3 are used. The spherical source is applied at
t=960 after saturation of the resistive ballooning turbulence, which is excited in the peripheral
region. Then, the source is switched off at t=1800. Figure 1 shows the time evolution ofinternal energy for each Fourier mode, which is defined by
Em,n =

1 1∫ r dr pm,n (r)
2 0

2

(2)

where m and n represent the poloidal and toroidal mode numbers. It is shown that the energy
is directly transferred into (0,0) and (1,0) modes as well as finite n modes from the source at
T=960. Figure 2 shows the time evolution of flux averaged total density profile at T=960,
1800, 2500, which is calculated byptot (r) = peq (r) + p0,0 (r)(3)

It is shown that the nonlocal transport appears after switching-on the source at T=1800. The q-profile is also shown in Fig.2, where q indicates the safety factor. Although the nonlocal
transport appears in the vicinity of q = 3/2 surface, it is not an exact location of q = 3/2
surface.
The simulation result is summarized as follows: (1) When the spherical source is applied in
the plasma edge, the energy is directly transferred into p0,0 and p1,0 modes as well as finite nmodes. (2) When the source is switched on, the dipole structure starts to deform and the spiralstructure is formed by p1,0 modes, which connects the core and the edge regions. (3) p1,0
modes interact with p0,0 via the three wave coupling and the toroidal coupling, whichproduces the nonlocal transport in the vicinity of q = 3/2 surface. We confirm this mechanism
not only by simulations without the convective nonlinearity in p0,0 evolution equation or
without the toroidal coupling (in the cylinder limit). Comparing with the cylindrical source, it
is found that the mechanism to produce the nonlocal transport is the same, however, it appears
after switching-on the source for spherical source, on the other hand, it appears just afterswitching-off it for cylindrical source[5]. The nonlocal transport obtained in this simulation isdriven by the macroscopic structure, which is consistent with experimental observation[6].
epre( 0/ 0)epre( 2/ -1)
epre( 3/ -1)
epre( 4/ -2)
epre( 22/ -6)
epre( 35/ -9)epre( 19/ -4)
epre( 1/ 0)

10-3
10-4
10-5
E

m,n

3 10-2

4
T=1800

2 10-2

10-6

3

P

q

T=960

tot

-7

10

T=1800

1 10-210-8

2T=2500

10-9
10-10

0

500

1000

1500

2000

2500

0
0.2

0.3

0.4

0.5

t

0.6

0.7

0.8

1
0.9

r

Fig.1 Time evolution of internal energy for eachFig.2 Time evolution of flux averaged

Fourier mode. The source is applied at t=960 andtotal density profile at t=960(green),is switched off at t=1800.

1800(red), 2500(blue). q-profile at
t=1800 is also plotted(black).

This work is partially supported by Grant-in-Aid for Scientific Research (No. 23246163). The
simulations are performed by using HELIOS super computer in IFERC-CSC, Aomori, Japan.
The Collaborative Research Program of Research Institute for Applied Mechanics, Kyushu
University is also acknowledged.
[1] U. Stroth, L Giannone, H-J Hartfuss, et al., Plasma Phys. Control. Fusion 38, 1087 (1996).
[2] N. Tamura, S. Inagaki, K. Ida, et al., Phys. Plasmas 12, 110705 (2005).
[3] T. Iwasaki, S.-I. Itoh, M. Yagi, et al., J. Phys. Soc. Jpn. 68, 478 (1999).
[4] M. Yagi, S. Yoshida, S.-I. Itoh, et al., Nucl. Fusion 45, 900 (2005).
[5] M. Yagi, A. Matsuyama, N. Miyato and T. Takizuka, 14th International Workshop on
Plasma Edge Theory in Fusion Devices (Cracow Poland, 23-25 Sep. 2013).
[6] S. Inagaki, T. Tokuzawa, K. Itoh, et al., Phys. Rev. Lett. 107,115001 (2011).

Integrated Heat Transport Simulation of High Ion Temperature
Plasma of LHD
S. Murakami, H. Yamaguchi, A. Sakai, N. Nagaoka∗ , H. Takahashi∗ , H. Nakano∗ , M. Osakabe∗ ,K. Ida∗ , M. Yoshinuma∗ , M. Yokoyama∗ A. Wakasa, A. Fukuyama and LHD Exp Group∗
Department of Nuclear Engineering, Kyoto Univ., Kyoto 606-8540, Japan∗ National Institute for Fusion Science, 322-6 Oroshi, Toki 509-5292, Japane-mail: murakami@nucleng.kyoto-u.ac.jp
High ion temperature plasma with carbon pellet injection of LHD is investigated by the integrated
heat transport simulation. NBI heating deposition of time evolving plasma is evaluated by GNETTD and the heat transport of multi-ion species plasma (e, H, He, C) is studied by TASK3D. Achievement of high ion temperature plasma is attributed to the 1) increase of heating power per ion due
to the temporal increase of Z-eﬀ, 2) reduction of eﬀective neoclassical transport with impurities, 3)
reduction of turbulence transport. The reduction of turbulence transport is most signiﬁcant contribution to achieve the high ion temperature and the reduction of the turbulent transport fromthe L-mode plasma (normal hydrogen plasma) is evaluated to be factor more than three by using
integrated heat transport simulation code.

19-3[x10 m ]20(a) #110599
Tang. x3
Introduction
10Perp. x2
0High ion temperature, Ti , experiments have been per(b)2formed applying the tangential and perpendicular NBI heating
1
C pellet0systems in LHD[1]. The high Ti plasma up to 7.3keV has been
2
7
obtained during the decay phase of the density after rapid in6
crease due to a carbon pellet injection. Simple heat transport
5
4
analysis of these high Ti plasma has been done but the rapid3change of density eﬀect on the NBI heat deposition and the21
multi-ion species eﬀects on the heat transport is not treated0
accurately. In order to analyze the transport property of the-6 -4 -2
0
2
4
6
time evolving plasma, we have to use NBI heat depositionv||/vth0
analysis code including the eﬀect of plasma time evolution.
3.5
P (r/a<0.5)
P (r/a<0.5)
3
Also the plasma contains suﬃcient impurities due to the He
2.5
gas puﬀ and C pellet injection and the heat transport simula2
tion should take care the multi-ion species.
1.5
In this paper we study the high Ti plasma with C pel1
let injection of LHD applying the integrated simulation code
0.5
GNET-TD and TASK3D. We perform time-dependent NBI
0
4
4.2
4.4
4.6
4.8
5
heating simulation using GNET-TD code, which is a modiﬁed
Times[s]version of the 5D drift kinetic equation solver GNET[2]. Next
we investigate the heat transport of high Ti plasma assum- Figure 1: (Top) Time history of NBI
ing multi-ion species plasma (e, H, He, C) by the integrated heating power and the density, (center)NBI beam ion distribution at t=4.7s,transport code, TASK3D [3].
[MW]

e
i

(bottom) Time history of NBI heat deposition.

NBI heating analysis of C-pellet injection plasma
In order to analyze the NBI heat deposition proﬁle of time evolving plasma we have developed
GNET-TD including the ﬁnite drift orbit and complex motion of trapped particles. We performed
NBI heating simulation of high-Ti discharge (# 110599) plasma, in which the plasma density changes
in time. In the simulation time change of densities and temperature proﬁles are from the experimental measurements. For simplicity we assumed T e = T i because Vbeam >> Vi . There are ﬁve NBI
injectors and one of them #4 is modulated in order to measure the ion temperature (Fig. 1-(top)).
We evaluate the beam ion distributions including time evolutions of density, temperature andbeam power modulations. Before the C pellet injection beam ion distribution reaches a steadystate slowing down distribution, and just after the pellet injection strong slowing down of energetic
1

beam can be seen at the edge region. But we can not see clear change in the central region. Afterthe density decay phase the beam distribution is recovered. Figure 1-(center) shows the beam ion
distribution in velocity space after the pellet injection in the core region (r/a < 0.5). We can see
the slowing down distribution of co and counter injected and perpendicular injected NBI beam ions.
Figure 1-(bottom) shows the NBI heat deposition in the core region (r/a < 0.5). The heat deposition
is jump up by the pellet injection and the start of the perpendicular #4 injector at t = 4.575s.
Heat transport simulation of high-Ti plasma
In order to perform the integrated heat transport simulation of high Ti plasma we apply the
improved TASK3D code assuming multi-ion species plasma (e, H, He, C). In TASK3D one dimensional diﬀusive transport is solved in the TR module. The heat source is loaded from the resultsof GNET-TD.We assume that the heat transport consists of the neoclassical and the turbulent
ones. In TASK3D the neoclassical transport database, DGN/LHD, evaluates the neoclassical heattransport for all species and the radial electric ﬁeld is determine by the am-bipolar condition ofneoclassical electron and ion ﬂuxes. For the turbulent transport we assume the gyro-Bohm( model
)( )
ρ
T
T
B
and the gyro-Bohm grad-Ti model for the electron and ion heat conductions as χe = Ce eB
L ,
(

)( )(

0

)

ρ
T
aT
χTi B = Ci eB
L
T , where Ce and Ci are the constant factors estimated by best ﬁtting of thehydrogen L-mode (normal) shot proﬁle data[3].
We perform the integrated heat transport simulation of highTi discharge (# 110599) plasma using the density proﬁles from(a)the experimental measurements. The electron temperature drops
rapidly by the pellet injection and backed to the previous temperature. This shows similar behavior with the experimental one.
The Ti increases about 0.5keV just after the pellet injection until
t=4.7s and gradually decreases. Also we found that the incrementof Ti after the pellet injection, ∆Ti , increases more if we change
the carbon density from the experimentally obtained values. Thisis due to the increase of heating power per ion by the temporaldecrease of the ion number density and, also, due to the reductionof eﬀective neoclassical transport with impurities by the C-pellet
injection. However the increment of ion temperature is lower than(b)that of experimental one (∆Ti ' 3.5keV). This indicates that the
reduction of the turbulent transport is necessary to explain theobserved high-Ti .
Next we introduce the conﬁnement improvement factor,
γT B , inheat transport model as χTi B =( the) ( ion( turbulent
)
)
0
ρ
T
aT
1
γT B Ci eBL
T . Also we introduce the carbon density factor, fC , as nCsim = fC × nCexp . Figure 2-(a) shows ∆Ti , as afunction of γT B with diﬀerent carbon densities. It is found thatabout 2keV of ∆Ti is obtained when we set γT B = 2.5 and that ∆Ti[s]increases about 1keV when fC = 2.0. Figure 2-(b) shows the time
evolution of ∆Ti changing the carbon density with γT B = 2.5. We
can see a rapid increase of ∆Ti and, then, the decrease gradually. Figure 2: (a) Increment of ion
temperature after C-pellet injecThis increase of ∆Ti is mainly due to the improvement of turbulent tion, ∆T , as a function of the imi
transport. It is also found that the increment of ∆Ti due to the proving factor, γT B with diﬀerent
increase of the carbon density.
carbon densities. (b) Time history[1] H. Takahashi et al., Nucl. Fusion 53 (2013) 073034.
[2] S. Murakami et al., Nucl. Fusion 46 (2006) S425.
[3] A. Wakasa, et al., Proc. 39th EPS Conf. and 16th Int. Cong.
Plasma Phys., 2012, P2.028.

2of ∆Ti changing the diﬀerent carbon densities with γT B =2.5.

Modeling of Heat Diffusivity for Ion Temperature Gradient Turbulencein Helical Plasmas
S. Toda, M. Nunami, A. Ishizawa, T. -H. Watanabe, H. Sugama

National Institute for Fusion Science, 322-6 Oroshi-cho, Toki 509-5292, Japan

a(ρ)

R/L

Tc

Turbulent transport is one of the most critical issues for plasma confinement in magnetic
fusion devices. This is because the turbulent transport induces a large amount of the particle
and heat loss in toroidal plasmas. Further study of the plasma turbulence and the transport isnecessary for the improvement of the plasma performance. In the integrated transport analysis
suite, TASK3D, which is the one-dimensional (the radial) transport code, the gyro-Bohm
scaling on the heat turbulent diffusivities is used in LHD [1]. In the case of ITG turbulence,
the form χ i ~ρ2tiv ti f(γ, ŝ, τ ZF )/R is taken for the ion heat diffusivity [2,3], where f is somefunction of γ, ŝ and τZF. Here, γ is the growth rate of the mode energy, ŝ is the magnetic
shear parameter and τZF is so-called the decay time of zonal flows [4]. We need to study thatwhich modes are destabilized to determine which models should be chosen for the turbulent
diffusivities in transport codes. The GKV-X code solves the gyro-kinetic equation to examine
Ion Temperature Gradient (ITG) modes and zonal flows in Large Helical Device (LHD) for
the turbulent transport in helical plasmas [2]. For the first step, we perform the linear
gyro-kinetic analysis by using GKV-X codes for ITG modes in the LHD high-Ti dischargecase of the shot number 88343, even if GKV-X codes retain the nonlinear characteristic. Thisis because the simulation cost of the linear analysis is extremely smaller than that of thenonlinear simulation [5]. The turbulent diffusivity is roughly determined by the linear growthrate and the wavelength according to the mixing length estimate in the gyro-kinetic simulation
results. The calculation of the growth rate by use of the gyro-kinetic codes at each time step inthe dynamical transport codes takes a high cost. Modeling of the ion heat diffusivity isnecessary in terms of the significant parameter for the plasma instability, which is available inthe dynamical transport codes. In this study, the significant parameter for the ITG instability
is considered to be L Ti (= −Ti / Ti ') . The other parameters such
20.0
(a)as Ln and the safety factor are fixed. The dependence of thezonal flow decay time τZF on the function f is neglected in
15.0
modeling χ i . The form for the ion heat diffusivity χi is derived.
This form is used in the transport code to compare the
10.0
simulation results with the radial profile for the ion temperature
in LHD experimental results.
5.0Firstly, the simulation of the linear analysis is done using
GKV-X codes. The ITG instability is examined in the LHD
0.0
0.0
0.2
0.4
0.6
0.8
1.0
ρ
high-Ti discharge case of the shot number 88343. This
0.7
simulation is for the three dimensional equilibrium field
0.6
configuration with the major radius R=3.75m. We set Te=Ti.
(b)
0.5
The saturation level of the linear growth rate for the
perturbation of the electrostatic potential is calculated. The
0.4
value of the ion heat diffusivity χ i / χGB is evaluated as
0.3
χ i / χ GB = C ∑ k γ / k 2y , where χ GB (= ρ2ti v ti / R) is the gyro-Bohm
0.2
diffusivity, ky is the poloidal wavenumber and γ is the linear
0.1
growth rate for the ITG mode. We use the numerical factor C
0.0
0.0
0.2
0.4
0.6
0.8
1.0
to match the simulation results with the experimental results
ρand set C=1. To find the critical ion temperature gradient for
1 The radial dependence
the ITG mode, the dependence of χ i / χGB on the normalized ion Fig.
of (a) the critical ion
temperature gradient R / L Ti is examined. The dependence of temperature gradient R/LTc andχ i / χ GB on R / L Ti is examined with all plasma parameters except (b) the slope a(ρ) in modeledyion diffusivity

2

diffusivity(m /s)i

T (KeV)the temperature gradient fixed at the stationary values for the density and the safety factor at
fifty radial points. We model χ i / χGB as χ i / χGB = a(ρ)(R / L Ti − R / L Tc ) , where L Tc is thenormalized critical ion temperature gradient for the ITG instability and ρ = r / a . The slopea(ρ) depends on the values of the normalized density gradient and the safety factor, whichchange due to the radial positions. The critical ion temperature gradient for the ITG mode,
R/LTc and the slope a(ρ) in terms of R/LTi are shown in Fig. 1.
Secondly, the transport dynamics is examined using the modeled ion turbulent diffusivity
above, when the transport code [6] is performed. The neoclassical diffusion coefficient isderived from DGN/LHD database in the case of low- β limit. The radial profile of the electricfield is derived from the ambipolar condition. The three solutions of the ambipolar radial
electric field are found in the radial region 0.25 < ρ < 0.80 . The#88343,2t=2.23s,2(almost2stationary)24
LHD2experimental2results (a)positive radial electric field is chosen from three solutions of
the ambipolar conditions in the radial region 0.25 < ρ < 0.80 .
3
The plasma radial profiles except the ion temperature, such as2
the density, the electron temperature, the radial electric field
and the safety factor are set to temporally be constant. The
1experimental result in #88343 at t=2.230s of LHD is used for
simulation2results2(stationary)
the radial profiles of the density, the electron temperature, and0
the safety factor. The dynamics of the Ti radial profile is
0.0
0.2
0.4
0.6
0.8
1.0
ρ
simulated, by use of the transport code. The Ti profile of the
6.0(b)experimental results is used as an initial state for the
5.0
dynamical transport simulation. We show the simulationχ
4.0
results for the stationary ion temperature profile with the solidline in Fig. 2(a). The dashed line indicates the radial profile of
3.0
Ti in #88343 at t=2.230s of LHD. The simulation results forχ
2.0
the radial Ti profile can show the good agreement with one
1.0
experimental result in #88343 at t=2.230s of LHD. The
stationary radial profiles of the turbulent and neoclassical
0.0
0.0
0.2
0.4
0.6
0.8
1.0
ρ
diffusivities are also shown in Fig. 2(b). The ITG mode isFig. 2 Radial profiles of (a)stabilized in two radial regions 0.0 < ρ < 0.25 and 0.8 < ρ < 0.92 .
the ion temperature and (b)the ion diffusivity
The solid and dashed lines represent ion turbulent andneoclassical diffusivities. The turbulent transport is dominantcompared with the neoclassical transport in the radial region where the positive electric field
is chosen.
In this study, modeling of the ion heat diffusivity is done in terms of the scale length of theion temperature gradient to introduce the ITG turbulence in the transport codes, especially in
the LHD high-Ti discharge case of #88343. Using the dynamical transport codes, we
compared the simulation results using modeled ion turbulent heat diffusivity with anexperimental result in the LHD high-Ti case of #88343. We should compare the simulation
results of the linear gyro-kinetic analysis with the results by use of the nonlinear analysis.
Furthermore, the effect of zonal flows is necessary to be included. Modeling of the heatturbulent diffusivity will be done in more general cases. These are for the future studies.
NEO

i

i

References[1] A. Wakasa et al., 39th EPS Conference on Plasma Physics joint with the 16th International
Congress on Plasma Physics, P2.028 (2012)[2] M. Nunami et al., Physics of Plasmas 19, 042504 (2012)[3] J. W. Connor and H. R. Wilson, Plasma Phys. Control. Fusion 36, 719 (1994)[4] H. Sugama and T. -H. Watanabe, Physics of Plasmas 13, 012501 (2006)[5] M. Nunami et al., Physics of Plasmas 20, 092307 (2013)[6] S. Toda and K. Itoh, Plasma Phys. Contr. Fusion 44, 325 (2002)

Numerical analysis of neoclassical flow by the moment method at the
NBI Plasma in Heliotron J
K. Nishioka1, Y. Nakamura1, S. Nishimura2, H. Y. Lee3, S. Kobayashi4, T. Mizuuchi4, K. Nagasaki4,
H. Okada4, T. Minami4, S. Kado4, S.Yamamoto4, S. Ohshima4,S. Konoshima4 and F. Sano4
1

Graduate School of Energy Science, Kyoto University 2 National Institute for Fusion Science3

Natural Science Research Institute, Korea Advanced Institute of Science and Technology,

4

Institute of Advanced Energy, Kyoto University

The investigation of the plasma flow and viscosity is important to understand the transport and stabilityof magnetically confined plasma. For the experimental effort, carbon parallel flow measurement has been
carried out in many tokamak and stellarator/Heliotron devices, also in Heliotron J [1], based on the CXRS
method. However, it is difficult for CXRS to measure the bulk ion flow velocity directly. Therefore,
numerical method has important roles for the investigation of parallel flow and viscosity.
When we targeted the NBI plasma, we should take into account the existence of external momentum.
When high energy neutrals are injected into plasma, parallel momentum was also injected. This externalsource affects the neoclassical transport matrix itself through the change in parallel momentum balance.
When fast ion exists, the general force balance equation is as follows [2,3];

 M a   0    00Mb
00   Λ aa
0   Λ ba

 
 
M N   Λ Na
Na
0
=


0

Λ ab
Λ bb
Λ Nb
0
Nb
0

Λ aN    U a Λ bN    U b     
Λ NN    U N 
0   Xa   Z a  Ca 





C 
0   Xb   Zb BE   b  BFf 1
          
N N   XN  ZN 
C N 

where M is  jmax  1   jmax  1 viscosity coefficient matrix, Λ is  jmax  1   jmax  1 friction
coefficient matrix, N is 2   jmax  1 matrix which indicates viscosity from radial potential gradient, X :drive force from pressure and temperature gradient, Z : momentum driven by parallel electric field, C :momentum transfer ratio to each particle specie and

BF f 1is the total parallel momentum input.

Thejmax is the order of Laguerre expantion. Our new term is the last term of right hand. We obtained parallel
particle and heat flow U a by solving above simulatenous equation considering the effect of external
momentum and finite radial electric field which satisfies ambipolar condition  ea a  0 .
aNumerical estimation by the moment method and CXRS measurement of the parallel C 6+ flow velocity
in the standard (STD) and high gamma (HG) field on Heliotron J are shown in Figure 1(a) and (b).

Fig. 1. Radial profiles of parallel C6+ flow velocity measured by CXRS (dotted-lines) and of calculated neoclassicalflow velocity (solid-lines) in STD configuration at (a) and HG configuration at (b).

These figures indicate that numerical estimation of carbon flow is consistent with the experimental result
at both co- and counter- NBI plasma. In addition to this, we observed that HG configuration reduces u|| bya factor of 1/2 that in STD configuration at central region both in measurement and calculation.
We show the parallel viscosity coefficient of C6+ ion both in STD and HG configurations in Fig.2. This

figure denotes that neoclassical C6+ viscosity coefficient in HG configuration is much larger than that in
STD configuration by the large parallel magnetic ripple depth. The detailed estimation of controllabilityof parallel plasma flow, and poloidal and toroidal viscosity of ions are the furure work.

Fig. 2. Radial profiles of parallel C6+ ion viscosity coefficients in STD and HG configurations

References[1] H. Lee, et al. Plasma Phys. Ctr. Fusion 55 (2013) 035012[2] N.Nakajima and M.Okamoto, J. Phys. Soc. Jpn. 61, (1992), 833[3] S. Nishimura, et al, Phys. Plasmas. 17, 082510 (2010) , Phys. Plasmas. 18, 069901 (2011)Integrated simulation of ELM control by Pellet Injection	
 N. Hayashi 1, N. Aiba 1, T. Takizuka 2, N. Oyama 1
1

JAEA, 2 Osaka Univ.

The energy loss caused by the edge localized mode (ELM) needs to be reduced for ITER operations
with ELMy H-mode plasmas. The pellet injection is considered as one possible method to reduce the ELM
energy loss by increasing the ELM frequency (so called ELM pacing). The ELM pacing requires thesignificant reduction of ELM energy loss, small impact on the plasma performance, and less particle fueling.
The reduction of ELM energy loss by the pellet injection for ELM pacing has been studied by an integrated
code suite TOPICS (it has been formerly referred to as TOPICS-IB [1]. Now it is simply referred to as
TOPICS, since the core transport module which was called TOPICS has evolved integrating other functions
and modules so as to represent the suite). TOPICS simulations showed that the energy absorption by pellet
ablation clouds and the subsequent ExB drift can steepen the pressure gradient locally in the pedestal, which
destabilizes high-n ballooning or intermediate-n peeling-ballooning modes and triggers an ELM, where n is
the toroidal mode number [1]. It was also found that the ELM energy loss can be significantly reduced by the
low-field-side injection of a pellet with about 1 % of pedestal particle content and a speed fast enough toapproach the pedestal top when the pedestal pressure is about 95 % of natural ELM onset [2]. In this paper,
we do simulations of JT-60U plasmas with various parameters and ITER plasmas of standard operation
scenario, examine the effectiveness of the above suitable conditions for ELM pacing and evaluate the pellet
particle content required for ELM pacing.
The pellet triggered ELM is investigated by TOPICS, which couples the core transport module withan MHD stability code MARG2D, a scrape-off-layer / divertor transport model D5PM and a pellet model
APLEX taking account of the ExB drift of pellet plasma cloud [1]. The core transport module is based on the
1.5D core transport code, which solves the 1D transport and current diffusion equations on the normalized
toroidal flux coordinate, and the Grad-Shafranov equation of MHD equilibrium on the 2D plane. Particle andthermal diffusivities consist of neoclassical and anomalous ones where the anomalous diffusivities are given
as empirical formulas. In order to produce the H-mode pedestal, the transport near the edge is reduced to theneoclassical value calculated by the matrix inversion method with a prescribed characteristic width ofpedestal. The ELM model is developed by coupling TOPICS with MARG2D. In the present simulation, thestability of n = 1-50 modes are examined at given time-intervals. When unstable modes appear, an ELM is
assumed to occur. ELM enhanced diffusivities are added on the basis of eigenfunction profiles of unstablemodes, where the maximum value is given as a parameter. The ELM enhanced transport is maintained for atime interval given as a parameter. D5PM is based on time-dependent fluid equations. Particle flux and heat
fluxes across the separatrix obtained in TOPICS are used as inputs and D5PM calculates the SOL density
and temperatures at the separatrix, which are used as boundary conditions in TOPICS. The particle balancein the whole regions is modeled by integrating neutral models with TOPICS and D5PM. APLEX assumes asphere pellet and simulates the pellet ablation, ExB drift of ablated and detached plasma cloud, cloud energy
absorption from the background plasma and deposition to the background plasma. The pellet model APLEX
takes into account 3D geometry, in which plasma clouds produced from a pellet drift along the major radiusand expand along magnetic field lines. It is worth noting that, during the pellet ablation and further drift, we
use the 1.5D transport code for the profile evolution of background plasma (1D transport and 2D MHD
equilibrium), and the linear MHD stability code MARG2D to check the stability of the 2D MHD equilibrium
and to model ELMs in the 1D transport. In our study, we fix the model parameters such as in the ELM modeland examine the dependence appearing in simulation results.
JT-60U parameters are the same as in previous studies [1,2], typical ones in ELMy H-mode
experiments. ITER parameters are for standard scenario and the ELM model parameters are the same as in
JT-60U simulations, except for four times larger diffusivities because the characteristic length, the pedestal
width, is roughly doubled in the real coordinate. The pedestal width almost agrees with the EPED prediction.
Figure 1 shows the ELM energy loss ΔWELM normalized by the pedestal stored energy Wped as a function ofpellet position at the ELM onset ρpel,onset for JT-60U and ITER simulations, where mode locations are also
shown. As the same as the previous result, with the suitable conditions, a pellet can significantly reduce theELM energy loss by penetrating deeply into the pedestal and triggering high-n ballooning modes localizednear the pedestal top. Without the suitable conditions, the pellet triggers lower-n modes with wide mode
width before the pellet approaches the pedestal top, and enhances the ELM energy loss, which becomescomparable to the natural ELM level at the maximum. In all simulations with the suitable conditions, asimilar level of reduced energy loss is obtained. Thus, the effectiveness of suitable conditions has been
confirmed. Figure 2 shows the relation between the pellet particle content Np (normalized by the pedestal
particle content Nped) required for ELM pacing and the variation of normalized pressure gradient Δα(normalized by the critical gradient αcr of ideal ballooning mode (IBM)) required to destabilize the high-nballooning modes. Here, the IBM stability can be used to judge the onset of high-n ballooning modes [3].
The pellet particle content required for ELM pacing is larger for the pedestal plasma with higher density andfarther from the IBM stability boundary near the pedestal top. For ITER, the required pellet particle content
is within the range of present empirically-designed value for the pellet to reach the pedestal top [4] and thus
our evaluation gives the physical background to the design.
References[1] N. Hayashi et al, Nucl. Fusion 51 (2011) 103030.
[2] N. Hayashi et al, Nucl. Fusion 53 (2013) 123009.
[3] N. Hayashi et al, Nucl. Fusion 52 (2012) 114002.
[4] A. Loarte et al, 24th IAEA FEC, San Diego, 2012

FIG.1. Dependence of ΔWELM/Wped on pellet position at ELM onset for JT-60U simulations (blue open) andan ITER simulation (red closed), where horizontal solid lines denote mode locations, circle symbols do
results with suitable conditions of pellet injection for ELM pacing, a square symbol does an example without
suitable conditions and shaded regions denote other results in JT-60U. Vertical and horizontal dotted lines
denote pedestal top position and ΔWELM of natural ELM level, respectively.

FIG.2. Dependence of Np/Nped on Δα/αcr for JT-60U (blue open) and ITER (red closed) simulations, whereshade region denote pacing pellet size of present ITER design [3]. In ITER simulations, normalized internal
inductance li is varied from a typical value to ITER lower limit (A→B).

ペレット入射プラズモイドのドリフト機構のモデリングと実験との比較	
 Study on Drift Displacement of Pellet Plasmoid and Comparison to Experiments
松山顕之
Akinobu MATSUYAMA
日本原子力研究開発機構
Japan Atomic Energy Agency
1. Introduction
Historically, an issue for pellet fuelling was to develop high-speed injectors that are capable of injectingpellets successively with the speed of 5-10 km/s. At present, a more favorable scheme is available with thehigh-field side (HFS) pellet injection. Injecting the pellet from the HFS was found to improve the fuelingefficiency as compared to a conventional low-field side (LFS) injection scheme. In such an injection
configuration, a plasma cloud surrounding the pellet drifts spontaneously towards the plasma core, andprovides deep deposition of pellet particles without fueling degradation due to high additional heating power.
Therefore, the HFS injection is now considered to be the most promising technique for fuelling tokamak
reactors, allowing use with the injection velocity of < 1 km/s, which is realized by the present technology.
To predict the pellet fueling properties in present tokamaks and future reactors, simulation models forthe pellet injection especially taking into account the drift
displacement of particle deposition have been developed. Using
the HPI2 code [1], the simulation of pellet particle deposition
has been compared to experimental results in present tokamaks,
such as JET, Tore-Supra, and DIII-D. However, the validationof existing drift models against experiments is still underdiscussion. As an attempt for extending such validation study,
we here investigate the pellet particle deposition in magnetic
configurations other than tokamaks such as RFPs and helicaldevices. This is a good test for validating if present drift models
can reproduce the drift displacement taking into accountdifferences of the magnetic geometries in these devices. For
this purpose, we extend the drift models implemented in the
HPI2 code and compare the simulation data against the pellet
Fig. 1 	  Schematic picture of theinjection experiments in the Large Helical Device (LHD) [2].
suppression of polarization current byplasmoid internal current in a straight
2. Drift Modeling of Pellet Plasmoid
Helical devices possess three-dimensional magnetic geometries, helical geometry.
which exhibit intrinsically different characteristics from that ofaxisymmetric tokamaks. We consider two important characteristics that affect the drift displacement: (1)steep magnetic field gradient that varies with the characteristic scale of a plasma minor radius, depending on
the distance from the external helical coils, and (2) a helical modulation of the magnetic field with shortconnection lengths, where the direction of the magnetic field gradient varies rapidly when plasmoids expandin the toroidal direction. Because of the property (1), the plasmoid drift acceleration at the pellet injection
location is much stronger than what is expected for tokamaks by up to a factor of the aspect ratio.
Nevertheless, the plasmoid internal current suppresses this strong drift displacement as soon as the plasmoid
expands over the connection length. The latter was originally considered for tokamak geometry [3], which
we extend to general magnetic configurations and implement in the HPI2 code. Figure 1 shows a schematic
picture of this effect in a straight helical configuration. At the pellet injection point on the low field side, themagnetic field gradient directs inwardly, but after the plasmoid expends at the sonic velocity over a
half-period of the helical magnetic field, the gradient at the end of the plasmoid is seen to direct outwardly.
For the divergence-free condition of the plasmoid internal current to be satisfied, the field-aligned current
connects these two opposite currents. As a result, field-aligned currents suppress the drift acceleration
significantly, which is understood in the analogy to the Pfirsch-Schlüter current in the MHD equilibriumtheory. We found that the extended HPI2 code taking into account this drift suppression effect well
reproduces the drift velocity and experimental deposition in the LHD [2].

3. High-Field Side Injection in Large Helical Device
A clear demonstration of the effect of different magnetic geometries between tokamaks and helical devices
has been observed experimentally for the HFS injection. Surprisingly, the HFS injection in the LHD did notexhibit the improved capability of fuelling, and even the outward drift displacement was observed [4]. Our
main result is that this unexpected result for the HFS injection can be understood from the effect offield-aligned current discussed in the previous section, if one takes into account realistic magnetic field
geometry including the magnetic hill in the outer plasma region of the LHD. Figure 2 shows the effect ofmagnetic well/hill structure on the particle deposition, simulated by the extended HPI2 model. If we change
the magnetic-field configuration from that with the magnetic well to that with the magnetic hill, it can beseen that for the HFS injection, the outward displacement as compared to the ablation profile is obtained.
Even when one injects the pellet from the HFS and the local drift velocity directs inwardly, the drift
suppression due to the field-aligned current occurs. After the compensation of a local drift velocity, the drift
direction is dominated by whether the magnetic field in the deposition region possesses magnetic well or hill
structure. This is because the drift due to the averaged magnetic curvature is unaffected by the field-aligned
current. We mention that the drift displacement due to the magnetic well discussed here is correlated to aparticular importance of the interchange mode in heliotron devices. As a whole, the LHD results are seen to
be in contrast to what is expected for tokamaks, where the direction of drift displacement is dominatedmainly by a local magnetic curvature(ballooning-likedriftmotion).
The
connection length should dominate a relative
importance between such interchange andballooning drift motion; this idea extends theoriginal model discussed by Rozhansky, et al.
[3], and is useful for understanding the drift
displacementin
different
magnetic
configurations. We applied the extended
HPI2 code for experimental parameters of
Fig. 2	  Effect of edge magnetic field structure on the[4] and found that the simulation data well
deposition profile: (1) LFS injection case and (2) HFS
reproduces observed outward displacement
injection case when the magnetic field configuration iswith respect to both the peak location and thevaried from that with the magnetic hill to that with thedensity increase.
magnetic well.

4. Conclusion
In this work, the pellet fueling models are extended for understanding the drift displacement in different
magnetic configurations, i.e., tokamaks, RFPs and helical devices. The extended HPI2 code was comparedwith the experimental results in the LHD. An essential point is to include the effect of field-aligned current
that compensates the local drift acceleration, which yields the long-term drift depending on whether themagnetic configuration possesses magnetic well or hill structure. The developed code has been confirmed,
for the first time, to well reproduce main characteristics of the pellet fueling in the LHD.
Further code validation will be important, in particular, concerning the effect of low-order rationalsurfaces [5]. In the previous work, the correlation between the deposition peak and the position of low-orderrational surfaces has been reported for tokamaks [6], and it will affect quantitative prediction of the fuelingthe reactors. For this reason, testing this effect in the LHD is of particular interest, where the strong drift
damping due to the q = 1 surface at the deposition region is suggested (q: safety factor).

References
[1]
[2]
[3]
[4]
[5]
[6]

F. Koechl, et al., Submitted to Nucl. Fusion.
A. Matsuyama, et al., Nucl. Fusion 52, 123017 (2012).
V. Rozhansky, et al., Plasma Phys. Control. Fusion 46, 575 (2004).
H. Yamada, et al., Plasma Phys. Control. Fusion 43, A55 (2001).
B. Pégourié and M. A. Dubois, Nucl. Fusion, 29, 745 (1989).
N. Commaux, et al., Nucl. Fusion 50, 025011 (2010).

11th Burning Plasma Simulation Initiative (BPSI) Meeting
December 19-20, 2013 (Kyushu University)

Radio-frequency sheath-plasma interactions in the ion cyclotron range of frequencies
Haruhiko Kohno1,*, James R. Myra2 and Daniel A. D’Ippolito2
1

Department of Mechanical Information Science and Technology, Kyushu Institute of Technology
680-4 Kawazu, Iizuka, Fukuoka 820-8502, Japan
2
Lodestar Research Corporation, 2400 Central Avenue P-5, Boulder, Colorado 80301, USA*Author to whom correspondence should be addressed E-mail: kohno@mse.kyutech.ac.jp

ABSTRACTA numerical code that solves self-consistent radio-frequency (RF) sheath-plasma interactions in the ion cyclotron
range of frequencies (ICRF) is developed based on a finite element technique. The present code solves for waves in the edge
plasma subject to a sheath boundary condition. In this two-dimensional study, it is found that rapid variation of conditions
along the sheath surface promotes coupling of the incident RF branch to a short-scale-length sheath-plasma wave (SPW).
The SPW propagates along the sheath surface in a particular direction dictated by the orientation of the magnetic field withrespect to the surface, and the wave energy in the SPW accumulates near places where the background magnetic field istangent to the surface. This phenomenon is expected to be relevant to the formation of localized hot spots on the bounding
surfaces, which are commonly observed on powered ICRF antennas.
Keywords: Radio-frequency sheaths, Plasma waves, Magnetic confinement fusion, Finite element method
1. INTRODUCTION

2. MATHEMATICAL MODELINGRadio-frequency (RF) waves in the ion cyclotron
range of frequencies (ICRF) have been successfully
applied to various heating experiments in fusion plasmas.
However, it has been revealed by experimental andtheoretical studies that deleterious edge plasmainteractions can occur due to various nonlinear
mechanisms [1]. One of the most important nonlinear
effects is RF sheath formation, in which the sheathpotential on the walls and limiters of the tokamak device
is enhanced by the ICRF waves. RF sheath formation
causes various problems including impurity generation
by enhanced sputtering and power dissipation in the
scrape-off layer (SOL), leading to local hot spots andreduced power. Because sheaths occur on the scale of the
Debye length, which is much smaller than domain sizes
or RF wavelengths of interest for fusion applications,approximate RF sheath boundary condition (sheath BC)
approaches have been developed [2,3]. These approaches
effectively treat the sheath as a thin vacuum layer, andthus capture the important effect of sheath capacitance on
the boundary. In addition, a finite element procedure
employing this sheath BC has recently been developedby the authors for the purpose of solving self-consistent
RF sheath-plasma interactions in the ICRF [4]. In thisshort paper, we will briefly summarize our recent findingof an important role in governing the distribution of wave
energy and induced potentials along the sheath surface
subject to rapid variation of conditions [5].

The governing equation for plasma waves in the
SOL is a combined form of Maxwell’s equations
described asÑ´Ñ´E -

w2
ε × E - iwm 0 J ext = 0,c2(1)where the electric field E and the external current J ext
vary on the RF time scale. Here, w is the applied ICRF
wave frequency, c is the speed of light, and m0 is thepermeability in vacuum. The dielectric tensor ε is given
based on the cold plasma model. At the metal wall, thesheath effect is taken into account by means of a sheathBC, which is written as follows:æDö
E t = Ñ t çç sh Dn ÷÷.
(2)
e
è sh
ø
Here, Dsh is the time-averaged sheath width, e sh is thedielectric constant in the sheath, Dn is the component of
the electric displacement normal to the sheath, and thesubscript t denotes the two components tangential to theboundary. For the linear case, the sheath width isexpressed as follows:
D sh = Cth lDe ,(3)

where lDe is the electron Debye length. The coefficient
Cth is a function of the surface contact angle of themagnetic field line (see [5] for more detail).

11th Burning Plasma Simulation Initiative (BPSI) Meeting
December 19-20, 2013 (Kyushu University)
3. RESULTS AND DISCUSSION
In the present finite-element numerical analysis,
we consider solutions of the sheath-plasma interactions
for spatially varying background magnetic fields whose
flux surfaces have tangency points at a boundary. Thismagnetic field configuration actually happens withlimiter surfaces in a tokamak. Figure 1 shows the filled
contour plot of E|| , which is the electric field component
parallel to the magnetic field line. Here, the antenna andmagnetic field lines are also superimposed on the plot.
Since this is a plot of E|| , the fast wave contribution isconsidered to be negligible. It is observed that the slow
wave (SW) propagates along the magnetic field lines
although it is evanescent from the antenna position. In
addition, the sheath-plasma wave (SPW) is generated on
the sheath surface between the “convex” flux-surface
tangency point (y = 0.3 m) and the intersection point (y =0.6 m) of a magnetic field line penetrating the antenna
and the wall (hereafter, the flux-surface tangency points
from which the magnetic field lines inside the plasma areseen as convex and concave curves are simply called theconvex and concave tangency points at the sheath surface,respectively). While the present geometry is too simple
to apply directly to a tokamak, in general the complex
magnetic fields of shaped diverted plasmas, and realistic
vessel walls can exhibit both types (convex and concave)of interactions. The generation of the SPW is also
confirmed in Fig. 2, which shows the variations of thereal part of the parallel electric field component along thesheath surface for two different boundary conditions at
the right wall; here, the result obtained by imposing the
conducting-wall BC is added for comparison. It is seenthat the SPW is generated only by imposing the sheath
BC. Furthermore, it is important to notice that the SPW isdriven in the narrow region within which the magnetic
field lines connect the antenna to the wall andaccumulation of its energy occurs at the convex tangency
point. It has been clarified through the analysis of theelectromagnetic SW-SPW dispersion relation that thewave mode propagates energy away from the concave
point and towards the convex point [5]. The dramatic
slowing of the group velocity near y = 0.3 m leads to an
accumulation of wave energy there.
In recent Alcator C-Mod experimental work, large
plasma potentials (~ 100 V) were found on field lines
connecting the emissive probe to a nearby plasma limiter[6]. This plasma potential rose rapidly in the vicinity of
the limiter tip and then decayed exponentially along the
major radius direction behind the limiter tip. The wave
energy deposition at the magnetic field tangency point on
the sheath surface, which has been discovered in thepresent work, may be the cause of these experimental
results. However, of course, further work is necessary toestablish this point.

FIG. 1. Filled contour plot showing the distribution of theelectric field component parallel to the magnetic fieldlines.

FIG. 2. Plot of the real part of the parallel electric fieldcomponent along the sheath surface with a comparison to
the result for the conducting-wall BC.
REFERENCES
[1]

[2]
[3]
[4]
[5]
[6]

J. R. Myra, D. A. D’Ippolito, D. A. Russell, L. A.
Berry, E. F. Jaeger, and M. D. Carter, Nucl. Fusion
46, S455 (2006).
J. R. Myra, D. A. D’Ippolito, and M. Bures, Phys.
Plasmas 1, 2890 (1994).
D. A. D’Ippolito and J. R. Myra, Phys. Plasmas 13,102508 (2006).
H. Kohno, J. R. Myra, and D. A. D’Ippolito,
Comput. Phys. Commun. 183, 2116 (2012).
H. Kohno, J. R. Myra, and D. A. D’Ippolito, Phys.
Plasmas 20, 082514 (2013).
R. Ochoukov, et al., J. Nucl. Mater. 438, S875
(2013).

One-Dimensional Transport Analysis of ECRH Assisted Start-up in JT-60SA
K. Hadaa, K. Nagasakib, K. Masudab, S. Idec and A. IsayamacaGraduate School of Energy Science, Kyoto University, Japanb
Institute of Advanced Energy, Kyoto University, Japanc
Japan Atomic Energy Agency, Naka, Ibaraki, Japan
E-mail address: hada.880@iae.kyoto-u.ac.jp

There is a considerable interest in plasma start-up assisted by electron cyclotron resonanceheating (ECRH) in superconducting tokamaks because available loop voltage is generally low
compared with conventional normal-conducting tokamaks. In JT-60SA superconducting
tokamak (R0 = 2.96 m, a = 1.18 m, and BT = 2.25 T), which is now under construction in
JAEA [1], the toroidal electric field for plasma breakdown is limited up to 0.5 Vm1. Thismay make it difficult to start up plasmas if the first wall is not in good condition. On JT-60SA,
a 7 MW ECRH system with multi frequency (110 GHz and 138 GHz) gyrotrons is equipped
for heating and current drive as well as for plasma start-up. The main purposes of this paper
are to examine the conditions for the reliable start-up and to understand the dominant physicalprocess in the plasma start-up phase. We have previously studied the plasma start-up assisted
by ECRH in JT-60SA using a zero-dimensional (0-D) model [2]. The calculation using
JT-60U parameters qualitatively reproduces the time evolution of the experimental results.
However, the 0-D model is not enough to investigate the plasma start-up assisted by 2nd
harmonic X-mode ECRH in large device such as JT-60SA. Hence, we are developing a
one-dimensional (1-D) model. In order to consider the change in ECRH absorption powerwith time evolution of electron temperature and density, we consider the ECRH absorption
efficiency. The 1-D model consists of energy transport equations for electrons and hydrogenions, particle transport equations for electrons and hydrogen atoms, and toroidal current
equation. In this model, we assume a cylindrically symmetrical plasma, a fixed plasmaboundary, and that peak of ECRH absorption is located on axis. Figure 1. shows that thethreshold ECRH power for plasma start-up is observed, the threshold ECRH power of around
1 MW is required to start up the plasmas for nH = 3.0 × 1018 m3, error field Berr of about 1
mT, and EC beam radius of aboutn (t = 0) = 1.0  10 [m ]5 cm. This indicates a promising1000
n (t = 0) = 2.0  10 [m ] B = 1.0 [mT]
outlook for a reliable start-up inn (t = 0) = 3.0 10 [m ] n /n = n /n = 0.1%
JT-60SA for a good wall condition.
800
n (t = 0) = 4.0 10 [m ]In case of low initial neutraln (t = 0) = 5.0 10 [m ]600
density, we can start up plasmas
without ECRH. However, another400
physical process is important in
low initial neutral density, that is,200
we have to consider the generation
100of runaway electron causing direct0loss. Figure 2. shows that the
0.0
0.5
1.0
1.5
2.0
2.5
3.0
threshold ECRH power increases
Injected ECRH power, PECRH, INJ [MW]with error field (Berr up to around7 mT). This results suggests thatFig 1. Dependence of electron temperature at  = 0.5
the optimized magnetic field and t = 0.15 s on initial neutral density and injected
control may be necessary for the ECRH power.
reliable start-up because the error
18

3

H

18

3

18

3

Te (= 0.5, t = 0.15 s) [eV]

H

err

C

H

18

3

18

3

H
H

e

O

e

Te (= 0.5, t = 0.15 s) [eV]800

Berr = 1.0 [mT]nH (t = 0) = 3.0  1018 [m3]Berr = 2.0 [mT]nC/ne = nO/ne = 0.1%

Berr = 3.0 [mT]Berr = 5.0 [mT]

600

Berr = 7.0 [mT]Berr = 10.0 [mT]

400
200
100
0
0.0

0.5

1.0

1.5

2.0

2.5

3.0

Injected ECRH power, PECRH, INJ [MW]Fig 2. Dependence of electron temperature at  = 0.5
and t = 0.15 s on error field and injected ECRH power.

4

OH start

3

2.0
ECRH start2ne (=0)

1.5
1.0

Te (=0)
Ip

1

0.5

0

0.0

T [keV], I [100 kA]

analysis provides an insight into
the dominant physical process asfollows. After the loop voltage isturned on, electron and ion
densities
gradually
increase;
therefore,
ionization
andequipartition
lossesbecomedominant. After ECRH injection,
electrondensitypromptly
increases, nevertheless electron
temperature still remains low.
When hydrogen atoms are fully
ionized,
ECRH
absorption
efficiency is improved andelectron temperature rises. This
leads to a positive feedback, as aresult,electrontemperature
drastically rises. The error field
loss becomes dominant due tohigh electron temperature until thepoloidal magnetic field generated
by the plasma current exceeds theerror field. This is considered to
be a reason for the ECRH power
threshold noted above to exist.

n [1018 m3]

evolution of electron density,
electron temperature, toroidal
current, and electron powersource and loss terms. This 1-D

1000

Power density (=0) [Wm3]

field generated by induced current
in passive structures such as
vacuum vessel, stabilizing plate,port, and cryostat is estimatedaround5-10
mT
without
optimized magnetic field control[3]. Figure 3. shows a time

106

104

ECRH
OH
IONIZ
EQU
ERROR
DRIFT

102

100
0.00

0.01

0.02

0.03

0.04

Time, t [s]Fig 3. Time evolution of electron density andtemperature, toroidal current, and electron powersource and loss terms.

[1] S. Ishida et al., Nucl. Fusion, 51 (2011) 094018[2] K. Hada et al., J. Plasma Fusion Res. 7, (2012) 2403104[3] M. Matsukawa et al, J. Plasma Fusion Res. 9, (2010) 264

Enhanced Neo-classical Resistivity Due to the m/n = 1/1 MHD-mode
Deformation of Central Core Plasma Column for Tokamak Hybrid Scenario
T. Takizuka
Graduate School of Engineering, Osaka University
Hybrid scenario is planned for ITER to test the various engineering components towards
DEMO with high reliability, high neutron fluence and long pulse length [1,2,3]. This scenario
is the combination of the inductive current drive and the non-inductive current drive, and theq profile is of the weak magnetic shear. The fraction of the non-inductive current including
the bootstrap current is expected to be 40~60%. Physical features of the hybrid scenario havebeen studied experimentally in the world. It was reported that the concentration of the current
density peculiar to the inductive discharges was much smaller than that of simple estimation.
DIII-D paper considered that this anomaly was caused by the hyper resistivity (see [4,5]) [6].
JT-60U paper showed similarly that the current concentration was smaller than the estimation[7]. Afterwards JT-60U report considered a possible reason, i.e., the enhanced electric
resistivity was increased in the central region by the peaking of effective charge Zeff with theaccumulation of tungsten impurity [8].
We recently presented a new model of the physical mechanism on the reduction of thecurrent concentration [9]. According to a JT-60U paper [7], there existed the n = 1 MHD
activity (different from the neo-classical tearing mode) in hybrid-scenario discharges. The
measured fluctuation of the electron temperature was large in the central region. This MHD
activity can be the internal kink mode with a poloidal mode number m = 1 and a toroidal mode
number n = 1, which becomes unstable for higher poloidal beta value (βp ~ 1) [10]. This modebrings the m/n = 1/1 kink deformation of the central plasma column between the magnetic axis
with q(0) < 1 and a magnetic surface with q ~ 1. We assume a simple form of the finite
normalized-amplitude of the deformation, ξ = ξmax – ρ/2 in a region 0 < ρ < ρb = 2ξmax, so thatmagnetic surfaces do not overlap. In the central core of a tokamak plasma with such an MHD
deformation, trapped particles can exist even on the magnetic axis. The effective inverse
aspect ratio becomes εeff = εa (ρ + ξ) at the normalized minor radius ρ, and the neo- classicalresistivity becomes as large as ηNC ≈ ηS /{1 – (2εeff)1/2} in the central region (εa : inverse
aspect ratio at the plasma surface, ηS : Spitzer resistivity).
In order to compare the present model with the experimental measurements, it is necessary
to take into account of the finite-deformation effect on the diagnostics of poloidal magnetic
field Bp. The measured Bp value oscillates in time due to the toroidal rotation of the mode. Itstime average <Bp> has been reported as the measured Bp, but it is smaller than the true Bp on amagnetic surface. Based on the result of E48158 at t = 27s in JT-60U paper [7], we revise themeasured data, poloidal magnetic field Bp and electron temperature Te, to the functions of ρ.

A trial function for the current density profile is introduced inside 0 < ρ < ρb, Jtf(ρ) = Jex(ρb)× (1 – k2ρ2) / (1 – k2ρb2), where Jex is the experimental function given in [7] and the k2 value isdetermined so that the total current is constant. By changing the magnetic-axis displacementξmax as a parameter, we reconstruct J(ρ), Bp(ρ), and <Bp>. When we put ξmax ≈ 0.28, we
obtain a reconstructed <Bp> profile very similar to the measured Bp profile. Figure 1 shows
the reconstructed J(ρ) of the standard convex profile (red solid line) in contrast to the hollow
profile of Jex(ρ) (blue circle chain). The central value J(0) is 1.35 times larger than Jex(0).
The neo-classical resistivity ηNC(ρ) is calculated from the reconstructed Te profile and theenhancement effect by εeff assuming the uniform Zeff. We estimate a current density profile
Jcal(ρ) ∝ 1/ηNC(ρ) immediately (purple broken line in Fig. 1), and compare it with the reconstructed J(ρ). They agree well with eachother. The enhancement of neo-classicalresistivity is important. If this enhancement
is not taken into account, the estimated J
profile is very peaked near the magnetic axis
as shown by a purple dotted line.
The q profile is of the normal shear with
q(0) ≈ 0.8 in contrast to the reported qex
profile of the weak reversed shear with
qex(0) ≈ 1.1. This result is consistent withthe assumption that the central core plasma
Figure 1. Radial profiles of current density J(ρ)
column is deformed by the m/n = 1/1
for JT-60U hybrid scenario. Experimentallyinternal kink mode. Present model can bereported Jex (blue circle chain), reconstructed J
valid to predict the J profile for the hybrid(red solid line), calculated Jcal (purple dashedscenario, and has to be applied to theline) and calculated one without enhanced neointegrated transport simulation.
classical resistivity (purple dotted line).
[1] C. Gormezano et al., Nucl. Fusion 47, S285 (2007).
[2] A.C.C. Sips et al., Plasma Phys. Control. Fusion 47, A19 (2005).
[3] X. Litaudon et al., 2013 Nucl. Fusion 53, 073024 (2013).
[4] A.H. Boozer, J. Plasma Phys. 35, 133 (1986).
[5] H.L. Berk et al., “Hyper-resistivity Theory in a Cylindrical Plasma”, LLNL Report UCRL-ID142741 (2001).
[6] T.A. Casper et al., Nucl. Fusion 47, 825 (2007).
[7] N. Oyama et al., Nucl. Fusion 49, 065026 (2009).
[8] T. Suzuki, private communication; T. Suzuki et al., “q-profile in Long-pulse Hybrid Discharge in
JT-60U”, 6th ITPA IOS TG Meeting (2011).
[9] T. Takizuka, proc. 30th annual meeting JSPF, Tokyo (2013) 5aE39P:
www.jspf.or.jp/jspf_annual2013/jspf30/pdf/05aE39P.pdf.
[10] S. Tokuda et al., Nucl. Fusion 22, 661 (1982).

Development of implicit MHD code

M. Sato, N. Nakajima, K.Y. Watanabe, Y. Todo and Y. Suzuki
National Institute for Fusion Science

As the beta value increases in LHD (Large Helical Device) plasmas, the Shafranov
shift becomes large and the magnetic well is formed in the core region and the Mercier
unstable region shifts form the core region to the peripheral region. Thus, the pressure
driven MHD instabilities are theoretically unstable in the peripheral region for high beta
LHD plasmas. On the other hand, high beta plasmas with about 5% of the volume
averaged beta value are routinely obtained in the LHD experiments. This implies that
nonlinear saturation level of the pressure driven modes is low enough to maintain suchhigh beta plasmas. In this study, the linear properties and nonlinear phenomena of thepressure driven modes for high beta LHD plasmas have been investigated by MHD
simulation in order to understand why high beta plasmas are stably obtained in LHD.
The MHD simulation has been carried out for three MHD equilibria with β0 =7.4%,
9.4% and 11% by using the MIPS code (MHD Infrastructure for Plasma Simulation)[1]which solves the full MHD equations in the cylindrical coordinates. Here β0 is the central
beta value. The MHD equilibira are constructed by HINT2 code[2]. For these MHD
equilibria, the ballooning modes extending into the plasma peripheral region with achaotic magnetic field are destabilized. Since the linear growth rate decreases as themagnetic Reynolds number increases, the unstable modes are considered to be resistive
modes. In the nonlinear phase, the core region comes under the influence of theinstabilities and the central pressure decreases (Fig.1). There is a tendency that modes
are suppressed as the beta value and/or magnetic Reynolds number increase, which isconsistent with a result that high beta plasmas enter the second stable region of the ideal
ballooning modes as beta increases [3] and remaining destabilized ballooning modes are
considered to be resistive type.
In the MIPS code, the fourth-order explicit Runge-Kutta method is used for the time
integration. The explicit time integration method has the problem of the CFL condition
where the interval of the time step is limited by the fastest physics phenomenon. For
extending the MIPS code to treat the free-boundary problem using the pseudo-vacuum
plasma model, the initial density profile should be small in the pseudo-vacuum region. In
the pseudo-vacuum region, the Alfven velocity increases so that the problem of the CFL
condition becomes increasingly serious. In order to solve this problem, we are developingnew MHD code using implicit time integration method. In the implicit time integration
method, the interval of the time step is not limited by the CFL condition so that thesimulation can be done using the large interval of the time step. For the first step of thedeveloping of implicit MHD code, linear MHD code in two-dimensional xy-space has beendeveloped in order to evaluate of the basic techniques for the implicit solver. In theimplicit MHD code, the fourth-order finite difference method is used for the spatial
derivatives and the first-order implicit Euler method with the Kryolv method is used forthe time integration. In order to improve the convergence of the Krylov metod, thephysics-based preconditioning[4] is introduced. The physics-based preconditioning
method can improve the convergence of the Kryolov method by transforming thehyperbolic problem into the parabolic problem using the Schur decomposition. In order tocheck the effectiveness of the physics-based preconditioning, the dependence of theiteration number of the convergence of the Krylov method on the peripheral plasma
density has been investigated using the two-dimensional implicit linear MHD code. Whenthe physics-based preconditioning is not used, the iteration number for the convergencesignificantly increases as the peripheral plasma density decreases. On the other hand, for
the implicit code with the physics- based preconditioning, the dependence of the iteration
number on the peripheral plasma density is weak so that the physics-based
preconditioning seems to be effective for the free-boundary problem. The extension of thetwo-dimensional implicit MHD code to three-dimensional is under development.

Fig.1. Time evolution of pressure profile on a poloidal section for β0 =9.4% and S=106.

[1] Y. Todo. et al., Plasma Fusion Res. 5 (2010) S2062.
[2] Y. Suzuki et al., Nucl. Fusion 46 L19.
[3] N.Nakajima et al., Nucl. Fusion 46 (2006) 177.
[4] L. Chacón, Phys. Plasmas 15 (2008) 056103.

乱流計測シミュレータを用いた非拡散輸送過程の数値診断
Numerical Diagnostics of Non-Diffusive Transport Process by Use of
Turbulence Diagnostic Simulator
N. Kasuya, M. Sasaki, S. Inagaki, K. Itoh 1, M. Yagi 2 and S.-I. Itoh
Research Institute for Applied Mechanics, Kyushu University, Kasuga, Fukuoka 816-8580,
Japan
1
National Institute for Fusion Science, Toki, Gifu 509-5292, Japan
2
Japan Atomic Energy Agency, Obuchi, Rokkasho-mura, Aomori 039-3212, Japan
E-mail: kasuya@riam.kyushu-u.ac.jp
A simulation of heat source modulation is carried out to understand plasma turbulent transport,
including non-local responses. In this simulation, the pressure source term is modulated in the
nonlinear saturation state of drift-interchange mode turbulence in helical plasmas.
Characteristic response to the modulation can be extracted by the conditional averaging. It isfound that the heat flux increases abruptly after the gradual change of the flux and pressure
gradient on the onset of the additional heating, which is associated with variations of global
modes propagating from core to edge. This process can cause immediate energy transfer todistant locations in magnetically confined plasmas.
Many experiments have revealed that local flux-gradient relation is insufficient to describe
plasma dynamics in torus, and non-local effects should be taken into account to understandtransient transport phenomena, e.g. L/H transitions. To clarify the transport dynamics,
response to active control with additional modulation are studied. In ECH modulation
experiment, conditional average shows sudden change of fluctuation amplitude in spite ofgradual change of temperature, and hysteresis in the flux-gradient relation [1]. There isanother experimental observation of global structures in torus plasmas [2], which is one of thecandidates to cause immediate transfer to separated positions. Therefore, the role of the global
structure on transport must be clarified.
Here, we carried out a global simulation of plasma turbulence with source modulation in torus,
using Turbulence Diagnostic Simulator (TDS) [3]. The TDS is combination of fluid
turbulence codes, measurement modules and analysis routines, to carry out numericalexperiments of turbulent structural formation. In this study, drift-interchange modes in helicalplasmas are analyzed using a reduced MHD model [3]. The pressure profile is sustained withthe applied source, which has the peak at the plasma center. Simulations are performed with
Fourier modes -32  m  32, -8  n  8, where m and n are the poloidal and toroidal mode
number, respectively. Low m, n modes are excited in the linear growing phase, and saturation is
obtained with energy exchanges between various modes by nonlinear couplings. In the
nonlinear saturated state, modes spreading broadly in the radial direction, such as (m, n) = (1,1) and (2, 1), and localized near their rational surfaces, such as (3, 2) and (5, 5), are both
excited. These fluctuations have the self-correlation time of 50 Alfvén time, though those ofhigher m modes are smaller than 10 Alfvén time. The former and latter corresponds to thecorrelation times of the global modes and micro modes, respectively.
Response to additional modulation is investigated to clarify the nonlinear dynamics. The
additional source is switched on and off cyclically. Characteristic response to the modulation

is extracted by the conditional averaging over 100 cycles of the modulation. The additionalsource is only applied near the plasma center. Figure 1 shows the time evolutions of thepressure, its gradient and the heat flux at the position without the direct modulation to see
how the pressure modulation propagates from the center. The additional source is switched on
at t = 0 and off at t = 500 in these figures. In this region, the pressure has time delay from t =0, as the distance from the additional source region becomes larger. The heat flux increases as
the pressure gradient increases, and vice versa. The heat flux also has the time delay, as shownin Fig. 1 (c), and it is accelerated from t ~ 60 to give the same rump-up rate at different radii.
The transition between two states takes place intup ~200, which is much smaller than thecharacteristic time for pressure evolution ~ 500.
The mechanism of this response can be explainedwith the response of the global mode as (m, n) =(2, 1) fluctuation (Fig. 2). The heat source
modulation affects the global mode in r / a < 0.3,
and the modulation propagates outward, as
indicated by the additional line. This propagation
corresponds to those in the pressure and the heat
flux. Typical time scale for the propagation
comes from the transit-time by the diamagnetic
velocity ~ 100 Alfvén time, which is muchsmaller than that of the collisional diffusion ~1000 Alfvén time. The former and latter is theturbulence response time and the global
confinement time, which is the two different timescales included in this simulation. The
modulation can propagate rapidly in the region
where the global mode spreads, which can be
identified with the large radial correlation [4].
When the modulation reaches to the position at r
/ a ~ 0.6, acceleration of the ramp-up takes place.
This is the region, where the nonlinear couplings
cause the variation of the other modes. The
profile modification arises in this region, andpropagates to the other regions with the shortertime scale characterized by the self-correlation
time of higher m modes, which is smaller than 10
Alfvén time. In this way, rapid propagation of theheat modulation mediated by the global modetakes place. The dynamics of their radial profiles
reveal the role of global mode on transport, andthe mechanism of propagation to separated radial
locations.

Fig. 1: Time evolutions of (a) thepressure, (b) absolute value of its
gradient, (c) heat flux at r / a = 0.6 and(d) the source amplitude. The values
are the deviations from the mean value
at each radius.

Fig. 2: Time evolutions of radial
profiles of the global mode (2, 1).
Deviation from the mean value at eachradius is plotted.

[1] S. Inagaki, et al., Nucl. Fusion 53 (2013) 113006.
[2] S. Inagaki, et al., Phys. Rev. Lett. 107 (2011) 115001.
[3] N. Kasuya, et al., Plasma Fusion Res. 8 (2013) 2403070.
[4] N. Kasuya, et al., Plasma Fusion Res. 6 (2011) 1403002.

Numerical simulations of structural formationin resistive drift wave turbulence
M. Sasakia,b*, N. Kasuyaa,b, K. Fukunagac, S. Inagakia,b, K. Itohb,d, M. Yagib,e, Y. Nagashimaa,b,and S.-I. Itoha,baResearch Institute for Applied Mechanics, Kyushu University,
Kasuga, Fukuoka 816-8580, Japanb
Itoh Research Center for Plasma Turbulence, Kyushu University,
Kasuga, Fukuoka 816-8580, Japanc
Interdisciplinary Graduate School of Engineering Sciences, Kyushu University,
Kasuga, Fukuoka 816-8580, Japand
National Institute for Fusion Science, Toki, Gifu 509-5292, Japane
Japan Atomic Energy Agency, Rokkasho, Kamikita, Aomori 039-3212, Japan* Email: sasaki@riam.kyushu-u.ac.jp

Importance of meso-scale structures such as zonal flows and streamers has been
widely recognized for turbulent transport in improved confinement states of magnetizedplasmas [1]. The elementary process for formation of the turbulent structures and theirselection rule is considered in simple cylindrical geometry [2]. Recently, a macro-scale
structure with a finite m, where m is the azimuthal mode number, has been identified inexperiments [3], and the nonlinear coupling with the turbulence has been confirmed [4]. A
streamer [5] and flute type structure [6] are examples of the structures with finite m. Thesetypes of nonlinear structures directly induce radial transport due to their poloidal
inhomogeneities, and have their own characteristic effects to determine the level of transport.
In addition, there can be formed a solitary vortex [7] in a rather less unstable case. Thus, it isimportant to identify their role on transport; whether their formation enhances or suppresses
the total transport [8]. In this view point, it is necessary to understand the fundamental
processes and the competition among these nonlinear structures.
Simulation code Numerical Linear
Device (NLD) is used, which is based ona three-field fluid model (density,
potential and electron velocity parallel tothe magnetic field) for the resistive drift
wave turbulence in cylindrical magnetizedplasmas [2]. A particle source term isintroduced, and self-consistent nonlinear
saturations of turbulence are obtained inthree-dimensional turbulence simulations.
Selective formation of a zonal
flow, streamer and flute structure is
obtained by changing the strength of theparticle source S0, which corresponds tothe driving force of the instability, as isshown in Fig. 1. Here, the definition of
the energies of the flute structure, drift
wave, mode with n = 0, and zonal flow

Fig.1. Source strength dependences of (a) the total
fluctuation energy, and (b) energy partitions among theflute structure, drift wave, mode with n = 0, and zonal
flow.

are in [9], where n is the axial mode number. The streamer appears in the weak turbulent
cases around 2 < S0 < 4. The flute structure, which is linearly stable, becomes dominant justbelow the marginal stability region of the zonal flow around 4 < S0 < 6. The turbulent
structures which enhance the transport, such as the streamer and flute structure, are dominant
in the smaller source cases. The energy partition of the zonal flow exceeds that of the flute
structure around S0 = 6, and the zonal flow is dominant in the larger source cases even withthe existence of the flute structure.
The particle confinement properties in the turbulent states are investigated. The
streamer and the flute structure enhance the transport, and affect the amplitudes of the drift
waves and the zonal flow by nonlinear couplings. The dependence of the total particle number
is shown in Fig. 2. The property of the confinement changes around S0 = 2, where the drift
wave is linearly destabilized, and around S0 = 6, where the zonal flow begins to be effectively
driven. In the cases with S0 < 2, the drift wave is linearly stable, and the diffusion term
determines the density profile, where the total particle number is proportional to the strengthof the source. In the cases with 2 < S0 < 6, where the streamer or the flute structure isdominant, the total particle number
follows the scaling as 0.72 power of thestrength of the source strength. Thisscaling is similar to the L-mode plasma,
where the energy confinement time
follows 0.5 power of the strength of the
input power [10]. In the cases with S0 > 6,
the total particle number increases inproportion to the strength of the source, so
the particle confinement is improved. In
this way, the spontaneous change in thetransport property is obtained even in theFig.2. Source strength dependence of the total particle
cylindrical magnetized plasmas.
number.

The condition for formation of the solitary vortex is also studied. This structure isformed by nonlinear couplings of harmonics with a single dominant unstable mode [7]. Linear
analyses show that a strong diffusion term and /or peaked density profile give the condition.
Based on the results, nonlinear simulations are carried out to form the structure selectively
with a steeper density gradient, which has convergence of waves at the wave front. The
strength of the convergence corresponds to the strength of the localized radial transport, so is
important to evaluate dynamics of the transport.

References[1] P. H. Diamond, et al., Plasma Phys. Control. Fusion 47, R35 (2005).
[2] N. Kasuya, et al., Phys. Plasmas 15, 052302 (2008).
[3] S. Inagaki, et al., Phys. Rev. Lett., 107, 115001 (2011).
[4] T. Kobayashi, et al., Plasma Phys. Control. Fusion 54, 115004 (2012).
[5] T. Yamada, et al., Nature Phys. 4, 721 (2008).
[6] D. Biskamp, et al., Phys. Rev. Lett., 30, 706 (1995).
[7] H. Arakawa, et al., Plasma Phys. Control. Fusion 51, 085001 (2009).
[8] K. Itoh, et al., Plasma Phys. Control. Fusion, 54, 115004 (2012).
[9] M. Sasaki, et al., submitted to Nucl. Fusion (2013).
[10] U. Stroth, et al., Nucl. Fusion, 35, 131 (1995).

Modeling of turbulence in the edge-core coupling region in
tokamak plasmas
Y. Kosuga[1,2][1] Institute for Advanced Study, Kyushu University, Japan[2] Research Institute for Applied Mechanics, Kyushu University, Japan

Predictive modelling of turbulent transport is essential to the success of ITER and
DEMO. While numerical simulation has been succeeding in reproducing experimental
observations, recent validation study identified the problem of short fall[1]. In that study,
turbulent ion and electron heat flux was observed in L-mode plasmas and compared tothe prediction from numerical simulations (TGLF, GYRO, GEM). The basic trend wasthat while both experiments and simulations agree in the core region, the disagreement
becomes prominent toward the edge region. In particular, in the region that connects thecore and the edge with ρ = 0.6 ∼ 0.8, the predicted level is lower by roughly factor of 2
than the observed value. This is the problem of short fall. The short fall issue indicatesthat in order to model turbulence in the edge-core coupling region, new effects are requiredin addition to the conventional turbulence models assumed in the simulations.
An effect that is missed from the conventional modelling is the effect of non-locallity.
There are several candidates that can contribute to the turbulence level in the couplingregion non-locally, as shown in Fig.1.[2] Among these, here we focus on the incoming
holes from the edge[2,3]. These holes are generated in the relaxation process associatedwith the steep gradient in the edge region. In this process, a pair of a hole and a blob isgenerated. Once generated, a blob propagates down the gradient while a hole propagatesagainst the gradient and releases a free energy. Since propagating holes can enter theedge-core coupling region and release free energy stored in the region, holes are importantelement in understanding the phenomenology in the edge-core coupling region.
In modelling the dynamics of holes or blobs, an important variable is the potential
vorticity(PV). For example, a basic model for describing drift wave turbulence, such as
Hasegawa-Mima or Hasegawa-Wakatani model, conserves potential vorticity dq/dt = 0,
where q = (1 − ∇2⊥ )φ for HM model or q = n − ∇2⊥ φ for HW model, φ is electrostatic
potential and n is density. PV conservation constraints the dynamics of porpagaing holes
as shown in Fig.2. Due to the PV constraint, hole growth must involve accelerationof flows. Thus flow dynamics is important to understand self-consistent dynamics ofpropagating holes. The overall coupled dynamics determines the level of hole-induced
transport, which can be a part of relaxation process in the edge-core coupling region.
Fruitful discussion with P.H. Diamond, O.D. Gurcan, G. Dif-Pradalier, S.-I. Itoh, K.
Itoh, S. Inagaki, M. Lesur is acknowledgeged. This work was supported by Grants-in-Aid
for Scientific Research of JSPF of Japan (21224014,25887041).
References: [1] T.L. Rhodes, et al.: Nucl. Fusion 51 (2011) 063022 [2] Y. Kosuga, et al.:Nucl. Fusion 53 (2013) 043008 [3] Y. Kosuga and P.H. Diamond: Plasma Fusion Res. 8(2013) 2403080

Density profilen(r)‘No Man’s Land’✛✲✲

Tokamak edge✲✲❑❆ Incoming holes❆

DW propagationfrom the core+ ‘beach effect’❆

Turbulence✛

spreading

✛

❆
❆

❆

from the edge

❆

❆❆
❆❆❯

✲r

Figure 1: Several processes that can contribute to turbulence in the edge-core couplingregion. These include, but not limited to, incoming holes from the edge, turbulence
spreading from the edge, drift waves propagating from the core with beach effect. These
processes can non-locally enhance turbulence level in the connecting region.

dn/dt = 0 (?)
and hn(x)idq/dt = 0 (!)
and hq(x)i

❍❍
✛❍

❍❍
✛❍

❍❍

❍❍

❍
❍

ñ grows→ hṽx ñi → transport❍❍q̃ grows
hṽx ñi → transport→ hṽx q̃i ✒❘ hṽx ∇2 φ̃i → flow❅⊥

Figure 2: A cartoon for the contrast between dn/dt = 0 and dq/dt = 0. If dn/dt = 0,
propagating holes grow by extracting free energy and cause transport. On the other
hand, if dq/dt = 0, not only holes drive transport, but also they drive flows. The coupling
between propagating holes and flows determine overall dynamics and transport level.

Anomalous configuration of hydrogen trapped in tungsten
vacancy and its application to plasma facing material
Kazuhito Ohsawa1 and Yuuji Hatano2

1Research Institute for Applied Mechanics, Kyushu University
2Hydrogen Isotope Research Center, University of Toyama

Introduction
Investigation of interaction between hydrogen (H) isotopes (D and T) and tungsten(W) is an important subject in the field of plasma facing materials (PFMs) because W
and its alloys are promising as divertor armor tiles installed in future fusion reactors.
We investigate stable configurations of H atoms trapped in bcc transition metals interms of first-principle calculations. Anomalous H configurations are found associatedwith multiple H atoms trapped in a W vacancy. We estimate correct binding energy ofeach H isotope to the W vacancy on the basis of the valid H configurations. It is enabled
us to study on isotope exchange mechanism and vacancy-hydrogen (V-H) cluster
concentrations at finite temperature.
Methodology
First-principle calculations are performed using Vienna ab initio simulation package(VASP) based on density functional theory. Binding energy, ek, of H or H isotopes to a Wvacancy is defined, as follows.

ek=E[Wn-1V]-E[Wn-1VHk]-k(E[WnHT]-E[Wn]),
where the function E is the cohesive energy of the supercell. For example, Wn-1VHk
indicates a supercell composed of n-1 W, a vacancy V, and k H atoms. V-H cluster
concentrations at finite temperature are calculated in an equilibrium thermodynamic
model.
Stable configuration of H atoms
Tetrahedral and octahedral sites (T site and O site) are typical interstitial positions in
bcc lattice, as shown in Fig. 1 (a). It has been reported that an H atom prefers to belocated at T site in bcc transition metals. While, a single H atom trapped in a vacancy islocated on the inner surface of the vacancy close to an O site, as shown in Fig. 1 (b).

Fig. 1: (a) Typical interstitial sites, T site
and O site, in BCC metals. (b) Stable
position of single H trapped in vacancy
close to O site.

Fig.2: Stable configurations of multiple H
atoms in W vacancy. Characteristic
configurations of 4, 6, 12H atoms in W
vacancy.

A maximum of 6 H atoms is believed to be trapped in a vacancy in bcc transition
metals because there are 6 O sites next to the vacancy. However, H atoms in a Wvacancy are shifted in the <100> directions from the O sites, as the number of H atoms
in the vacancy increases. Figure 2 shows some characteristic stable configurations ofmultiple H atoms in the W vacancy. (a) Stable configuration of 4 H atoms is tetrahedral
structure but planar one. (b) Stable configuration of 6 H atoms is complicated. 4 H
atoms out of the 6 are displaced from O sites by about 18 % of the lattice constant in the<100> directions, and the remainder of the H atoms are slightly shifted from O sites.
(c) 12 H atoms in a W vacancy are located close to T sites.
It has been reported that a vacancy in usual bcc transition metals is capable ofaccommodating 6 H atoms. However, a maximum of 12 H atoms can be accommodated
in a W vacancy.
Application to study on plasma facing material
Tritium absorption and desorption are an important subject for PFMs in fusion
reactors because T is a radioisotope whose half-life is 12 years. Isotope exchange is
expected to be an effective method to induce T desorption from the PFMs. The binding
energy of T to a W vacancy is supposed to be smaller than that of other H isotopes. The
isotope exchange is caused by the energy difference. It is necessary to calculate correct
zero-point vibrational energies on the basis of the valid H configurations. We estimate
binding energies of H and T to a W vacancy and confirm the expectation, as shown in
Fig. 3.
V-H cluster concentrations are calculated in an equilibrium thermodynamic model atfinite temperature and H concentration in a bulk W. A maximum of 12 H atoms can beaccommodated in a W vacancy at 0K, according to first-principle simulations. However,
the major V-H clusters are composed of 6 H atoms in a wide range of temperature and H
concentration. Figure 4 is an example of V-H cluster distribution.

Fig3: Binding energy of H and T to W
vacancy as a function of the number of H
atoms trapped in W vacancy.

Fig. 4: Number ratio of V-H cluster
distribution at H concentration 10-3 and
T=700K.

Particle Simulation of Peripheral Plasmasin Tokamaks using PARASOL
1 Shuhei AZUMA, 1 Atsushi FUKUYAMA, and 2 Tomonori TAKIZUKA
1 Graduate School of Engineering, Kyoto University
2 Graduate School of Engineering, Osaka University

1

Introduction

make two divertor conﬁgurations possible; anupper null (UN) and lower null (LN) conﬁguration.
We consider that plasma consists of electrons and single-species ions, and any ion isassumed to be charged singly for simplicity.
The motion of electrons is approximated bytheir guiding-center motion, while that of ionsis fully followed. The electric ﬁeld is calculated self-consistently with the PIC method, on
the assumption that boundary walls surrounding plasma are electrically conductive. Anomalous transport is simulated with a Monte-Carlo
random-walk model, and the eﬀect of Coulomb
collisions is added by a binary collision model.
A particle is assume to be absorbed whenreaching divertor plates, and to be reﬂected
when arriving at the other side walls. Given
that ∆Nloss is the number of lost ions in one time
step, in the following time step, (1 − R)∆Nloss
pairs of an electron and an ion are supplied in
the central region H, and the R∆Nloss pairs in
the recycling area C1 and C2, at a prescribedtemperature (see Fig. 1 for H, C1, and C2). HereR represents a recycling rate and have the values of 0, 0.5, 0.7, 0.8, and 0.9 in this study. As
for other paramters, the number of simulation
ions Ni is 106 and the number of spatial cells
MR × MZ is 320 × 512. The pitch of the magnetic ﬁeld Θ is 0.2, and the aspect ratio is set
as A = 5.5. The mass ratio mi /me is 1800, The
normalized ion Larmor radius ρi /a ∼ 0.02, thecollisionality Lk /lmfp ∼ 1, and the normalized
diﬀusion coeﬃcient Danom /aCs ∼ 10−5 (Cs is asound speed).

In tokamak reactors, peripheral plasma is expected to play an important role as an acceptor of high heat ﬂux and particles, such as fuel,helium ash, and impurities, from core plasma.
The ﬂow in peripheral plasmsa is also expected
to take them to the divertor, where they are removed [1]. Their permformance have been so
far examined by ﬂuid simulations, while thesesimulations couldn’t give a satisfactory exlanation for some results of experiments conducted
to investigate the ﬂow [2].
On the other hand, recent studies have demostrated, based on particle simulations using
PARASOL, that the eﬀects of ﬁnite-orbit-size ofions are essential for the formation of the plasma
ﬂow [3], and have developed a new model of theﬂow, called “ion-orbit-induced ﬂow”, for comprehensive divertor simulation with ﬂuid models[4]. Since the simulations for the previous studies did not take into account recycling or radiation cooling, however, little is still known about
their eﬀects on the peripheral plasma ﬂow.
We examined in this study, by using full particle simulation code PARASOL, how recyclingin the vicinity of the divertor aﬀects peripheral
plasma ﬂow. Also note that we do not consider
here radiation cooling in order to focus on theeﬀect of recycling.

2

PARASOL simulation models

The two demensional PARASOL (PARASOL2D) is a full particle simulation code which
simulates a scrape-oﬀ layer (SOL) and divertor plasma in the toroidal system with cylindrical coordinates, employing the particle-in-cell(PIC) method, along with a binary collisionmodel [3]. The spatial conﬁguraion is shownin Fig. 1. It should be noted that adjusting theratio of upper- and lower-divertor coil currents

3

Results

Figure 2 shows, for the values of recycling rates,
simulation results of the patterns of the SOL
ﬂow parallel to the magnetic ﬁeld in the cases1

of (a) the LN and (b) the UN conﬁguration,
where the direction of ion ∇B drift is downward. In the ﬁgure, red expresses the co-ﬂow
to the plasma current, while blue refers to thecounter-ﬂow. The red co-ﬂow and blue counterﬂow represent anti-clockwise and clockwise ﬂow
in the poloidal cross-section. It has proved evident from these results that recycling has little
impact on SOL plasma ﬂow patterns.

tionally called calculated Mth here.) We notehere that FΓ = 1/(1 − R), and the index p denotes a quantity measured in front of diverter
plates. Figure. 3 shows the comparison of measured Mth with calculated Mth by FΓ . It follows
from these data that Mth is less subject to therecycling eﬀect in the 2D toroidal system where
the ion-orbit-induced ﬂow dominates.

References[1] A. Loarte et al, Progress in the ITER
Physics Basis Chapter4:Power and particle
control Nucl. Fusion 47 S203 (2007)[2] N.Asakura and ITPA SOL and Divertor
Topical Group, J. Nucl. Mater. 363-365
41 (2007)[3] T. Takizuka, K. Shimizu, N. Hayashi, M.
Hosokawa, and M. Yagi, Nucl. Fusion 49
075038 (2009)[4] T. Takizuka, K. Hoshino K. Shimizu, andM. Yagi, Contrib. Plasma Phys. 50 No.3-5,
267―272 (2010)[5] T. Takizuka, M. Hosokawa, K. Shimizu, J.
Nucl. Mater. 753 (2001) 290-293

We also examined the ion Mach number parallel to the ﬁeld at the divertor throat near the
X point, Mth , to gain a quantitive understanding of the recycling eﬀect on√the ﬂow. Considering Mth = |Vi |/Cs (Cs ≡ T /mi ), and using
Vi and T from simulation results, we can obtain Mth (called measured Mth for convinience
here ). On the other hand in the study of the
one dimensinal SOL plasma system, it have been
veriﬁed that Mth can√be approximated by using
the relation, Mth ∼ Tp /Tth /2FΓ , in the caseof high recycling rates [5]. (Also it is convenzbw

particle source by recycling
divetor plate

R=0

θ

C1

∂/∂θ = 0

R =0.9400

300

particle source-aw

R =0.5

500

C2

aw R

H

(a) LN
200

100

Upper-null case-bw

R = R0φ=0
0Figure 1: Space conﬁguration for theupper null case

0

100

200

300

0

100

200

300

0

100

200

300

0

100

200

300

0

100

200

300

0

100

200

300

500

400

measured
calculated

300(b) UN

Mth

200

100

FΓ

Figure 3: FΓ vs. measured and calculated Mth0

Figure 2: Simulation results of plasma ﬂow patterns

2Progress of ITER Integrated Modeling Activity
A. Fukuyama
Department of Nuclear Engineering, Kyoto University, Kyoto, JapanFor the success of the ITER project, the role of integrated modeling of burning plasmas
is essential. It is necessary for operation support (preparation, execution, reconstruction)and plasma research (model validation, scenario development, upgrade of equipments).
The modeling activity is a part of activities in Science Division of Plasma Operation Department of ITER organization. In order to promote collaboration with ITER participating
parties, the Integrated Modelling Expert Group (IMEG) takes the role of adviser for thedevelopment of integrated modelling codes and topical groups of International Tokamak
Physics Activity (ITPA) contribute construction of experimental database, development
and validation of physics models, and operation scenario development.
The present main activity of ITER integrated modeling is to complete the conceptual
design of Integrated Modelling Analysis Suite (IMAS). The conceptual design has been
carried out by a contract with a research group from EU. IMAS will be constructed basedon a common infrastructure and various components implementing physics models are
combined with each other for specific purpose. The components will be developed mainly
by participating parties conforming with the ITER Modelling standard.
The infrastructure includes• Physics data model (Data structure, Data specification)• Component model (Application program interface)• Workflow tool (Compilation of components, data flow)• Monitoring (Real time data glancing)• Life cycle (How to upgrade data models)For physics data mode, rules and guidelines have been documented and a preliminary
set of Interface Data Structure (IDS) are available. Table 1 indicates the list of IDS in
Physics Data Model v2.0. Each IDS includes various data. As an example, the data
definition of IDS core profile is shown in Table 2. The structure data ids properties isTable 1 List of IDS in v2.0
IDS name
actuator
controllers
core profiles
core sources
core transport
em coupling
equilibrium
magnetics
pf active
pf passiveschedulesdn
simulation
temporarytf

Description
Generic simple description of a heating/current drive actuator
Feedback and feedforward controllers
Core plasma radial profiles
Core plasma source terms (for the transport equations)Core plasma transport
Description of the axisymmetric mutual electromagnetics
Description of a 2D, axi-symmetric, tokamak equilibrium
Magnetic diagnostics for equilibrium identification and plasma shape control
Description of the axisymmetric active poloidal field (PF) coils and supplies
Description of the axisymmetric passive conductors, currents flowing in them
Description of Pulse Schedule
Description of the Synchronous Data Network parameters and the signals on it
Simulation parameters for a general simulation
Storage of undeclared data model components
Toroidal field coils (all assumed to be identical)

1

a common data for all IDS. Since the data includes the time history, most of data is
FLT 2D, two-dimensional floating data with radial and time index. The component model
is similar to the interface of BPSD. As a workflow tool, KEPLER workflow engine is used
and the transport code CORSICA was implemented on it. The procedure to update theinfrastructure, especially the date model, is under discussion.
As an BPSI activity, the interface between BPSI (data model and application interface)
and the ITER physics data model will be developed. The collaboration of the integrated
codes in Japan, TOPICS, TASK, TASK3D, and TOTAL, will be promoted, and futurecontributions to IMAS should be discussed.
Table 2 Definition of IDS core profiles
Full path nameids properties
rho tor norm
vacuum toroidal fieldpsitet i averagenen e fast
ions
ions/aions/z ion
ions/z nions/label
ions/n iions/n i fastions/t iions/p iions/v tor iions/v pol iions/charge statesn i total over n emomentum tor
zeﬀ
pep i totalpressure thermalpressure perpendicularpressure parallel
j totalj tor
j ohmic
j non inductive
j bootstrap
conductivity parallele field parallelqmagnetic shear
global quantities
code

Description
Interface Data Structure properties
Normalised toroidal flux coordinate
Characteristics of the vacuum toroidal field
Poloidal magnetic flux
Electron temperature
Ion temperature
Electron density (thermal+non-thermal)
Density of fast (non-thermal) electrons
Quantities related to the diﬀerent ion species
Mass of atom
Ion charge
Nuclear charge
String identifying ion
Ion density (thermal+non-thermal)
Density of fast (non-thermal) ions
Ion temperature
Ion pressure
Toroidal velocity
Poloidal velocity
Quantities related to the diﬀerent charge states
Ratio of total ion density over electron density
Total plasma toroidal momentum
Eﬀective charge
Electron pressure
Total ion pressure
Thermal pressure (electrons+ions)
Total perpendicular pressure
Total parallel pressure
Total parallel current density
Total toroidal current density
Ohmic parallel current density
Non-inductive (includes bootstrap) parallel current density
Bootstrap parallel current density
Parallel conductivity
Parallel electric field
Safety factor
Magnetic shear
Various global quantities derived from the Core Profiles
Generic description of the code specific parameters

2

Data typestructure
FLT 1D
structure
FLT 2DFLT 2DFLT 2DFLT 2DFLT 2D
structure
FLT 0DFLT 0DFLT 0D
STR 0DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2D
structure
FLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2DFLT 2D
structure
structureProgress of integrated transport analysis suite, TASK3D-a (analysis
version), and its extensive application to LHD experiment
M.Yokoyama, R.Seki, C.Suzuki, S.Satake, M. Osakabe, K.Ida, S. Murakami*, H.Yamaguchi*,A.Fukuyama*, and LHD Experiment Group
National Institute for Fusion Science, 322-6, Oroshi-cho, Toki-city, 509-5292, Japan*Department of Nuclear Engineering, Kyoto University, Kyoto 615-8540, JapanThe Large Helical Device (LHD) has steadily extended its parameter regime [1]. More
importantly, not only the parameter extension, but also the increased diagnostic capability has
provided a wide-range database to be exploited for the accurate physics discussion. Integrated
transport analysis suite, TASK3D-a (Analysis version) [2], has been developed for conducting
automated energy confinement analyses for NBI-heated plasmas in LHD. The capability of
TASK3D-a has been gradually increasing. In this abstract, significance and usefulness of
TASK3D-a are described from the viewpoints of its close link to LHD experiment data andthen making systematic analyses possible.
Calculation flow of TASK3D-a consists of the three-dimensional equilibrium specification,
NBI deposition and energy balance calculations, so that they are sequentially executed in an
automated manner. It has an LHD experiment data interface part, which has a direct link to aso called LHD Kaiseki Data Server [3] and TSMAP (real-time coordinate mapping system)[4]. Processed diagnostics data are registered onto the Kaiseki Data Server in a common
format (ASCII file, so called eg-format), and then mapping of profile data such as temperatureand density (which are required for transport analyses) from the real geometry to the effective
minor radius, reff, defined by TSMAP. Then mapped profile data are transferred to TASK3D-a
as inputs. Some of TASK3D-a results (NBI deposition profile, ion and electron heatdiffusivity etc.) are then registered to the Kaiseki Data Server for common use. Users justneed to login to the Kaiseki Data Server
for utilizing these results. The dataflow
described above is conceptually drawn in
Fig.1. The words written with small
characters are names of registered data (socalled eg data). The logo near the“TASK3D-a” is newly introduced logo for
TASK3D collaborative research.
In such a way, TASK3D-a, fully
integrated into LHD experiment database,has smoothly utilized experimental data forincreasing cases of analyses. So far, highion temperature (Ti) plasmas [5] (Ti~5-7
keV, Ti>Te at the core with the density of
1~2x1019 m-3) and medium-to-high density
plasmas [6] (Ti~Te~3-4 keV at the corewith the density of 3~5x1019 m-3) havebeen extensively analysed. High-Te
plasmas will also be analysed once ECH
modules are fully integrated into
TASK3D-a. So far, the order of 104

Fig.1

The close link between TASK3D-a(logo
isattached)andLHD
experiment data (centralized to LHD
Kaiseki Data Server) is schematicallyshown.

timings has been analysed from multiple shots, including cases without Ti profile
measurement (only Te and ne are available).
Figure 2, as an example of a systematic plot, shows the ion and electron heat diffusivity as afunction of temperature ratio, Te/Ti, for reff/a99~0.4 based on dynamic transport analyses. Herea99 is the plasma minor radius in which 99 % of the total stored energy is confined. Generally,
data for Te/Ti <1 correspond to high-Ti plasmas, and those for Te/Ti>1 to medium-to-highdensity plasmas. The diffusivity is normalized by Gyro-Bohm temperature scaling, T1.5. The
tendency is recognized that the normalized ion (electron) heat diffusivity increases (decreases)
as Te/Ti is increased. It can be considered that high-Ti plasmas in LHD are situated in a Te/Ti
regime with smaller ion diffusivity. From this systematic plot, it is implied that plasmas with
Te~Ti are favourable from the viewpoint ofsimultaneously small ion and electron heatdiffusivity. This implication is consistent withrecent trials for increasing Te in high-Ti
plasmas in LHD (from Ti>Te towards Ti~Te)
ion
through the increase of available ECH power.
electron

1.5

1

i,dy

/T

i

,χ

e,dy

/T

e

10

χ

The important feature of a huge TASK3D-a
analysis database is that it stores plasma
parameters, experimental power balance,
neoclassical diffusion and dimensionlessparameters etc., all in “PROFILES”. Thus,
probably based on the statistical approach, it
is anticipated to deduce heat diffusivity
scaling with radial dependence, also
distinguishing ion and electron channel. This
will be one of research targets in the near
future and will be reported in a separate
opportunity. It should be far beyond theglobal energy confinement scaling [7,8], andthus significant increase of the predictive
capability is expected. Such deduced scaling
will be implemented into TASK3D-p, thepredictive version of TASK3D.

1.5

2[m /s]

dytrns6-gsrake_Comparison_reffa99=0.4_HighTi_Med-High-Density

0.1
0.4

0.6 0.81
1.2 1.4 1.6
Te/Ti (@reff/a99~0.4)
Fig.2 Example figures describing
systematic dependence based on
TASK3D-a analysis database (multiple
shot-timings). Ion and electron heatdiffusivity evaluated with the dynamictransport as a function of temperature
ratio, Te/Ti at a specific radius(reff/a99~0.4).

Acknowledgements
The authors are grateful to LHD experiment group for making the experimental data available.
This work has been supported by the NIFS Collaborative Research Programs,
NIFS11KNTT008 and NIFS11UNTT006. One of authors (M.Y.) appreciates a grant-in-aid fromthe Future Energy Research Association.
References[1] O.Kaneko, H.Yamada et al., Nuclear Fusion 53 (2013) 104015.
[2] M.Yokoyama et al., Plasma and Fusion Res. 8 (2013) 2403016.
[3] M. Emoto et al., Fusion Eng. Des. 81 (2006) 2019.
[4] C.Suzuki et al., Plasma Phys. Control. Fusion 55 (2013) 014016.
[5] K.Nagaoka et al., Nucl. Fusion 51 (2011) 083022.
[6] A.Dinklage, M.Yokoyama et al., Nucl. Fusion 53 (2013) 063022.
[7] For example, ITER Physics Basis, Nucl. Fusion 39 (1999) 2175.
[8] H.Yamada et al., Nucl. Fusion 45 (2005) 1684.

Particle transport simulation in tokamak plasma by
TASK/TX
KANAMORI Ryosuke, FUKUYAMA Atsushi
Graduate School of Engineering, Kyoto Univ.
In a magnetic confinement nuclear fusion reactor, fusion output power is strongly influenced
by fuel ion density. The density profile of fuel ions is govered by ionization and transport
though the particle transport in tokamak plasmas is not yet well-known. In previous particle
transport analyze, diﬀusion equation of particle flux has been solved for only ion particle
transport and electron density is obtained by quasi-neutrality condition. In this report, we
simulate particle transport by the TASK/TX code and try to resolve various issues in particle
transport. TASK/TX is a 1D dynamic transport simulation code based on two-fluid equations.
It solves particle transport of electron and ion simultaneously, and solves Gauss’s law without
using quasi-neutrality condition. It is able to describe the time evolution of density, rotationand temperature for each species. Our goal is to construct particle transport modeling which
agrees with experimental observations by improving physics models and various parameter.
The set of 1D two-fluid equations is composed of equation of continuity for density ns ,equations of motion for radial velocity usr , poloidal and toroidal rotation usθ , usφ , and heattransport equation for internal energy 32 ns Ts and written as
(
)
1 ∂
1 ∂
∂ns
∂ns
=−
(rns usr ) +rDm vT s
+ Ss
∂t
r ∂r
r ∂r
∂r
∂
1 ∂
1
∂
(ms ns usr ) = −(rms ns u2sr ) + rms ns u2sθ −(ns Ts ) + es ns (Er + usθ Bφ − usφ Bθ )
∂t
r ∂r
r
[ ∂r
]
∂
1 ∂
1 ∂
∂ ( usθ )(ms ns usθ ) = − 2 (r2 ms ns usr usθ ) + 2
r3 ms ns µs+ es ns (Eθ − usr Bφ )
∂t
r ∂r
r ∂r
∂r r
[
]
1 ∂
∂
NC
C
W
L
N
CX
+
rDm vT s (ms ns usθ ) + Fsθ
+ Fsθ
+ Fsθ
+ Fsθ
+ Fsθ
+ Fsθ
r ∂r
∂r
(
)
∂usφ
∂
1 ∂
1 ∂
(ms ns usφ ) = −(rms ns usr usφ ) +rms ns µs+ es ns (Eθ + usr Bθ )
∂t
r ∂r
r ∂r
∂r
[
]
∂
1 ∂
C
W
L
N
CX
+
rDm vT s (ms ns usφ ) + Fsφ
+ Fsφ
+ Fsφ
+ Fsφ
+ Fsφ
r ∂r
∂r
(
)
(
)
∂ 3
1 ∂ 5
3
∂Ts
ns Ts = −rusr ns Ts − rns χs+ es ns (Eθ usθ + Eφ usφ )∂t 2r ∂r 2
2
∂r
[
]
1 ∂
∂
+
rDm vT s (ns Ts ) + PsC + PsL + PsR + PsRF
r ∂r∂r
where terms including Dm represent diﬀusive eﬀects due to magnetic surface destruction. The
perpendicular viscosity µs and thermal conductivity χs represent anomalous transport due tothe turbulence. The particle source Ss , neoclassical viscous force FsNC , classical collisional momoentum transfer FsC , forces due to the interaction with turbulence electric field FsW , paralleltransport loss in the SOL region FsL , friction force with neutrals FsN , charge-exchange force
FsCX , collisional energy transfer power PsC , collisional energy loss power PsL bremsstrahlung
power PsR , and direct RF heating power PsRF are calculated from local quntities.
In the SOL region, a < r < b, where a is the plasma radius and b is the wall radius, losses
along with magnetic field line is dominant. The loss rate of particle and ion heat are given by1νL = Cs /L where Cs is the sound speed and L is the length of magnetic field line. Electron
heat loss is modeled by collisional heat conductivity
5/2

νLTe =

κ0 Te
L 2 newhere κ0 =

(4π0 )2
1/2

me Zeﬀ e4 ln ΛeWhen the electron heat loss flux rate exceeds the electron thermal velocity, the heat limiting
model νLlim = vTe /L is introduced to giveν̂LTe =

νLTe · νLlimνLTe + νLlim

Magnetic surface destruction by magnetic perturbation in plasma edge region causes particle and heat diﬀusion. In the case of low collisionality, the radial diﬀusion coeﬃcient Ds is
proportional to the thermal velocity vT s .
D v(rmin < r < rmax )
(∆r)2
m Ts
=
Ds ≡
0
τ
(otherwise)
where the magnetic diﬀusion rate in length λ is given by
Dm ≡

(∆r)2
λ

First, we analyse the influence of heat limiting model. The dashed curve is the result withheat limiting model, and the solid one is classical model in Fig.1. With the heat limiting
model, abrupt change of electron temperature at the plasma edge disappears.

0.6
0.4
0.2
0.0
0.90

0.95

1.00

1.05

1.0
0.8

2.5
2.0

0.6

1.5

T e[keV]n e [10 19/ m3 ]T e [keV]

0.8

0.4
0.2

r/a

0

1.0
0.5

-4
-8
-12
-16
0.0 0.2 0.4 0.6 0.8 1.0
r/ a

0.0
0.0 0.2 0.4 0.6 0.8 1.0
r/ a

0.0
0.0 0.2 0.4 0.6 0.8 1.0
r/ a

1.10

Radial electric field

Electron temperature

Electron density

E r [kV / m]Electron temperature

Fig.1

Next, we analyze eﬀect the of magnetic surface destruction. The solid curve in Fig.2 isthe case without destruction, and the dashed one is the case with Dm = 7.0 × 10−10 m−1 for
0.85 < r/a < 1.05. The radial electric field changes strongly though there is little change inthe density profile.

1.5

0.2
0.0
0.0 0.2 0.4 0.6 0.8 1.0
r/ a

1.0
0.5
0.0
0.0 0.2 0.4 0.6 0.8 1.0
r/ a

1.2
1.0
0.8
0.6
0.4
0.2
0.0
0.0 0.2 0.4 0.6 0.8 1.0
r/ a

Radial electric field0

E r [kV / m]

0.6

T i [keV]

2.5
2.0

T e[keV]n e [10 19/ m3 ]

1.0
0.8
0.4

Ion temperature

Electron temperature

Elesctron density

-4
-8
-12
-16
0.0 0.2 0.4 0.6 0.8 1.0
r/ a

Fig.2

We have examined some parameter dependences in particle transport simulation by TASK/TX
to improve the agreement with experimental observations. Electron temperature becomessmooth at the plasma edge by using the heat limiting model. More improvement of themodeling is required to describe the eﬀect of magnetic surface destruction.
2

DIII-D におけるアルフベン固有モードと高速イオン輸送のシミュレーション
Simulation of Alfvén eigenmodes and fast ion transport in DIII-D experiment
Y. Todo1), M. A. Van Zeeland2), A. Bierwage3), W. W. Heidbrink4)
1)
National Institute for Fusion Science, 322-6 Oroshi-cho, Toki, Gifu 509-5292, Japan
2)
General Atomics, P.O. Box 85608, San Diego, California 92186-5608, USA
3)
Japan Atomic Energy Agency, 2-166 Omotedate, Obuchi, Rokkasho 039-3212, Japan
4)
University of California, Irvine, California 92697, USA
E-mail contact of main author: todo@nifs.ac.jp
A multi-phase simulation that is a combination of classical simulation and hybrid simulationfor energetic particles interacting with a magnetohydrodynamic (MHD) fluid is developed tosimulate the nonlinear dynamics on the slowing down time scale of the energetic particles. The
hybrid simulation code MEGA is extended with realistic beam deposition profile, collisions,and losses, and is used for both the classical and hybrid phases [1]. The multi-phase simulation
is applied to DIII-D discharge #142111 where the fast ion spatial profile is significantly
flattened due to multiple AE modes [2]. In the simulation, the stored fast ion energy is saturated
due to AE modes at a level lower than in the classical simulation. After the stored fast ion
energy is saturated, the hybrid simulation is run continuously. It is demonstrated that the fast
ion spatial profile is flattened due to the interaction with the multiple AE modes with amplitude
v/vA~δB/B~O(10-4). The fast ion pressure profile is close to that in the experiment. The dominant
AE modes are toroidal Alfvén eigenmodes (TAE modes). The amplitude of the temperaturefluctuations brought about by the TAE modes is of the order of 1% of the equilibrium
temperature. These are consistent with electron cyclotron emission measurement in theexperiment. Temperature fluctuations due to two of the toroidal Alfvén eigenmodes (TAE
modes) are compared in more detail with the experimental measurement. The radial profiles
and phases are in very good agreement with the measurement, and the amplitude is also inagreement within a factor of 2.
We have investigated nonlinear evolution of Alfvén eigenmodes and beam ions in DIII-D
discharge #142111 at t=525ms. The equilibrium is reconstructed with the EFIT code, and arealistic beam ion deposition profile is given by the TRANSP code. Collisions (slowing down,
pitch angle scattering, energy diffusion) with realistic physics parameters are taken into account.
Particle losses are assumed to take place at the plasma boundary (r/a=1). 8 million particles are
injected with a constant time interval in 150ms, although both the classical and multi-phasesimulations are terminated before t=150ms. The beam injection power is 4.95MW with the fullenergy component. The half and third energy components are not included in this simulation.
In the experiment, AE modes with toroidal mode number n=1-5 are observed at t~525ms [2].
Then, we restrict the energetic particle drive in the simulation to n=1-5 in order to reduce
numerical noise. The numbers of grid points are (128, 128, 256) for (R, , z) coordinates,
respectively. Figure 1 shows the time evolution of stored fast ion energy. In the multi-phasesimulation, the classical simulation was run for 4ms and alternately the hybrid simulation was
run for 1ms. This set was repeated until stored fast ion energy is saturated at t=75ms. Beyond
t=75ms, the hybrid simulation with MHD is continued, and the MHD fluctuation reaches to asteady level. Figure 2 compares the fast ion pressure profiles among the classical and multiphase simulations, and the experiment. We see in Fig. 2 that the fast ion pressure profile in themulti-phase simulation is significantly flattened due to AE modes. At t≥75ms, the beam
injection is balanced with transport and energy loss due to collisions and the multiple AE modes.

We analyzed the bulk temperature fluctuation spectra with toroidal mode number n=1-5 at
r/a=0.4 for t≥75ms. We found the temperature fluctuation normalized by the equilibrium
temperature is δT/T~O(10-2) that is comparable to the experiment. We have compared in detailthe temperature fluctuation profile in the simulation with the ECE measurement for toroidal
mode number n=3 and 4. For the comparison with the electron temperature fluctuation in theexperiment, we simply divide the MHD temperature fluctuation by a factor of 2. We found very
good agreement in spatial profile, amplitude, and phase between simulation and experiment.
Figure 3 compares temperature fluctuation profiles for n=3. We have agreement in amplitude
within a factor of 2 for both n=3 and n=4. Figure 4 demonstrates excellent agreement inoscillation phase profile of temperature fluctuation for n=3. The interaction with energetic
particles relaxes the restriction of ideal MHD model where oscillation phase is 0 or  becauseof the Hermitian property. The agreement in fast ion pressure profile and temperaturefluctuation profiles validates the multi-phase simulation based on MHD model.

Fig. 1. Evolution of stored fast ion energy forclassical and multi-phase simulations with MHD
perturbations.

Fig. 3. Comparison of temperature fluctuation
profile between simulation and experiment for
n=3. The left vertical axis is for simulation,
while the right for experiment.

Fig. 2. Comparison of fast ion pressure profile
among multi-phase simulation with MHD
perturbations, classical simulation, andexperiment.

Fig. 4. Comparison of temperature fluctuation
phase profile between simulation andexperiment for n=3.

References[1] Todo Y. et al submitted to Nucl. Fusion[2] Van Zeeland M. A. et al 2012 Nucl. Fusion 52 094023

Scheme-independent error in entropy conservationin turbulent kinetic simulations
M. Lesur

Itoh Research Center for Plasma Turbulence, Kyushu University
Resonant interactions between particles and waves are relevant to a wide range ofcollisionless laboratory and space plasmas. Particle trapping in the troughs of electric
potential is a fundamental kinetic nonlinearity, which leads to the formation of vortex-likestructures in phase-space [1,2], especially in regimes of strong resonance. These structures
have important impacts on wave stability, turbulent transport, intermittency, anomalous
resistivity, heating, flows, etc. [3,4]Accurate numerical simulation of such kinetic nonlinearities is challenging. Indeed,trapping involves the filamentation of phase-space (phase-space mixing). In thecollisionless limit, numerical discretization replaces collisions to dissipate the smallest
scales. A concern is that numerical dissipation and other simple models of collisions are
artificial, and may impact the physics of interest. In this work, I show that the evolution ofobservables is surprisingly robust to changes of numerical treatment, which suggests thatthe details of small-scales dissipation are irrelevant. In particular, a systematic error of 15%
in entropy conservation is found, regardless of the numerical treatment, without impacting
observables.
In this work, as a paradigm for kinetic models in the presence of strong resonances, I
restrict the analysis to one-dimensional ion-electron plasma with an initial velocity drift.
The mass ratio is mi /me = 4, the ion and electron temperatures are equal, and the initial
velocity drift is vd = 4.2vT,i (slightly above linear instability threshold), where vT,i is the ion
thermal velocity. No collision operator is included. The simulation code is COBBLES [5],which is based on the Constrained-Interpolation-Profile, Conservative Semi-Lagrangian(CIP-CSL) scheme [6]. In this scheme, the evolutions of space- and velocity-integrals of the
distribution function are computed from separate kinetic equations, along with theevolution of the distribution function itself. The implementation guaranties the local
conservation of density. Fig.1 shows
the redistribution of electrons andions. Redistribution is due mainly tothe evolution and interaction ofphase-space structures (phase-space
turbulence) [3].
Of common knowledge is that
accurate simulations require careful
treatment of conserved physicalquantities, such as total mass, totalenergy and total entropy. For
example, to obtain the turbulent
steady-state
accurately,
spurious
Fig.1. Snapshots of the velocity distributions.
heating must be avoided. It is less

known that the conservation properties strongly depend on the simulated physics. With
the COBBLES code for two-species ion-acoustic turbulence, compared to one-species
simulations of the dissipative bump-on-tail instability, the conservation of energy andentropy are degraded by several orders of magnitude. However, the quality of conservationis still satisfying, even in the long-time, fully nonlinear evolution of ion-acoustic turbulence.
Indeed, the relative error in total energy remains below 0.1%.
Let me now focus on the entropy, which is trivially conserved in the physical model.
The Fig.2 (a) shows the time-evolution of the error in electron entropy conservation, in
COBBLES for four different numbers of grid points, and in a Particle-In-Cell (PIC) code,
PICKLES, for two different numbers of particles. The error in total electron entropy
conservation is well below 0.0001%, until phase-space filamentation occurs. Then, it quickly
increases to 15%. Surprisingly, this error does not depend on the number of grid points,
the number of particles, the time-step width (not shown in this figure), or even the type ofnumerical approach (Semi-Lagrangian or PIC). Although I don’t have any clear explanation
for this observation, one hypothesis is that the error in entropy conservation is actually ameasure of infinitely small phase-space structures that are created during the nonlinear
saturation.
A 15% error in entropy(a)conservation is troubling. However, I
didn’t find any impact on the timeevolution of electric field, anomalous
COBBLES,N
xN=
128x1024xv
resistivity, or velocity distribution. As
128x2048
an example, Fig.2 (b) shows the
256x1024
time-evolution of electric field
1024x1024
amplitude. In particular, let me
PICKLES,compare two COBBLES simulations,
Np = 4.2x106
with 256x1024 and 1024x1024 grid
16.8x106
points. The time traces of electricfield amplitude are almost perfectly(b)0.5
overlapping, even though there is a0.4
significant discrepancy in the error in
0.3
entropy conservation (14% and 15%).
0.2[1] H. Schamel, Phys. Scr. 20, 336 (1979)

0.1[2] T.H. Dupree, Phys. Fluids 15, 334(1982)

00.6[3] M. Lesur,

submitted

to

0.7

0.8

0.9

1

PlasmaPhys. Contr. Fusion[4] Y. Kosuga et al., Nucl. Fusion 53,
043008 (2013)[5] M. Lesur, Y. Idomura, and X. Garbet,
Phys. Plasmas 16, 092305 (2009)[6] T. Nakamura et al., J. Comput. Phys.
174, 171 (2001)

Fig.2. Time evolution of (a) error in entropy
conservation and (b) electric field amplitude. Insets:zoom on a smaller timescale. The legend is shared.

ITER Remote experimentation Centre
Takahisa Ozeki
Fusion Research and Development Directorate
Japan Atomic Energy Agency
Plan of ITER remote experimentation center (REC) based on the broader approach (BA)
activity of the joint program of Japan and Europe (EU) is presented. Objectives of REC
activity are 1) to identify the functions and solve the technical issues for the construction of
the REC for ITER at Rokkasho, 2) to develop the remote experiment system and verify thefunctions required for the remote experiment by using the Satellite Tokamak (JT-60SA)
facilities in order to make the future experiments of ITER and JT-60SA effectively, and 3) totest the functions of REC and demonstrate the total system by using JT-60SA and existingother facilities in EU. Preliminary identified items to be developed are 1) Functions of theremote experiment system, such as setting of experiment parameters, shot scheduling, realtime data streaming, communication by video-conference between the remote-site and on-site,
2) Effective data transfer system that is capable of fast transfer of the huge amount of data
between on-site and off-site and the network connecting the REC system, 3) Storage systemthat can store/access the huge amount of data, including database management, 4) Data
analysis software for the data viewing of the diagnostic data on the storage system, 5)
Numerical simulation for preparation and estimation of the shot performance and the analysisof the plasma shot. Detailed specifications of the above items will be discussed and thesystem will be made in these three years in collaboration with facilities of JT-60SA and EU
tokamak, experts of informatics, activities of plasma simulation and ITER. Finally, the total
system of REC will be demonstrated by the middle of 2017.
Remote experiment system: The concept of shot schedule of remote experiment and therequired overall function are identified. The sequence of the procedure of the remoteparticipation will be like that in Fig. 1. The right side of the figure is the main site of ITER in

Fig.1. Concept of the shot schedule of the remote experimentCadarache (or JT-60SA in Naka), and the left side of the figure is a remote site (Rokkasho).
Shot parameter prepared by the on-site researcher will be sent to the shot schedule storage,
and the experiment will be executed by using this shot parameter after the validation checkfrom the viewpoints of the safety and the operation of the facility. After the execution ofexperiment, the measured data are acquired and stored. Remote-site researcher will prepare
the shot parameters in the same manner as the on-site researcher, and the shot schedule will
be sent to the same shot schedule storage through a remote experiment server. To prepare theshot schedule, the remote-site researcher will first access the data storage of the remote site or
the main site, analyze the previous shot data and then prepare the next shot parameters withthe aid of simulations. On the remote site, it is important to perform fusion experiments under
the same condition as the on-site researchers. Main functions to achieve this are preparationand setting of shot parameters, viewing the status of control data, streaming of the plasma
status, data-exchange function of shot events, and monitoring of the facility operation. In
addition, for the remote-site researchers, some tool that provides the atmosphere of being at
the control room of the real facility is needed. Also, we need tools for smooth communication
between the remotes and main sites.
Fast data transfer and Storage: In ITER, the amount of data to be experimentally measured
can be expected to increase up to more than 1TB per discharge, and the amount of theaccumulated data will be an order of 1PB in a year. Therefore, the fast data transfer and thelarge storage are critical issues in the remote site. For the network connection, broadband
network system with larger than 10 Gbps is necessary for REC. The network in Rokkasho site
has been connected to the nearest data center (Hirosaki data center) of SINET4 by 10 Gbps
line in 2012 in collaboration with National Institute of Informatics (NII). Upgrade of the
network bandwidth would be expected to more than 100 Gbps in future. To obtain the fast
data transfer, investigations of effective data transfer method, development of accelerators, etc.
are needed. For the data storage system, high reliability (redundancy), high accessibility(high-speed storage network), extendablity (scalability) and easy maintenance are required.
However, for analysis and reference of experimental data from the remote site, it is notefficient to access data of the on-site storage each time. Therefore, the storage system on the
remote site is necessary for efficient data analyses. Also, investigations of data storage
method and data management in order to find the best way of data mining from the hugeamount of data are also needed.
Simulation and Data analysis: In experiments of large facility such as ITER, the cost of per
discharge is high, and the machine-time must be shared by many researchers. Therefore, inorder to make effective use of the limited machine time, it is important to improve theenvironment for the shot preparation and the data analyses. This can be performed by detailed
data analysis and more efficient plasma performance prediction by the simulation. For thedata analysis software, analyses of experimental data should be done from the remote siteusing the on-site software and also the remote-site software. Functions of the data viewing ofplasma diagnostics, statistical analyses of diagnostic data and mapping of diagnostic data arenecessary. Front-end program to prepare the detailed and complicated data analysis andsimulation codes, i.e. interface to the plasma analyses codes, such as the high energy particle
phenomena, MHD stability, plasma turbulence, etc. will be needed. Numerical simulation isuseful for the experiment. Simple simulation code of the plasma discharge, such as plasma
current, plasma shape, etc. for the pulse preparation is needed as a first step. In order toevaluate a plasma performance, the integrated simulation including the plasma model of thetransport, MHD, divertor, etc. will be needed in the next step. To realize the simulation, it isnecessary to collaborate with the activities of the integrated code, etc.