24D04
Comparison of turbulent transport models in tokamak plasmas
トカマクプラズマにおける乱流輸送モデルの比較
Tomoaki IKARI, Atsushi FUKUYAMA碇 知朗，福山 淳
Department of Nuclear Engineering, Kyoto University,
Yoshida-Honmachi, Sakyo-Ku, Kyoto 606-8501, Japan
京都大学工学研究科原子核工学専攻 〒 606-8501 京都市左京区吉田本町

In tokamak plasmas, the turbulent transport that causes anomalous transport is a non-linear phenomena.
The strong non-linearity of turbulent transport models sometimes causes numerical instability in transport
simulation as a stiff problem. In this case, we have to take a very small time step for stable calculationand the simulation requires very long computation time. We have introduced a numerical scheme proposedby Pereverzev and show that the convergence is strongly improved for a simple diffusion equation with a
non-linear diffusion equation. This scheme is implemented into the transport simulation code TASK/TR
and simulation results for several transport models are compared with experimental observations.

1 Introduction

In tokamak plasmas the transport phenomena isclassified into two categories, neoclassical transport
and anomalous transport. Neoclassical transport iscaused by coulomb collisions and magnetic mirror due
to the tokamak torus configuration. Many experiments
have shown, however, that the transport exceeds theamount which is predicted by the neoclassical theory.
This increased transport is called anomalous transportand dominant in actual tokamak devices. Anomaloustransport is considered to be caused by turbulent transport, which can be driven by the microscopic electmagnetic fluctuations coming from pressure and current density gradients. Turbulence transport have notbeen figured out fully, and a number of models describe it have been proposed.
The turbulent transport that causes anomalous transport is a non-linear phenomena. The strong nonlinearity of turbulent transport models sometimes
causes numerical instability in transport simulation asa stiff problem. In this case, we have to take a verysmall time step for stable calculation and the simulation requires very long computation time. Pereverzev
has proposed [1] a numerical scheme strongly improves the numerical convergence.
In our simulation of tokamak plasmas we use TR
module of integrated transport analysis code TASK[2]. The TASK/TR module assumes 1-D diffusivetransport model on the ground of flux surface average. The purpose of this paper is to introduce the improved numerical scheme into the TASK/TR moduleand compare several transport models including stiff
ones with experimentally observed profiles.

2 Basic Equation

The basic equations of TASK/TR are following

three diffusive equations:∂∂(ns )V ! = − (V ! Γs ) + SsV ! ,(1)
∂t
∂ρ!
"
∂ 3
∂!5/3
ns TsV= −V !2/3
(V ! Qs ) + SEsV !5/3 ,∂t 2∂ρ(2)!
"#!$%"I∂ ∂ψ∂ η#∂|∇ρ |2 ∂ ψ=∂t ∂ρ∂ ρ µ0 V ! $R−2 % ∂ ρR2∂ρ&η#−$(JCD + JBS )Bφ % ,(3)I$R−2 %
Γs = $|∇ρ |%nsVs − $|∇ρ |2 %Ds

∂ ns,
∂ρ(4)

∂ (ns Ts )Qs = VEs ns Ts − $|∇ρ |%χs∂ρ
#!
"&5∂ ns+ − $|∇ρ |%Ds − χs Ts, (5)2∂ρwhere s denotes the particle species, ns is the particle
density, Ts the temperature , Γs the particle flux, Qs
the heat flux, Ds the particle diffusivity, χs the thermaldiffusivity, VEs the sum of the particle pinch velocity
Vs and the heat pinch velocity, ψ the poloidal magnetic flux, η# the parallel neoclassical resistivity, µ0
the permeability of free space, I = Bφ R where Bφ thetoroidal magnetic field, R is the major radius, jCD andJBS are the induced current density and the bootstrap
current density, V ! the derivatives of the volume enclosed by a flux surface with respect to the normalizedminor radius ρ and $|∇ρ |% and $|∇ρ |2 % are geometricquantities.

Some turbulent transport model such as GLF23 [3]are numerically stiff and when using these models intransport simulations we have to decrease the size oftime step considerably, which is very time-consuming.
We, however, can avoid this problem if the virtual diffusive term is added to transport equation.
We explain this concept using a simple diffusion
equation in cylindrical coordinates!
"∂u 1 ∂∂u=
ρD+ S.
(6)∂tρ ∂ρ
∂ρ

Flux q

3 Stable numerical schemeq = D0 η + D1 (η - ηcr)

,q0q = Deff (η 0 ) ηq∞η η η∞η, η=-uIteration

x0 10Fig.1 shows the gradient-flux relation of a typical stiff
transport model. The solid line shows that the diffusive coefficient and the flux increase abruptly when the Figure 1: The gradient-flux relation of a typical stiff
amplitude of the gradient exceeds the threshold ηcr . transport model. The dash line shows the relation ofTransport models evaluate the diffusion coefficients as improved equation.
Deff illustrated by the dotted line in Fig.1. Provided
that the gradient is η0 at a certain time step, the transport model returns the value Deff (η0 ). If the size of100
Modifiedtime step in calculation is not sufficiently small, then!
Original
80
the transport code evaluates the gradient as η0 , andnext the transport model returns much higher Deff than60
previous value. This leads to numerical instability andthe calculation doesn’t come to converge numerically.
40
Now we consider the modified equation20!
"∂u 1 ∂∂u=ρ D̃− ρ V̄ u + S,
(7)0∂tρ ∂ρ∂ρ
-2-1
010 -510-4

10-3

101010where D̃ = Deff + D̄，V̄ = D̄(∂ u/∂ ρ )/u, and D̄ is thearbitrary virtual diffusive coefficient. Eq.(7) is idenFigure 2: Comparison of the number of iteration.
tical to Eq.(6) mathematically, but the non-vanishingterm remains if D̃ is calculated at t + τ implicitly whileV̄ is calculated at t, where t is certain time step and τ
is the size of time step. The gradient-flux relation of 5 Comparison of Transport Models
Eq.(7) is shown by the dash line in Fig.1, and obviWe will present the simulation result of actual tokaously the instability will not be caused and the calcu- mak plasmas with several turbulent transport models
lation will converge straightforwardly.
using improved TASK/TR.
Δt

4 Numerical Results

Acknowledgment

This work is supported by Grant-in-Aid for Scientific
We calculated Eq.(6) and Eq.(7) as preliminary calResearch (S) (20226017) from JSPS
culation by means of finite difference method. The
transport model is assumed to be stiff as is shown
Fig.1, and we set D0 = 0.01，D1 = 0.1，S = 1.0， ReferencesD̄ = 0.05. Then we found that the iteration using [1] G.V. Pereverzev, G.Corrigan, Computer Physics
Eq.(6) does not converge at the size of time step larger
Communication 179 (2008) 579-585
than ∆t = 10−2 , while that using Eq.(7) does converge [2] A.Fukuyama et al., Proc. of 20th IAEA Fusion
even when the size of time step is ∆t = 100 . Fig.2
Energy Conf. (Villamoura Portugal,2004)shows the comparison of the dependence of the numIAEA-CSP-25/CD/TH/P2-3
ber of iterations on the size of time step, from which [3] R.E.Waltz et al., Phys. Plasmas 4 (1997) 2482
we found that this scheme is validated in the preliminary calculation.