Plasma and Fusion Research: Regular Articles

Volume 14, 3403127 (2019)

First Application of 3D Peripheral Plasma Transport Code
EMC3-EIRENE to Heliotron J∗)
Ryota MATOIKE, Gakushi KAWAMURA2,3) , Shinsuke OHSHIMA1) , Masahiro KOBAYASHI2,3) ,Yasuhiro SUZUKI2,3) , Kazunobu NAGASAKI1) , Suguru MASUZAKI2) , Shinji KOBAYASHI1) ,Satoshi YAMAMOTO1) , Shinichiro KADO1) , Takashi MINAMI1) , Hiroyuki OKADA1) ,Shigeru KONOSHIMA1) , Toru MIZUUCHI1) , Hirohiko TANAKA4) , Hiroto MATSUURA5) ,Yuhe FENG6) and Heinke FRERICHS7)
Graduate School of Energy Science, Kyoto University, Gokasyo, Uji 611-0011, Japan
1)
Institute of Advanced Energy, Kyoto University, Uji 611-0011, Japan
2)
National Institute for Fusion Science, National Institutes of Natural Sciences, Toki 509-5292, Japan3)
Department of Fusion Science, Graduate University for Advanced Studies, Toki 509-5292, Japan
4)
Graduate School of Engineering, Nagoya University, Nagoya 464-8603, Japan
5)
Graduate School of Engineering, Osaka Prefecture University, Osaka 599-8531, Japan
6)
Max-Planck Institute for Plasma Physics, Greifswald, Germany
7)
University of Wisconsin-Madison, Wisconsin, U.S.A.
(Received 9 January 2019 / Accepted 5 April 2019)
The 3D peripheral plasma and neutral transport code, EMC3-EIRENE was applied to the Heliotron J witha wide and flexible controllability of magnetic configuration. This code requires a three-dimensional (3D) gridwith high resolution in the peripheral plasma region to reproduce the fine plasma structure. The grid generation
tool, FLARE, was utilized to create the grid in conjunction with the code developed to arbitrarily set the outerboundary of the peripheral grid. After setting up the 3D grid, we carried out the EMC3-EIRENE calculation
successfully for the first time in the Heliotron J’s standard configuration. In addition, the convergence of theiterative calculation and the eﬀects of diﬀerent grid resolutions upon the calculation were investigated. A good
numerical convergence was obtained, and the influence of resolution was observed in the electron density in thedivertor region.
c 2019 The Japan Society of Plasma Science and Nuclear Fusion ResearchKeywords: EMC3-EIRENE, Scrape-Oﬀ Layer, divertor, transport, modeling, Heliotron J
DOI: 10.1585/pfr.14.3403127

1. Introduction
Design of divertor is an essentially important issue tocontrol extremely high heat and particle fluxes in the divertor of fusion reactors [1]. For instance, at ITER, the heatload is estimated to be 10 MWm−2 under steady-state discharge [2]. Many eﬀorts have been dedicated to reducing
heat and particle fluxes to the divertor [3, 4], and magnetic
field topology is one of the keys for controlling the heatload to the surfaces. The heat and particle deposition profiles depend strongly upon the connection length distributions of the magnetic field on the divertor plates [5]. Helical devices inherently have a three-dimensional (3D) magnetic field configuration, and the divertor has a 3D structure as well. The 3D eﬀects of the peripheral magnetic
field are also discussed in tokamak devices by introducing
resonant magnetic perturbation coils to control the edge localized modes [6,7]. Therefore, analysis of the 3D eﬀect ofauthor’s e-mail: matoike.ryota.62c@st.kyoto-u.ac.jp∗) This article is based on the presentation at the 27th International Toki
Conference (ITC27) & the 13th Asia Pacific Plasma Theory Conference(APPTC2018).

magnetic field structure and modeling of peripheral plasmatransport is needed. The peripheral transport code, EMC3EIRENE [8, 9], has been used to model scrape-oﬀ layer(SOL) -divertor plasmas and plasma-surface interactions
on stellarators including W7-X [10–12] and LHD [13–17],and has also been applied to tokamak devices [18,19]. Heliotron J [20] is a helical-axis heliotron device with a highcontrollability of configuration. The dependence of electron temperature, density and heat/particle flux upon themagnetic configuration of divertor plasmas has been discussed using the divertor probe array (DPA) [21–23]; however, no comparison between experiment and simulation of
SOL plasma has yet been performed. In this study, we apply the EMC3-EIRENE code to Heliotron J for 3D modeling of peripheral plasma for the first time. The grid generation procedure, typical calculation results, and dependenceof calculation results upon diﬀerent spatial grid resolutions
are discussed.

3403127-1

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

Plasma and Fusion Research: Regular Articles

Volume 14, 3403127 (2019)

2. Numerical Modeling of Heliotron J
SOL Plasma
2.1

EMC3-EIRENE code

The peripheral transport code EMC3-EIRENE is used
to model the Heliotron J plasma. This code solves the Braginskii’s fluid equations below using a Monte Carlo (MC)method in a 3D grid with a field-aligned structure. The
conservation of mass, conservation of momentum, conservation of energy for electrons and ions respectively read∇ · (nv ) + ∇⊥ · (−D∇⊥ n) = S p ,(1)∇ · (mi nv v − η ∇ v )+∇⊥ · (−mi v D∇⊥ n − η⊥ ∇⊥ v ) = ∇ p + S m ,5∇ · − κe ∇ T e + nT e v25+∇⊥ · − χe n∇⊥ T e − T e D∇⊥ n
2= −k(T e − T i ) + S ee ,5∇ · − κi ∇ T i + nT i v2
5+∇⊥ · − χi n∇⊥ T i − T i D∇⊥ n
2= k(T e − T i ) + S ei ,

(2)

(3)

(4)

where p = n(T e + T i ) and η⊥ = mi nD. The variables p,
n, m and T represent pressure, density, mass, and temperature, respectively. The subscripts i and e indicate ions andelectrons, respectively. The subscripts ⊥ and  indicate
the perpendicular and parallel directions of the magnetic
field, respectively. The parallel transport coeﬃcients, η , κeand κi are considered to be classical, while the cross-fieldtransport coeﬃcients, D, χe and χi are anomalous and usually determined according to experiments. Source terms,

S p , S m , S ee and S ei are the particle, momentum and energy source arising from plasma-neutral interactions suchas ionization, excitation and charge exchange. EMC3 isself-consistently coupled with EIRENE, which solves the
kinetic Boltzmann equations for neutral particles. EMC3EIRENE solves a steady-state distribution of plasma andneutral particles through iterative calculations, and each
iteration uses values obtained from the previous iteration
step.

2.2 Generation of the 3D grid system
The magnetic field of Heliotron J has l/m = 1/4 helical periodicity and the poloidal cross-sections at toroidal
angles φ = 0◦ and 45◦ have up-down symmetry. Figure 1
shows contour plots of the connection length on thepoloidal cross-sections. Quantities of the vacuum magnetic field are obtained using the KMAG [24,25] code. The
simulation box size for the EMC3-EIRENE calculation is45◦ in the toroidal direction, i.e., a half helical section, andthe full-torus is modeled eﬀectively using the helical periodicity and an assumption of up-down symmetry at φ = 0◦
and 45◦ .
As shown in Fig. 1 (b), there are seven “points” likethe X-point of a tokamak divertor. A pair of divertor legsconnected to the X-point on (φ, θ) = (0, 0) is in contact withthe chamber wall. The structure of the connection length
in the SOL region is determined by the vacuum wall in Heliotron J, where the wall works like a limiter in tokamaks.
As shown in the contour plot in Fig. 1 (b), magnetic fieldlines in the SOL are connected to a local region in boththe poloidal/toroidal directions (see blue arrows shown in
Fig. 1 (b)). In the other region, the connection length on
the wall is quite small. Note that it seems that only a pairof divertor legs is connected to the wall, but all the magnetic field lines in the divertor legs are also connected toFig. 1 (a) Top view of Heliotron J plasma. (b) Contour plots of connection length on poloidal cross-sections every 22.5◦ . The plot on
φ = 90◦ is equal to φ = 0◦ . The white line represents the chamber wall. The green line on poloidal cross-section φ = 0 represents
the outermost boundary for plasma modeling.

3403127-2

Plasma and Fusion Research: Regular Articles

Volume 14, 3403127 (2019)
Table 1 Cell number of the grid system we prepared. The grid1 is the coarsest and the grid 4 is the finest. The cellnumber in the toroidal direction is 90 in all grids.

CORE
grid 1
grid 2
grid 3grid 4

Fig. 2 The computational grid at φ = 0◦ generated based on
the vacuum magnetic field for the standard configuration.
The bold red line represents the LCFS, and the bold blackline represents the vessel wall.

.
the chamber wall at the same poloidal/toroidal location indicated by the blue arrows in Fig. 1 (b).
The 3D grid of Heliotron J has been generated for the
EMC3-EIRENE calculation; it has four domains: CORE,
EDGE, SOL and VACUUM. These domains are shown in
Fig. 2 and labeled (1) - (4), respectively. The confined region of the plasma is covered by the CORE and EDGE
domains. The boundary between the CORE and EDGE
domains is set at r/a  0.8 and the outer boundary of the
EDGE domain is at the last closed flux surface (LCFS).
The area outside of the LCFS is covered by the SOL and
VACUUM domains. The outer boundary of the SOL domain is defined to cover the region with connection length
Lc > 10m, see Fig. 1 (b). The VACUUM domain is defined
to cover chamber wall. Neutral particle transport is solvedin the entire domains and plasma transport is solved only in
the EDGE and SOL domains. The grid points in these domains have field-aligned structures in the toroidal direction
to reduce aliasing error arising from the finite resolution of
the grids.
In the grid generation procedure, we combined thegrid generator tool FLARE [26] and an additional code developed by ourselves. The 3D grid is generated by the following FLARE code procedures;(a) generate a 2D base grid of the EDGE domain on the
cross section at φ = 0◦ ;(b) expand the outer boundary of the EDGE domain andgenerate a SOL domain of the base grid above;(c) trace field lines in the toroidal direction from each nodeof the base grid and generate a 3D grid;(d) generate the VACUUM and CORE domains.
In the procedure (b), the SOL grid between LCFS, i.e.,the outer boundary of the EDGE domain, and the outer

5
5
5
5

radialEDGE SOL
15
15
15
15

poloidal
VAC

15
45
15
45

10
10
10
10

120
120
360
360

boundary of SOL domain was generated manually by ouradditional code. This is because FLARE code cannot complete the procedure (c), due to a problem when a field line
goes to outside of the prepared magnetic field data when
we cover whole peripheral plasma on the cross section inprocedure (b). The separation of the plasma domains at
the LCFS with our additional code realizes the flux-surfacealigned grid in EDGE domain and the suﬃciently large
SOL domain simultaneously.
We generated multiple grids of the standard configuration of Heliotron J with a variety of resolutions using theprocedure above. In this article, as shown in Table 1, weshow typical four examples with diﬀerent cell numbers in
the radial/poloidal resolutions to check the eﬀect of grid
resolution on calculation results. The toroidal resolution iscommon to all domains. A large cell number means a highresolution in the directions.

3. Results and Discussion
3.1 Typical calculation results
We carried out EMC3-EIRENE simulations with thefollowing parameters; input heating power: P = 200 kW,
electron density upstream (blue dashed line in Fig. 2) of
the plasma grid: ne = 0.8 × 1019 m−3 . The perpendicular transport coeﬃcients are assumed as follows; D =0.5 m2 /s and χ = 2 m2 /s. Figure 3 shows the calculation results of electron density (Fig. 3 (a)), electron temperature (Fig. 3 (b)), Mach number (Fig. 3 (c)), and hydrogen atom density (Fig. 3 (d)). The distributions of electron
density and temperature reflect the structure of the connection length shown in Fig. 1. The Mach number distribution
suggests a poloidal modulation of plasma flow along thedivertor legs. The hydrogen atom density is locally highbecause of recycling by the divertor legs.

3.2 Convergence of calculation
The EMC3-EIRENE code solves the equations iteratively. The number of MC particles is increased step bystep to save calculation time. Figure 4 (a) shows the number of Monte Carlo particles during the iterative calculations in the EMC3-EIRENE operation. In this figure, we
compared two methods to increase the number of MC particles up to 2 million. As shown in Fig. 4 (b), the divertor

3403127-3

Plasma and Fusion Research: Regular Articles

Volume 14, 3403127 (2019)

Fig. 3 Contour plot of (a) electron density, (b) electron temperature, (c) Mach number, and (d) hydrogen atom density on the plane φ = 0◦
with the standard magnetic configuration. The dashed green line represents the LCFS, and the bold black line represents the vesselwall.

Fig. 4 (a) Number of Monte Carlo (MC) particles during iterative calculations. (b) Convergence of electron density
during calculations. The same color of plots represents
the same calculation between (a) and (b).

electron density converged within 1 to 2 hours when thenumber of MC particles was increased step by step. On theother hand, when the number of MC particles is increased
with one step, the convergence of electron density takes 4
or 5 hr. We use the acceleration method with MC particles
increased step by step in the calculations given below.
We checked the convergence of the calculation results
at each calculation step. To check the status of convergence, we plotted (a) relative change of ne , (b) relativechange of T e , (c) downstream ne , and (d) downstream T e at
diﬀerent radial and poloidal resolutions in Fig. 5. Here nds
eand T eds are the pressure-weighted average for the downstream cells just in front of the plasma-facing walls, andhave the following form:np(l) (l) (l)l=1 P V ne
np (l) (l) ,l=1 P Vnp (l) (l) (l)l=1 P V T e
T eds = .
np(l) (l)l=1 P Vndse =Fig. 5 Eﬀect of radial resolution in the SOL domain during iterative calculations about (a) relative change of electron
density, (b) relative change of electron temperature, (c)downstream electron density and (d) downstream electron temperature.

The variable P is local plasma pressure, of the form
(5)

(6)

P = ne (T e + T i ),

(7)
np

V is a volume of a cell. The summation l=1 is takenover downstream cells, where “downstream” is basically3403127-4

Plasma and Fusion Research: Regular Articles

Volume 14, 3403127 (2019)

Fig. 7 (a) connection length, (b) electron density, and (c) electron temperature along the line connecting the whitemarkers in Fig. 6.

Fig. 6 Electron temperature of the φ = 22◦ plane calculated with(a) grid 1, (b) grid 4. The green dashed line indicates the
LCFS.

defined as the outermost cells of the SOL domain grids,
where SOL plasma exists. The variable np indicates thecell number in the poloidal direction. Relative change isthe change ratio compared with previous iterative calculation results.
The relative error apparently decreases with the increase in MC particles in Figs. 5 (a) and (b). Statisticalerror increases slightly with the grid resolution, because
there are fewer MC particles per cell with larger resolutions. Even in the finest resolution case, i.e., grid 4 inTable 1, the relative errors of ne and T e maintain values
ds
below 1% after convergence. Similarly, ndse and T e also
converge to certain values during iterative calculations in
Figs. 5 (c) and (d). These results show that influence of grid
resolution upon the calculation result is small enough practically, for example, compared to the accuracy of experimental measurements. However, the downstream electron
density tends to increase with a finer grid, and it is necessary to choose a suﬃcient resolution carefully to solve
particle transport correctly.

3.3

Comparison with the experimental setup
All calculation results with resolutions shown in Sec.

3.2 converge to certain values properly, and the grid resolution should satisfy the spatial resolution required for thecomparison with experimental measurements.
For example, the electron temperatures of divertor
legs with diﬀerent resolutions are compared in Figs. 6 (a)and (b), which show the electron temperature in a part ofthe φ = 22◦ plane with grids 1 and 4, respectively. It can beseen that a much smoother spatial distribution structure is
obtained in Fig. 6 (b) compared to that in Fig. 6 (a). Whitecross-markers shown in Fig. 6 are the measurement pointsof the divertor probe array of Heliotron J [23].
The distributions of connection length, electron density and electron temperature along the white markers in
Fig. 6 are shown in Fig. 7. The structures of the electron
density and temperature are roughly captured even with thelower resolution grids, without changing the peak value or
spatial gradient. However, smoother and more natural distributions can be seen with the higher resolution grids.

4. Summary
The EMC3-EIRENE code has been successfully applied to Heliotron J. A 3D grid with a field-aligned structure is required for the EMC3-EIRENE calculation. A grid
generation tool, FLARE, is utilized to create the grid in
conjunction with the additional code we developed. We
generated multiple grids with diﬀerent radial and poloidal
resolutions. Calculations converged with all grids, regardless of resolution. The relative change of electron density
and temperature compared with the previous iterative calculation is reduced less than 1%. The detailed structure of
the peripheral plasma is obtained by the fine grid. Elec-

3403127-5

Plasma and Fusion Research: Regular Articles

Volume 14, 3403127 (2019)

tron density is slightly changed by changing grid resolution, and we should explicate the cause of this throughfurther research. However, the eﬀect of grid resolution issmall enough to practically compare with the experimental
measurements. Comparison with measurements should betaken into consideration when discussing the grid resolution.

Acknowledgments
The authors are grateful to the Heliotron J staﬀs
for useful discussion.
This work was supported
partly by JSPS KAKENHI Grant Number 16K18340.
This work was performed with the support and under the auspices of the NIFS Collaborative Research
Program (NIFS10KUHL030, NIFS17KUHL081, NIFS
18KUHL084, NIFS18KNTT047), the Collaboration Program of the Laboratory for Complex Energy Processes,
Institute of Advanced Energy, Kyoto University, Future
Energy Research Association and JSPS Core-to-Core Program, A. Advanced Research Networks.
[1] P.C. Stangeby and G.M. McCracken, Nucl. Fusion 30, 1225
(1990).
[2] A.S. Kukushkin et al., Nucl. Fusion 49, 075008 (2009).

[3] Y. Feng et al., Nucl. Fusion 49, 095002 (2009).
[4] A. Loarte et al., Nucl. Fusion 47, 203 (2007).
[5] S. Masuzaki et al., Nucl. Fusion 42, 750 (2002).
[6] T.E. Evans et al., Nature Physics 2, 419 (2006).
[7] M. Kobayashi et al., Nucl. Fusion 55, 104021 (2015).
[8] Y. Feng et al., Contrib. Plasma Phys. 44, 57 (2004).
[9] D. Reiter et al., Fusion Sci. Technol. 47, 172 (2005).
[10] Y. Feng et al., Nucl. Fusion 56, 126011 (2016).
[11] F. Eﬀenberg et al., Nucl. Fusion 57, 036021 (2017).
[12] H. Frerichs et al., Nucl. Fusion 57, 126022 (2017).
[13] M. Kobayashi et al., Fusion Sci. Technol. 58, 220 (2010).
[14] G. Kawamura et al., Contrib. Plasma Phys. 54, 437 (2014).
[15] G. Kawamura et al., Plasma Phys. Control. Fusion 60,
084005 (2018).
[16] S. Dai et al., Plasma Phys. Control. Fusion 59, 085013(2017).
[17] S. Dai et al., Nucl. Fusion 58, 096024 (2018).
[18] T. Lunt et al., Nucl. Fusion 52, 054013 (2012).
[19] J. Huang et al., Plasma Phys. Control. Fusion 56, 075023(2014).
[20] T. Obiki et al., Nucl. Fusion 41, 833 (2001).
[21] T. Mizuuchi et al., Fusion Sci. Technol. 50, 352 (2006).
[22] T. Mizuuchi et al., J. Nucl. Mater. 363-365, 600 (2007).
[23] T. Mizuuchi et al., J. Nucl. Mater. 313-316, 947 (2003).
[24] Y. Nakamura, J. Plasma Fusion Res. 69, 41 (1993).
[25] M. Wakatani et al., Nucl. Fusion 40, 569 (2000).
[26] H. Frerichs et al., Comput. Phys. Commun. 181, 61 (2010).

3403127-6