22F06
Momentum transport in full-f gyrokinetic simulations
full-fジャイロ運動論シミュレーションにおける運動量輸送	 1

1Yasuhiro Idomura and 2Sebastien Jolliet
井戸村	  泰宏, 2ジョリエ	  セバスチャン

12

Japan Atomic Energy Agency, Rokkasho, Kamikita, Aomori 039-3212, Japan
École Polytechnique Fédérale de Lausanne, CH-1015 Lausanne, Switzerland
1
日本原子力研究開発機構	  〒039-3212	  青森県上北郡六ヶ所村
2
ローザンヌ連邦工科大学	  CH-1015 Lausanne, Switzerland

Momentum transport in full-f gyrokinetic simulations is discussed. Accuracy issues of momentum
transport in gyrokinetic simulations are studied using a full-f gyrokinetic Eulerian code GT5D. Test
calculations in the axisymmetric limit and in turbulent tokamak show that toroidal angular momentum isconserved and that momentum transport and related radial electric fields are correctly computed.
Comparisons of ion temperature gradient driven turbulence with and without adaptive momentum sources,
which cancel a flow generation, show that momentum transport tends to suppress heat transport byenhancing radial electric field shear, which is connected with rotation profiles through a force balance
relation. An existence of residual stress is also identified in the latter simulation.

1. Introduction
Momentum transport induced by tokamak
micro-turbulence is one of key issues in
predicting the performance of ITER, and has beenstudied intensively using gyrokinetic simulations.
However, based on an ordering argument, Parraand Catto [1] raised a serious concern about theaccuracy of momentum transport and related
radial electric fields Er in gyrokinetic simulations,in particular, with so-called full-f approaches,
where both momentum transport and rotation
profiles develop self-consistently. On the otherhand, a recent work [2] clearly demonstrated thattoroidal angular momentum is conserved at anyorder, provided that the equation system isderived based on the modern gyrokinetic theory[3] with keeping an energetic consistency. In thiswork, we address accuracy issues of momentum
transport in a full-f gyrokinetic Eulerian code
GT5D [4] by studying basic properties ofmomentum transport in the axisymmetric limit
and in turbulent tokamak. After clarifying theaccuracy issues, we discuss influences ofmomentum transport on heat transport in the ion
temperature gradient driven (ITG) turbulence.
2. Calculation model
GT5D is based on the modern gyrokinetic theory,in which the gyrokinetic equation is simply given as∂f+ {f, H} = C(f ) + S,∂t(1)

where f is a full-f ion distribution function, { , } andH are the Poisson bracket operator and theHamiltonian in the gyro-center coordinates, C is alinear Fokker-Planck collision operator, and S is asource term, which involves fixed on-axis heating, aheat and momentum sink near the plasma boundary,
and an adaptive momentum source. Turbulent
fluctuations are determined by solving thegyrokinetic Poisson equation with linearized ion
polarization density and adiabatic electron response.
This equation system is derived from a single field
Lagrangian and keeps an energetic consistency that
leads to the energy and momentum conservation [2].
Eq. (1) is implemented using the toroidal angle andnon-dissipative and conservative finite difference,
which is a variant of centered finite difference, andtherefore, the scheme keeps the toroidal symmetry
numerically.
3. Accuracy of momentum transport
The first test case is in the axisymmetric andcollisionless limit. If one initializes plasma with alocal Maxwellian distribution fLM, its nonequilibrium component excites initial electric fields,
while plasma stays quiet when one uses a Vlasov
equilibrium fCM given as a function of the canonicalangular momentum !! . This shows a conservationof !! in the axisymmetric limit. The second test
case is still in the axisymmetric limit but withcollisions, i.e., in the neoclassical case. In this case,the system initialized with fLM shows a collisional
relaxation towards a neoclassical quasi-steady state,where the particle flux vanishes due to theambipolar condition, the heat flux converges toneoclassical levels, and Er is determined to satisfy aforce balance relation,��
Ti Id ln Ti
d ln nieV� =(k − 1)−+ Er , (2)mi Ωi ψ �
drdrTimomentum transport. Numerical experiments witha fixed on-axis heating are performed for
Cyclone-like parameters with plasma sizes of
1/ρ*=a/ρti=100, 150, 225, where a is the minor
radius and ρti is the thermal ion gyro-radius. Figure2 (a) shows radial profiles of the heat diffusivity χi,which shows worth-than-Bohm scaling [6]. In all
the plasma sizes, it is shown that transport levels
observed with the adaptive source (solid lines)
exceed those without the adaptive source (dashedlines). In Fig. 2 (b), it is shown that a spontaneous
flow generation due to residual stress enhances Er
shear through the force balance relation (2), whichis approximately satisfied even with turbulent
fluctuations. The enhancement of Er shear leads tosuppression of heat transport. Another important
finding in the numerical experiment with theadaptive source is the observation of finite
momentum flux induced at V||~0. The result clearlyshows an existence of residual stress, which isdriven by the asymmetry of k|| turbulent spectra.
Although various symmetry breaking mechanisms
such as Er shear [7] and profile shear [8] have been
proposed, the result show similar trend between
radial profiles of momentum flux and Er shear,
which suggests that Er shear is a dominantsymmetry breaking mechanism in the present
numerical experiment.

r s/!GB!!ii/!
GBieE "" /T/TeE
r s iwhich shows a balance among Er, the parallel flow
V||, and the gradients of density ni and temperature
Ti. In this test, a coefficient k is estimated using Er
determined by the gyrokinetic Poisson equation,and other plasma profiles, ni, V||, Ti, and it is foundthat k agrees with the neoclassical theory, whichsuggests that Er is correctly computed when V|| isgiven. We then examine the accuracy of momentum
transport in turbulent tokamak. By taking !!
moment of Eq. (1) and applying gyro-average andflux-surface average operators, one obtains a radialtransport equation or a conservation law of toroidal
angular momentum, which gives a balance among
evolutions of flows, momentum transport, andmomentum sources as shown in Fig. 1. It is notedthat a conservative nature of GT5D, which leads tocancellation of all the terms related to the poloidal
flux ψ, plays an important role in deriving thetoroidal momentum conservation. In the core region,a bipolar flow generation (green) is balanced withmomentum transport (blue), while in the edge
region, momentum transport (blue) is canceled byco-current momentum injection from a sink term
(magenta), which keeps a no-slip boundary
condition. In the whole region, the momentum
19
Plasma size
size scalingscaling ofof avalanche-likeavalanche-like heat
heat transporttransport in
in tokamaks
tokamaks
19
conservation is satisfied with good accuracy. In the Plasma(b) 0.04
1.5(a)
1.5
0.04
final test case, we examine influences of higher
100,ad100,ad100,ad100,ad150,ad150,ad150,ad150,ad0.02
0.02
order drift terms on heat and momentum fluxes by225,ad
225,ad225,ad225,ad1100
implementing higher order Hamiltonian [5]. The
0.5
0.5
−0.02
−0.02
result shows that influences of higher order drift
0−0.040−0.04
terms are negligibly small, as is expected from the00
0.2
0.4
0.6
0.8
1100
0.2
0.4
0.6
0.8
110.2
0.4
0.6
0.8
0.2
0.4
0.6
0.8
r/ar/ar/ar/agyrokinetic ordering. All the above test results
88 Radial profiles of (a) ion heat
0.2
0.2 diffusivity χ andFig.2(b)100,ad100,adi
suggest the correctness of turbulent momentum* 150,ad150,ad
66 electric field E in ITG turbulence with 1/ρ =100,radialr225,ad225,ad
0.1
0.1
transport in GT5D.

<!tRv||b"f>gf (A.U.)

100
00

0.5
0error<!tRv||b"f>gf<!V#">gf$<Rv||b"S>gf

-0.5
-1
0.2

0.4

0.6

0.8

r/a

Fig.1 Momentum conservation in ITG turbulence.

4. Impact of momentum transport on heattransport
In order to study influences of momentum
transport on turbulent transport, we compare ITG
turbulence with and without an adaptive momentum
source, which cancels a flow generation due to// ti

V /v/v
V// tiR0/L/LTi
R0 Ti

150,44 225. Solid (dashed) lines show results with100,ad100,ad(without)the adaptive
momentum00 source.
22150,ad150,ad
0.2
0.2

0.4
0.4

225,ad225,ad
0.6
0.6
r/ar/a

0.8

1

−0.1
−0.1
00.2

0.4

0.6

0.8

1

0.8
1
0
0.2
0.4
0.6
0.8
1
r/ar/aReferences[1] F. I. Parra and P. J. Catto, Plasma Phys. Control.
FigureFigure 8.
8. Radial
Radial profile
profile ofof χχi /χ(top, left),Rleft),R00/L/LTTii (bottom,(bottom, left),left), eρsE
eρsErr/T/Ti i (top,(top,i /χGB
GB (top,Fusion
50,(2008).
right)right)right) andand VV!065014/vtiti (bottom,(bottom,right) forfor thethe ρρ∗∗ == 1/100
1/100 (blue),(blue), ρρ∗∗ == 1/150
1/150 (red)(red) andand!/v∗∗ρρScott== 1/225
withmomentumsource.
are
1/225 (green)(green) simulations
simulations
with adaptive
adaptive
momentumsource. All
All profiles
profilesare[2] B. averaged
andJ.
Smirnov,
Phys.
Plasmas17,
averaged between
between tt11 == 1300a/c
1300a/css andand tt22 == 5130a/c
5130a/css.. OnOn thethe figures
figures ofof χχii andand EErr
112302(2010). aregraph,
dashedgraph,thethedashed lines
lines
are thethe simulations
simulationswithout
without thethe adaptive
adaptive momentum
momentum source.
source.
[3] A. Brizard and T. S. Hahm, Rev. Mod. Phys. 79, 421
AsthoseAs anan example,
example,those terms
terms areare displayed
displayed onon thethe left
left plot
plot ofof figurefigure 99 forfor thethe ρρ∗∗ == 1/150
1/150(2007).
simulation
simulation withwith andand without
without adaptive
adaptive source.
source. The
The FBEFBE describes
describes mostly
mostly aa balancebalance[4] theY.radial
Idomura,et al.,
Comput.
Phys.
Commun.
179,
betweenbetween theradial electric
electric field
field andand thethe temperaturetemperature profile.
profile. In
In thethe simulation
simulation without
without
391source,(2008).parallel velocity
thethethe adaptive
adaptivesource,the parallel
velocity profile
profile peaks
peaks at
at thethe same
same position
position than
than thetheelectric
ititgoes
theedge.
region,turning
electric
field
while
goestotozerozero close
close
theplasma
plasmaedge. In
Inthisthisregion,turning on
onthe
the[5]field
A.whileMishchenko
andtotoA.
J.
Brizard,Phys.
Plasmasadaptivefield
adaptive source
source does
does notnot modify
modify thethe electricfield much.
much. InIn Eq.
Eq. (23),(23), thethe kk coefficient
coefficient isis18,022305(2011). electriccomputed
computedfrom
fromthethesimulation
simulationresults
resultssuch
suchthat
thatthetheFBEFBEisisstrictlystrictlysatisfied.
satisfied. In
Inthetheright
rightS. Jolliet
andisisY.
Idomurasubmitted
to Nucl.
plotof9,
close
toestimate[31]
allplot[6]of figurefigure9, itsits value
value
close
to thethe neoclassicalneoclassicalestimate[31] forforFusion.
all simulations.
simulations.
ItIt isis[7]thenreasonabletosimulations
sourceobey
thenR.
reasonable
to assume
assume thatthat
simulations
with thethe adaptive
adaptivesourceobey aBa forceforce
R. Dominguez
andG. M. withSteabler,
Phys.
Fluidsbalancebalance equationequation TTfbfb ++nnfbfb ++EEfbfb == 0.
0. When
When ρρ∗∗ isis varied,varied, TTfbfb andand nnfbfb areare constantconstant andand5, 3876
(1993).
soso should
should EEfbfb be.
be. InIn otherother words,
words, itit means
means thatthat γγEE should
should be
be proportional
proportional toto ρρ∗∗.. TheThe[8] rates
Y. inCamenen,
Y.flatIdomura,andG.
shearing
shear
aredisplayed
table
1.
shearing
rates
in thethe positive
positive andandflatshear regions
regionsS.
areJolliet,
displayed onontableA.
1. Contrary
Contrary
towithoutadaptivesource,theinto thethe casecase
without thetheadaptive
source,51,
the shearing
shearing rates
rates
in thethe positive
positive shear
shear regions
regions
Peeters,
Nucl.
Fusion073039(2011).

approximatively
approximatively scale
scale like
like ρρ∗∗.. This
This demonstratedemonstrate thatthat thethe γγEEa/ca/css ∝∝ ρρ∗∗ scalingscaling cancan bebebrokenbroken by
by thethe momentum
momentum profile.
profile. This
This broken
broken scalingscaling isis inin fact
fact beneficial
beneficial asas cancan bebeseenseen fromfrom thethe χχi /χ
plot ofof fig.
fig. 8.
8. TheThe shearshear values
values fromfrom table
table 11 show
show thatthat inin thethe flat
flat
i /χGB
GB plot
shear
shear region,region, thethe ρρ∗∗ dependance
dependance isis lessless obvious
obvious between
between thethe ρρ∗∗ == 1/150
1/150 andand ρρ∗∗ == 1/225
1/225
simulations:
simulations: this
this cancan be
be attributed
attributed toto thethe smallsmall value
value ofof γγEE.. ThoseThose values
values also
also reveal
revealthatthat thethe shearshear isis globally
globally smaller
smaller inin thethe scanscan with
with thethe adaptive
adaptive momentum
momentum source,
source,givinggiving aa 20%
20% higher
higher heat
heat transporttransport inin thethe source
source free
free region.
region. This
This suggests
suggests thatthat forfor