GNETコードによるプラズマ加熱解析の現状と今後の発展

村上定義
京都大学大学院工学研究科原子核工学専攻

共同研究者：福山淳、阿久津拓、中島徳嘉、LHD実験グループ
 V. Chan, M. Choi, S.C. Chiu

村上定義、統合コード研究会、2004.3.18

Energetic Particle Confinement in Heliotrons
Confinement of energetic particle is an
important issue for fusion reactor.
α-particle
α-particle heating efficiency (20% ofthermonuclear power)∗ Enegetic particle due to plasma
heatingheating efficiency, broadening ofheating profile
Because of three dimensional magnetic
configuration, the behaviors of trapped
particles in helical ripples are complicatedand would enhance the radial transport inheliotrons.
Thus the confinement of energetic particle
is an important issue for a future reactor
based on the helical system.
∗♦♦

Physical understanding
Quantitative estimation
村上定義、統合コード研究会、2004.3.18

2πpoloidal angle♦

mod-B on flux surface (r/a=0.6)

π00

toroidal angle

2πSimulation Model♦

We solve the drift kenetic equation as a (time-dependent) initialvalue problem based on the Monte Carlo technique.

∂f + (v + v )⋅∇f + v˙ ⋅∇ f −C coll ( f , f ) − Lorbit ( f ) = S( f )//Dv∂t♦Writing the gyrophase averaged distribution function asf (x,v//,v⊥,t)= fbg(r,v2 )+δf (x,v//,v⊥,t)€the linearized drift kinetic equation can be given with initial€condition δf(x,v,t=0)=0
steady state solution (t=∞)∂δf + (v + v )⋅∇δf + v˙ ⋅∇ δf −C coll (δf )− Lorbit (δf ) = S( f )//Dvbg∂t€♦collorbit

C ,Land S are the linear collision operator, orbit loss andthe energy and particle source, respectively.

村上定義、統合コード研究会、2004.3.18

Simulation Model (II)♦It is convenient to introduce the Green function G(x,v,t | x’,v’)
which is defined by the homogeneous F-P equation∂G +(v + v )⋅∇G+ v˙ ⋅∇ G− C(G)− S(G)− L(G)= 0//
Dv∂t♦with the initial condition G(x,v,t = 0|x',v' ) = δ (x − x' )δ (v− v' )Then, the solution of the inhomogeneous problem is given by theconvolution with G;tδf (x,v,t) = ∫ dt' ∫ dx' ∫ dv'S( f bg )G(x,v,t −t'| x',v').
0♦€

In this approach, only the Green function G has to be determinedby the Monte Carlo technique.

村上定義、統合コード研究会、2004.3.18

Monte Carlo Simulation for G♦♦

Magnetic configuration
Finite β effects (magnetic configuration change due to Shafranov shift)
3D MHD equilibrium (VMEC+NEWBOZ)
Complicated particle motion
guiding center motion (µ is conserved.)
Hamiltonian of charged particle

q2 2
H=ρc + µB(ψ ,ϑ ,φ )+qΦ(ψ )2meq. of motion in the Boozer coordinates (ψ,θ,φ,ρc)♦ €Coulomb collisions

Liner Monte Carlo collision operator [Boozer and Kuo-Petravic]
energy and pitch angle scatteringCcoll

1 ∂  2 T ∂δf  ν d ∂2 ∂δf(δf ) = 2 v ν E vδf ++1−λ,()v ∂v m ∂v  2 ∂λ∂λ村上定義、統合コード研究会、2004.3.18

λ=v//vS. Murakami: 01/11/29

ECH Simulation Model♦

The drift kinetic equation for ECH can be written with initialcondition δf(x,v,t=0)=0
steady state solution (t=inf.)∂δf + (v + v )⋅∇δf + v˙ ⋅∇ δf −C coll (δf ) − Lorbit (δf ) = S ql ( f )//Dv
Max∂tcollorbitql♦

C ,Land S are the collision operator, orbit loss and thequasi-linear diffusion operator
for the ECRH, respectively.
Isosurface Plot of Sql♦Isofurface Plot of QL diffusion
term in W7-AS plasma(2nd X-mode ECRH).
n0=2.0x1019m-3, T0=2keV
3D space (r, v//, vperp)
red : positive region
blue : negative region€32vperp
1-2-1v//
村上定義、統合コード研究会、2004.3.18

001

2 0.2

0.1
r/a

Distribution of δf (W7-AS)
power deposition

Isosurface plots of δf

ray-tracing
44V22⊥402

0.2

0 V//r/a-2

0.4
0.6 -4

T=0.06msec
T=0.06msec
村上定義、統合コード研究会、2004.3.18

0GNETEXP

Beam Ion Transport Simulation♦♦

We solve the drift kenetic equation as a (time-dependent) initial value
problem in 5D phase space using the Boozer coordinates(VMEC+NEWBOZ code) .
The drift kinetic equation for a beam source plasma can be given
with the initial condition f(x,v,t=0)=0 [steady state solution (t=inf.)]∂f + (v + v ) ⋅∇f + v˙ ⋅ ∇ f = Ccoll ( f ) + Lparticle ( f ) + S NBI//Dv
∂t♦
♦

coll

C is the collision operator, where we assume a linear collision
operator (energy and pitch agnle scattering).
particle∗ L
is the sink of the distribution; the orbit loss at theoutermost flux surface and the CX loss.
NBI∗ S
is the heat and particle source by the NBI heating.
The beam distribution f is evaluated through a convolution of SNBI
with a characteristic time dependent “Green function”.
In this approach, only the Green function has to be determined by theMonte Carlo technique.
∗村上定義、統合コード研究会、2004.3.18

Trapped Particle Orbit
Rax=3.75m

Rax=3.6m

Rax=3.53m

Toroidalprojection

Classical

E=3.4MeVθp=0.47π
r0=0.5aθ0=π/2φ0=0.0σ-optimized

NC optimized◆

Inward shift of Rax improves the trapped particle orbit.

村上定義、統合コード研究会、2004.3.18

Simulation Results (GNET)
NBI beam ion distributionπ

Rax=3.75mπθpitch Rax=3.53m

θpitchπ/2π/20022406r/a1.0 8

beam injection♦♦v/vthnw
ogdniwslo

406r/av/vth

1.0 8T
Teoeo=T=Tioio=4.0keV,=4.0keV, nne0e0=2.0x10=2.0x101919mm-3-3,, E
Ebb=120keV=120keVThe beam ion distribution by GNET shows increase of the distribution of trappedparticle.
The energy loss rate by orbit loss decreases in the Rax=3.53m config.

村上定義、統合コード研究会、2004.3.18

Comparisons with Simulation Results

104

GNET simulation

410

Rax=3.53m(s37355)
Rax=3.60m(s37350)
Rax=3.75m(s37323)

Rax=3.53m
Rax=3.6m
Rax=3.75m
1000

ne=1.0x1019m-3
1.5-2.0s10010401000

ne=1.0x1019m-3Count [A.U.]♦The count rates are evaluated using a flux averaged beam ion distribution
by GNET simulation.
We can see the similar tendency of beam ion distributions.

Count♦

1005060708090

Energy [keV]

村上定義、統合コード研究会、2004.3.18

100 110 120

10405060

708090

Energy [keV]

100 110 120

Count (70keV<E<80keV) [A.U.]

Exp.

1000b

Count (70keV<E<80keV)/ I [A.U.]

τs Dependence of NDD Count◆◆
◆

100

10
Rax=3.53m
Rax=3.6m
Rax=3.75m
10246
τs [A.U.]

810GNET simulation

1000

100

10
Rax=3.53m
Rax=3.6m
Rax=3.75m
10246τs [A.U.]

8We can see the decrease of the count number as the increase of τs in the
Rax=3.75m case.
The count number of Rax=3.53m is lower than that of Rax=3.6m case.
GNET simulation results show similar tendency with experimental results.
S. Murakami et al., Joint Meeting of US-Japan Workshop and Kyoto Univ. 21st COE Symposium, 3 March 200410τs Dependence of Energy Loss Rate

Energy loss rate (orbit) [%]◆Rax=3.53m
Rax=3.6m
Rax=3.75m

20151050◆

Energy loss rate (orbit + CX) [%]GNET simulation

250246τs [A.U.]

810GNET simulation

252015
10Rax=3.53m
Rax=3.6m
Rax=3.75m

500246
τs [A.U.]

8Energy loss rate is higher in the Rax=3.75m case.
No clear difference of the total loss rates between the Rax=3.53 and 3.6m.
=> Because of low electron temperature.
S. Murakami et al., Joint Meeting of US-Japan Workshop and Kyoto Univ. 21st COE Symposium, 3 March 2004

10Radial Distribution of Beam Ionsf(70keV<E<80keV; v//=0)

10

5 10

Rax=3.53m
Rax=3.60m
Rax=3.75m

300f(v//=0) nf(v //=0)

200
150

2 105
3 10102 1010

00◆◆
◆

0.2

0.4

r/a

0.60.81

1.5 105
1 105

10050

Rax=3.53m
Rax=3.60m
Rax=3.75m

2.5 105

4 1010n250

total f(70keV<E<80keV)

3 105f(total)

350

Rax=3.53m
Rax=3.6m
Rax=3.75m

1 1010

5 104

00
0

0.2

0.4

0.6
r/a0.8

1

1.2

00.2

0.4r/a0.6

0.8

1

Radial profile of trapped ions (70keV < E < 80keV) shows a good confinementin the Rax=3.53 configuration.
A larger distribution near the edge in the Rax= 3.6m configuration.
Neutral weighted distribution indicates a higher count rate in the Rax=3.6m.

S. Murakami et al., Joint Meeting of US-Japan Workshop and Kyoto Univ. 21st COE Symposium, 3 March 2004

GNET Simulation Model♦

We solve the drift kinetic equation as a (time-dependent) initialvalue problem in 5D phase space based on the Monte Carlotechnique.

∂f min+(v// +vD )⋅∇f min +a⋅∇v f min −C( f min )−Q ( f min )−L =S
∂t
ICRF♦

C(f) : linear Clulomb Collision Operator€QICRF : ICRF heating term
wave-particle interaction modelSparticle : particle source=> by ionization of neutral particle(AURORA code)
Lparticle : particle sink (loss)=> Charge exchange loss=> Orbit loss (outermost flux surface)
The minority ion distribution f is evaluatedthrough a convolution of Sparticle with a characteristic
time dependent Green function.
村上定義、統合コード研究会、2004.3.18

particleparticleminority ionsource
n0=2.0x1019m-3
n0=1.0x1019m-3

Energetic Ion Distribution by ICH (LHD)asco psingtedr.
c n t singspakeV[ygreEn]

500prap

pminority = 2π ∫400300200100πPitch π/2angle00.5r/a0 1.0LHD(Rax=3.6m)
19m-3-3 B=2.75T@R=3.6m,
T=T io=1.6keV,=1.0x1019
Teo=1.6keV, nne0
eo=Tioe0=1.0x10 m B=2.75T@R=3.6m,
ffRF=38.47MHz, k =5m-1-1, k perp=62.8m~2.5MW=62.8m-1-1,, PPICH
RF=38.47MHz, k////=5m , kperp
ICH~2.5MW村上定義、統合コード研究会、2004.3.18

12
mv f (v //,v⊥ )dv //v ⊥ dv ⊥2

Simulation, Modeling and Validation of ICRH Physics
ICHシミュレーションコード
GNET: S. Murakami
ORBIT-RF: GA理論グループ♦ 共通の部分∗ Monte-Carlo codes following fast ion drift orbits in 2- and 3-D
magnetic confinement topologies, under going Coulomb collisions andICRF quasilinear heating modeled by Monte-Carlo techniques.
∗ Both codes can be used to study fast ion generation by minority ICRF
heating, self-consistent inclusion of neoclassical transport,enhancement of transport by 3-D magnetic fields, plasma rotation andmomentum transport, and interaction with NBI fast ions.
♦ 異なる部分：∗ GNET uses a particle source term to achieve steady-state∗ ORBIT-RF uses particle re-injection.
∗ GNET includes charge exchange losses and its collision operator
includes energy scattering. ORBIT-RF is implementing higher
harmonic ICRH.
村上定義、統合コード研究会、2004.3.18

Task 1♦♦♦Since the two codes are have been developed independently, it is useful to
do detailed benchmarking of the two codes before proceeding with jointly
development of new capabilities.
Task 1. Benchmarking between GNET and ORBIT-RF
Timescale: 6 months∗ Simulate minority ICRH at the fundamental frequency for aselected DIII-D discharge∗ Simulate NBI and evaluate NBCD for a selected DIII-D discharge∗ Simulate minority ICRH and NBI interaction for a selected DIII-Ddischarge
In all cases, we will compare absorption power (location and amplitude),fast ion spectrum, particle flux, and detailed convergence studies. S.
Murakami will implement an interface into GNET to use
magnetic/equilibrium data files from EFIT as currently done in ORBITRF. M. Choi (GA) will provide DIII-D kinetic profiles andmagnetic/equilibrium data files to S. Murakami.

村上定義、統合コード研究会、2004.3.18

Task 1♦Longer-term Topics:∗ Both GNET and ORBIT-RF have plans to improve the RF wave-fieldcalculation by coupling to a full wave-code. Benchmarking will bedone following the implementation.
∗ Identify momentum transport and plasma rotation problem for joint
study∗ ORBIT-RF plans to implement multiple N-parallel wave capability
and CX losses for comparison with GNET

村上定義、統合コード研究会、2004.3.18

Simulation, Modeling and Validation of ECRH Physics♦♦♦♦High power ECRH and ECCD for steady-state operation, profile control
and tokamak performance improvement is a major element of the DIII-D
experimental program.
The research is highly relevant to ITER advanced operation. ECRH is also
important for stellarator both for heating and controlling transport.
The GNET code operated in the quasilinear mode can also be used for
ECRH study. At GA, a state-of-the-art ray-tracing code TORAY-GA andthe CQL3D Fokker-Planck code are being used for modeling and planning
experiments.
GA is interested in adding the GNET code to the suite of codes for study
ECRH and ECCD. The following collaboration is proposed.

村上定義、統合コード研究会、2004.3.18

Task 2♦♦

Task 2. Preparation of GNET code for study of DIII-D ECRH
experiments.
Timescale: 6 months∗ S. Murakami and S.C. Chiu (GA) will formulate and implement
momentum conserving e-e collision operator in GNET for CDcalculation∗ S. Murakami will add Ohmic heating capability to GNET
∗ S. Murakami and V. Chan will complete calculation of self-consistentelectrostatic potential generation by ECRH
∗ GNET will be installed at GA after completion of the first threesubtasks, and M. Choi will do comparative calculations against
CQL3D∗ S. Murakami and GA will jointly investigate the effect of error field
on fast electron transport using GNET
Longer-term Topics:∗ Study effects of magnetic fluctuations and magnetic islands onelectron transport and bootstrap current∗ Extend capability to EBW and IBW

村上定義、統合コード研究会、2004.3.18

Simulation Model♦

We solve the drift kenetic equation as a (time-dependent) initialvalue problem based on the Monte Carlo technique.

∂f + (v + v )⋅∇f + v˙ ⋅∇ f −C coll ( f , f ) − Lorbit ( f ) = S( f )
//
Dv
∂t
♦

€

Writing the gyrophase averaged distribution function asf (x,v//,v⊥,t)= fbg(r,v2 )+δf (x,v//,v⊥,t)
the linearized drift kinetic equation can be given with initialcondition
δf(x,v,t=0)=0
steady state solution (t=∞)€∂δf + (v + v )⋅∇δf + v˙ ⋅∇ δf −C coll (δf , f ) − S(δf ) − Lorbit (δf )//Dvbg∂t= S( f bg )+ S neo ( f bg )+C coll ( f bg,δf )♦€collorbit

C ,Land S are the linear collision operator, orbit loss andthe energy and particle source, respectively.
Sneo is the usual driving term for neoclassical transport.
Sneo

∂f bg
∂f bg
˙
= −(VD ) r−v∂r∂v

村上定義、統合コード研究会、2004.3.18

Energy and Momentum Conservation♦♦

Conservation of Energy and momentum is important inconsidering the current drive and the neoclassical transport.
Energy and momentum conservation (locally) can be possible bysolving the drift kinetic equation iteratively. δf = δf0+δf1 +δf2+…
for l = 0 (s-species)(v // + v D ) ⋅∇δf0s + v˙ ⋅ ∇vδf0s − C coll (δf0s ) − Sql (δf 0s ) − Lorbit (δf0s)s ) + S neo ( f s )= S ql( f Max
Max

for l = n(electron)(v // + v D ) ⋅∇δfne + v˙ ⋅∇ vδfne − Ccoll (δf ne ) − Sql (δf ne ) − Lorbit (δf ne )e , δf e ) + Ccoll ( f e , δf i )= Ccoll ( f Maxn −1Maxn −1

for l = n(ion)(v // + v D ) ⋅∇δfni + v˙ ⋅∇ vδfni − Ccoll (δf ni ) − Sql (δfni ) − Lorbit(δfni )i , δf i ) + Ccoll ( f i , δf e )= Ccoll ( f Maxn −1
Maxn −1
村上定義、統合コード研究会、2004.3.18

Global Transport Simulation by GNET
We have developed a GNET(Global NEoclassical Transport) code solving
drift kinetic equation in 5D phase-space.
♦ We study the energetic particle transport in non-axisymmetricconfigurations by GNET code.
♦ ECRH generated suprathermal electron transport
W7-AS, CHS, LHD(collaboration with Max-Planck IPP)
♦ NBI generated beam ion transport
LHD, CHS♦ ICH generated energetic tail ion transport
LHD, W7-X♦ Future Plan (collaboration with GA group)♦ Benchmark with ORBIT-RF♦ Momentum and energy conserving collision operator♦ Investigation of 3D problem in tokamaks♦

村上定義、統合コード研究会、2004.3.18

S. Murakami: 01/11/29