Proceedings of ITC/ISHW2007

Development of hierarchy-integrated simulation code for toroidal
helical plasmas, TASK/3D
Masahiko SATO, Shinichiro TODA, Yuji NAKAMURA1) , Kiyomasa WATANABE,Atsushi FUKUYAMA2) , Sadayoshi MURAKAMI2) , Masayuki YOKOYAMA, Hisamichi FUNABA,
Hiroshi YAMADA and Noriyoshi NAKAJIMA
1)National Institute for Fusion Science, Toki 509-5292, Japan
Graduate School of Energy Science, Kyoto University, Kyoto 611-0011, Japan
2)
Graduate School of Engineering, Kyoto 606-8501, Japan

The present status of the development of the hierarchy-integrated simulation code for troidal helical plasmas,
TASK/3D, is reported. The TASK/3D is based on the integrated modeling code for tokamak plasmas, TASK. Inorder to extend the TASK for two dimensional conﬁguration to the TASK/3D for three dimensional conﬁguration,
new modules for the radial electric ﬁeld and rotaional transform for a general toroidal conﬁguration have beendeveloped. Numerical simulations for the time evolution of temperature and electric ﬁeld are done for an LHD
experiment by the combination of the diﬀusive transport module and the electric ﬁeld module. Improvement ofthe rotational transform module is also reported.
Keywords: integrated simulation, helical plasma, radial electric ﬁeld, plasma current

1Introduction

In order to systematically clarify conﬁnement physics
in toroidal magnetic conﬁnement systems, a hierarchyrenormalized simulation concept is being developed under domestic and international collaborations with universities and institutes. The hierarchy-renormalized simulation model in toroidal magnetic conﬁnement systems consists of a hierarchy-integrated simulation approach and ahierarchy-extended simulation approach. The former approach, which is mainly based on a transport simulation
combining various simpliﬁed models describing physicalprocesses in diﬀerent hierarchies, is suitable for investigating whole temporal behavior of experimentally observed
macroscopic physics quantities, and the latter approach,which includes ﬂuid core plasma description, kinetic coreplasma description, and peripheral ﬂuid/kinetic description, is focused on the description of mutual interaction
among neighboring hierarchies in a more rigorous way.
The hierarchy-integrated simulation code (TASK/3D) isbased on the integrated modeling code for tokamak plasmas, TASK[1], developed in Kyoto University. In order toextend the TASK code developed for two dimensional conﬁguration to three dimensional, the transport equations forthe rotational transform and the radial electric ﬁeld havebeen reformulated in a general toroidal conﬁguration. With
this new formulation, temporal evolution of the net current
in LHD has been analyzed [2,3]. In this research, present
status of the ﬁrst-stage development of TASK/3D is reported.
This paper is organized as follows. In section
2, numerical model equations and module structure of
author’s e-mail: masahiko@nifs.ac.jp

TASK/3D are described. In section 3.1, we show numerical simulation results obtained by the combination of thediﬀusive transport module and the electric ﬁeld module.
where the proﬁle of the radial electric ﬁeld is determinedby the ambipolar condition. In section 3.2, we show simulation results of time evolution of current proﬁle by using
rotational transform module. The improvement of the rotational transform module is also reported. Finally, section
4 is devoted to summary.

2

Numerical model in TASK/3D

In order to extension the TASK for two dimensional conﬁguration to three dimensional, we are developing andadding new modules into TASK as indicated by red character in Fig.1. As a ﬁrst step, we will carry out simulations by using the diﬀusive transport module TR, rotational transform module EI and electric ﬁeld module ER.
In TASK/3D, the TR module is used for solving the following particle transport equation and heat transport equation;
1 ∂(n s V ′ )V ′ ∂t∂ ( ′V 〈|∇ρ|〉n s V s∂ρ)∂n s+ S s,−V ′ 〈|∇ρ|2 〉D s∂ρ)
(1 ∂ 3′5/3nTVs sV ′ 5/3 ∂t 2(1 ∂3= − ′V ′ 〈|∇ρ|〉 n s T s VEsV ∂ρ2)
∂T s−V ′ 〈|∇ρ|2 〉n s χ s+ Ps ,∂ρ
=−(1)

(2)

Proceedings of ITC/ISHW2007
where n is density、T is temperature、ρ is normalized radial coordinate、V is volume enclosed by a magnetic surface, s expresses species of particles, and the prime denotes the derivative with respect to ρ. 〈〉 is average value
on magnetic surface. The rotational transform module EI
can be used for analyzing inductive plasma current, where
the equation of the rotational transform is solved;()−1∂ΦT∂ι ∂ΦT∂ι
=∂t∂ρ
∂ρ ∂t
 (()−1)−2∂ΦT∂ ∂V ∂ΦT+η ∥ ∂ρ∂ρ∂ρ ∂ρ
{∂〈B2 〉 [(S 11 ι + S 12 )Φ′T ]∂ρ}}∂ΦT ∂P∂ΦT+(S 11 ι + S 12 ) −〈J s · B〉 ,(3)∂ρ ∂ρ∂ρ

steady state of eq.(4). In the electric ﬁeld module, replacing the ambipolar condition to the equation of time evolution for electric ﬁeld is a future work. On the other hand,
the rotational transform module EI has not been added toTASK/3D yet, although temporal evolution of the plasma
current in LHD has been analyzed [3] by using the combination of the EI module, the VMEC MHD equilibrium
module, and BSC/FIT module, where BSC/FIT modules
are used for calculating non-inductive current. Adding theEI module into TASK/3D is underway.

where ΦT is toroidal magnetic ﬂux, η∥ is resistivity, J s isnoninductive current, B is magnetic ﬁeld and P is plasmapressure. The elements of susceptance matrix, S 11 and S 12are given by〈 〉V ′ gθθS 11 =(4)4π2 g〈〉
V ′ gθζS 12 =.
(5)
4π2 g
In eq.(3) 〈B2 〉 and P′ are entered instead of other elementsof susceptance matrix, S 21 and S 22 , by using the MHD
equilibrium condition. In helical plasmas, oﬀ diagonal element of susceptance matrix, S 12 appears due to nonaxisymetry.
For electric ﬁeld module ER, the equation of the electrostatic potential Φ0 is solved;
ϵ0 ϵr 〈|∇ψ|2 〉

∂ dΦ0
∂t dψ

=∑∑

〈B x · Sa1 〉a=e,i+e f 〈n f 1 u f 1 · ∇ψ〉,

3

Numerical Results

3.1 Radial electric ﬁeld module〈B x · ∇ · Πa1 〉

a=e,i−Fig. 1 Module structure of TASK/3D.

(6)

where ψ = ΦT /2π. The second term and the third term ofR.H.S. are source term and fast particle term due to NBI
etc., respectivly. A relative permittivity ϵr is given by〈〉 〈〉c2 (g2 )2〈g2 〉2  |∇ψ|2, (7)ϵr = 1 ++− 2 〈|∇ψ|2 〉  v2Av2A〈vA 〉 where vA is Alfvén velocity. g2 , which describes PﬁrschSchluter current, is given by( )(g )12B · ∇ 2 = B × ∇ψ · ∇ 2 .
(8)BBIn the current stage of the development of theTASK/3D, the electric ﬁeld module has been added to theTASK/3D, where the radial electric ﬁeld is determinedby ambipolar condition which is obtained by considering

In order to check the combination of TR module and electric ﬁeld module, transport simulations are done for anLHD plasma, where the plasma is heated by NBI and thesteady state is achieved. In the simulations, the time evolution of the temperature proﬁle is calculated by using theheat transport equation (2) in the TR module and the radial
electric ﬁeld is determined by ambipolar condition in theER module. The density and rotational transform proﬁles
are ﬁxed. For neoclassical transport model, the single helicity model proposed by Shaing[4] is used. We carried out
simulations for two diﬀerent anomalous transport models,a constant model and a parabolic model. Figure 2 showstime evolution of electron and ion temperatures obtained
from the simulation, where the anomalous transport coeﬃcient is chosen as χa = 1 + ρ2 (m2 /s). In this simulation, theproﬁles of the density and rotational transform are ﬁxed.
The density proﬁle is also plotted in Fig.2. The positiveelectric ﬁeld initially appears around ρ = 0.6 in Fig.3. Thisis mainly due to the positive gradient of the density. As theelectron temperature decreases, the positive electric ﬁeld

Proceedings of ITC/ISHW2007

2.0
1.0m( ytisned2.53-1.5))Vek( eT erutarepmet nortcele

3.0Te (t = 0s)t = 1.0s

100-10-20
0.0

0.2

0.4

0.6

0.8

1.0Fig. 3 Time evolution of radial proﬁle of electric ﬁeld corrseponding to Fig.2.

0.5

Te (t = 0.1 s)

t = 0.2s

minor radius

1.5

1.00s

t = 0.1s

dleif cirtcele

2.0x1019

t =rE

3.5(A)

20)m/Vk(disappears. We also carried out the simulation for the constant anomalous model, where the anomalous coeﬃcientis chosen as χa = 1(m2 /s). The behavior of teperatureand electric ﬁeld for χa = 1(m2 /s) is qualitatively same as
the case for χa = 1 + ρ2 (m2 /s) The dependence of temporal evolution of proﬁles on the transport model is currently
under investigation.

Te (t = 0.2 s)
Te (t = 1.0 s)

0.5

density

0.0

0.0

0.0

0.2

0.4

0.6

0.8

1.0

minor radius
3.5(B)

2.0x10

19

3.0

1.0

3-1.5m(ytisned2.0))Vek( iT erutarepmet noi

1.5

2.5

1.0

Ti (t = 0s)

0.5

Ti (t = 0.1 s)
Ti (t = 0.2 s)

0.5

0.00.0Ti (t = 1.0 s)density

0.2

0.4

0.6

0.8

0.0
1.0

minor radiusFig. 2 Time evolution of radial proﬁles of (A) electron and (B)ion temprature calcurated by the combination of the diffusive transport module TR and the electric ﬁeld module ER. The anomalous coeﬃcent is chosen as χa =1 + ρ2 (m2 /s). The proﬁles of the temperature and density observed in an LHD experiment are used for initial
proﬁles in the simulation. In this simulation , the density
proﬁle indicated by red line is ﬁxed.

3.2

with 2 neutral beam (NB) injection. The 1st NB is injectedto the counter-direction from t ∼ 0.7s and the 2nd NB to
the co direction from t ∼ 1.2s. The blue line denotes theobserved net toroidal current by Rogowski coils. Numerical simulations are done for the LHD neutral beam heatedplasma shown in Fig.6. In this simulation, non-inductivecurrent (bootstrap current and Ohkawa current) is calculated by BSC/FIT module. We use the time evolution of
the plasma density and temperature proﬁles observed in
the experiment. The equilibrium proﬁle is calculated by
the VMEC code at some time interval which is longer than
the time interval for the EI module.

Rotational transform module

The rotational transform module EI can be used for analyzing the plasma current in helical plasmas. As is shown in
refs.[2] and [3], Fig.6 shows a waveform of the discharge

The earlier version of the EI module uses the squareroot of the normalized toroidal magnetic ﬂux as radial coordinate ρ and the rotational transform is discretized at halfmesh points. In this case, the calculated noninductive current becomes negative around r = 0.9 at t=4.5s as shown
Fig.7(A). The problem is mainly due to the disagreement
between the plasma boundary and the half mesh point. Inorder to solve this problem, the EI module has been updated. The new version of the EI module uses the normalized toroidal magnetic ﬂux as radial coordinate s and therotational transform is discretized at full mesh points. Figure 7(B) shows numerical results obtained by the updated
EI module. In this case, the noninductive current does not
become negative around r = 0.9 at t=4.5s.
It is important to decide appropriate time interval forrecalculation of equilibrium proﬁle by VMEC. We are investigating the dependence of numerical results on the timeinterval for recalculating equilibrium proﬁle.

Proceedings of ITC/ISHW2007
1.0s

60

1.5

4.5s

㧔#Ohkawa
Ohkawa+BSObservation

40

120400

,$ 0+Ǳ

1Co-NB0200

4.5s

0Cntr.-NB-200

0.5-1-400246time(s)81.0s1000Fig. 4 Waveform of the discharge with 2 NB-injection. blue isobserved current, red is calculated Ohkawa and bootstrap
current , green is calculated Ohkawa current.Dotted line
and dashed line are absorbed beam power.

-2000.2

0.4

0.6

1rho
1.5

6 1054 1052 105

㧔$,$ 0+Ǳ

140.8

Summary
00.5

The hierarchy-integrated simulation code for toroidal helical plasmas, TASK/3D, is being developed by based on
TASK for two dimensional conﬁguration. In the current
stage, radial electric ﬁeld module has been combined with
the TR module, where the electric ﬁeld is determined by
the ambipolar condition and the single helicity model proposed by Shaing is chosen as neoclassical transport model.
In order to check the electric ﬁeld module ER, we carried out numerical simulations for an LHD experiment, assuming anomalous transport coeﬃcient is a constant or aparabolic function. In the simulations, the density and rotational transform proﬁles are ﬁxed.
For analyzing temporal evolution of the current proﬁlein LHD, numerical simulations are done for an LHD experiment by using the combination of EI, VMEC, BSC/FITmodules, where experimental data is used for time evolution of the density and temperature proﬁles. By using fullmesh points for discretization points, the updated EI module allows us to carry out reliable simulations.
In the next step, we will combine the EI modulewith the TASK/3D. Then, we can calculate temperature,density, radial electric ﬁeld and rotational transform, simultaneously. Next, for neoclassical transport model,
DCOM/NNW[5] will be added to TASK/3D, which isa neoclassical transport database for LHD plasmas constructed using the neural network method. Next, the ambipolar condition for determining electric ﬁeld is replaced
to the equation of time evolution of electric ﬁeld, eq.(4).
In this stage, nonlinearity of transport equations becomesstronger so that the present solver in TASK may not besuitable for solving such equations with strong nonlinearity. Hence, new stable and fast solver will be needed.
For solving a set of nonlinear equations, Newton method
is generally used. In the Newton iteration procedure, itis needed to solve large matrix problem, where the CPU
time is consumed. We are trying to solve the problem bydomain-decomposition method. After the stable and fastPGVEWTTGPV-2 10VU5000.20.4

0.6

0.8

1rho

Fig. 5 Numerical results obtained from (A) earlier version of EImodule and (B) updated version of EI module. In (A),solid lines are rotational transform proﬁles and dashedlines are non-inductive current proﬁles. The blue corresponds to result at t = 1.0s and the red at t = 4.5s. In(B), the red solid line is rotational transform proﬁle att = 4.5s. The green dashed line and the blue solid lineare non-inductive current proﬁle and net current proﬁleat t = 4.5s, respectively.

Newton solver is completed to develop, impurity moduleand heating modules will be added to TASK/3D in orderto further develop the hierarchy-integrated simulation codefor toroidal helical plasmas.
[1] A.Fukuyama et al., Proc. of 20th IAEA Fusion Energy Conf.
(Villamoura, Portugal, 2004) IAEA-CSP-25/CD/TH/P2-3.
[2] Y. Nakamura et al., Fusion Sci. and Tech. 50 (2006) pp.457463.
[3] Y. Nakamura et al., Proc. of 21th IAEA Fusion Energy Conf.
(Chengdu, China, 2006) IAEA-CN-149/TH/P7-1.
[4] K. C. Shaing, Phys. Fluids 27 (1984) 1567.
[5] A. Wakasa et al., Jan. J. Appl. Phys. 46 (2007) 1157.