Evolution of Nonthermal Particle Distributions in Radio Frequency Heating of Fusion Plasmas

1
Evolution of Nonthermal Particle Distributions in
Radio Frequency Heating of Fusion Plasmas

• Paul Bonoli on behalf of the SciDAC Center for
  Simulation of Wave-Plasma Interactions

SciDAC 2007 Conference Boston, Massachusetts June
25-28, 2007
2
Participants in the Center for Simulation of
Wave Plasma Interactions
L.A. Berry, D.B. Batchelor, E.F. Jaeger, E.
DAzevedo, M. Carter
D. Smithe
C.K. Phillips, E. Valeo N. Gorelenkov, H. Qin
P.T. Bonoli, J.C. Wright H. Kohno
R.W. Harvey, A.P. Smirnov N.M. Ershov
M. Brambilla R. Bilato
Politecnico di Torino R. Maggiora V. Lancellotti
M. Choi
D. DIppolito, J. Myra - Lodestar Research
3
Background and Motivation for Work
4
The role of RF power in fusion plasmas has
evolved over the years

• Early applications involved bulk plasma heating
  and non-inductive maintenance of the entire
  plasma current
• Ion cyclotron resonance heating (ICRH), lower
  hybrid heating (LHH), and electron cyclotron
  resonance heating (ECRH).
• Lower hybrid current drive (LHCD) was especially
  successful and efficient.
• More recent applications have utilized RF waves
  for localized control of the plasma current and
  pressure profiles
• LH and EC current drive for sawtooth control, NTM
  control.
• Mode converted ICRF waves for shear flow
  generation and current generation.
• Current profile control for access to high
  confinement (advanced tokamak) regimes with
  high bootstrap current fraction (gt 70).
• Most of the RF applications listed above are
  planned for ITER
• A predictive capability is needed to insure
  success in ITER applications.

5
So what are some of the challenges ?

• Complete description of the wave-particle
  interaction involves integrating several
  separate calculations
• Antenna coupling (linear and nonlinear response)
• Wave propagation
• Wave absorption
• Multiple spatial scales can occur
• ICRF Mode conversion and current drive
• RF sheath formation in the edge
• RF waves can interact in a nonlinear fashion with
  thermal and nonthermal particles in the plasma
  (long timescale)
• Nonthermal electron and ion distributions
  produced by the RF wave itself.
• Nonthermal ions already present in the plasma due
  to neutral beam injection (NBI) and fusion
  reactions (alpha particles).

6
ICRF Heating Involves Two Important
ProcessesIon Cyclotron Absorption with
Nonthermal Ion Tail Production Conversion of
Fast Wave to Short Wavelength Modes
7
Physics of Nonthermal Ion Tail Evolution

• If a minority ion species (5 hydrogen) is
  present in a majority ion species plasma (95
  deuterium) then an RF wave at the cyclotron
  frequency of the minority ion will have an
  electric field component with the same
  polarization as the minority ion
• Secular interaction - wave damps its power via
  cyclotron absorption on the minority ion.
• Nonthermal, anisotropic minority ion tail is
  generated that slows down and heats background
  electrons (drag) and majority ions (collisions).
• But this process is greatly complicated by a
  phenomenon known as mode conversion.

8
Dissection of the minority heating scheme in a
controlled fusion device the tokamakSimulation
and detection of mode converted ion cyclotron RF
waves
9
Mode Conversion Changes Polarization of Incoming
Fast Wave, Thus Modifying Resonant Ion Cyclotron
Absorption
10
Two full-wave solvers have been advanced within
the RF SciDac Center

• All Orders Spectral Algorithm (AORSA) 1D, 2D
  3D (Jaeger)
• Spectral in all 3 dimensions
• Cylindrical coordinates (x, y, ?)
• Includes all cyclotron harmonics
• No approximation of small particle gyro radius r
  compared to wavelength l
• Produces huge, dense, non-symmetric, indefinite,
  complex matrices
• TORIC 2D (Brambilla/Bonoli/Wright)
• Mixed spectral (toroidal, poloidal), finite
  element (radial)
• Flux coordinates
• Up 2nd cyclotron harmonic
• Expanded to 2nd order in r/l
• Sparse banded matrices, with dense blocks

Blowup region
Slow ion cyclotron wave
Electrostatic ion Bernstein wave
11
Understanding the mode conversion aspect of ICRF
heating required theory, experiment, and
computation

• Initial observations of mode converted ICRF waves
  in Alcator C-Mod presented a scientific
  conundrum
• Waves were detected on the tokamak LFS and at
  kR7 cm-1
• This was the wrong location and wavenumber to
  be the anticipated ion Bernstein wave (IBW)
• But our full-wave simulations (TORIC and AORSA)
  also revealed the presence of these waves at the
  wrong location and wavenumber.

12
Experimental Observation of a new type of mode
conversion- the ion cyclotron wave (ICW)
A. Mazurenko, PhD thesis, MIT (2001). Nelson-Melby
et al, PRL 90 (2003) 155004
13
Both the TORIC and AORSA Solvers also predicted
the ICW wave field feature
ICW
IBW
FW

• TORIC at 240Nr x 255 Nm AORSA at 230Nx
  x 230 Ny
• In fact the mode converted ICW had been predicted
  to exist years ago F. W. Perkins, Nuclear Fusion
  (1977) but had been forgotten.

14
Predicted RF Electric Field from the TORIC field
solver has been used in a synthetic diagnostic
code for the PCI
TORIC E field is used in synthetic PCI on C-Mod
that measures the perturbed density due to mode
converted ICRF waves.
Y. Lin, A. Parisot, J. Wright PoP, 2005, PPCF,
2005

• Shear flow generation for pressure profile
  control is possible in fast wave to ion cyclotron
  wave mode conversion (Myra, Jaeger et al Phys Rev
  Lett, 2003 and Melby et al Phys Rev Lett. 2003)

15
Calculations on the Cray XT3 Jaguar have allowed
the first simulations of mode conversion in ITER
ITER with DTHE3 202030 with NR NZ 350,
f 53 MHz, n 2.5x1019 m-3 (4096 processors for
1.5 hours on the Cray XT-3)
Blowup (E_parallel)
E_perp
Mode converted Ion Cyclotron Wave (ICW)
Calculation of flow or current drive for full
antenna spectrum in ITER requires petaflop
capability
16
Scaling of Full-wave ICRF solvers to gt 20,000
processors demonstrated for ICW Mode Conversion
in ITER in preparation for ITER mode conversion
studies
ITER with DTHE3 202030 with NR NZ 500,
f 53 MHz, n 2.5x1019 m-3
17
Dissection of the minority heating scheme in a
controlled fusion device the tokamakSimulation
and measurements of velocity space structure of
nonthermal ion distributions.
18
Wave propagation and the plasma response are
governed by the Maxwell-Boltzmann system of
equations
For time harmonic (rapidly oscillating) wave
fields E with frequency ?, Maxwells equations
reduce to the Helmholtz wave equation
Wave Solvers (AORSA) (TORIC)
The plasma current (Jp) is a non-local, integral
operator (and non-linear) on the rf electric
field and conductivity kernel
SIGMAD Module gives ? (f0,s)
The long time scale response of the plasma
distribution function is obtained from the bounce
averaged Fokker-Planck equation
Plasma Response (CQL3D)
0
where
Need to solve this nonlinear, integral set of
equations for wave fields and velocity
distribution function self-consistently. This
requires an iterative process to attain
self-consistency.
19
Calculation for C-Mod minority H, NR 128, NZ
128,256 processors for 3 hrs on Cray XT3 ORNL
20
2D Field and Dissipation Contours Show that
Heating is Concentrated at Ion Turning Points on
the Minority Resonance Chord at r/a ? 0.45
Heating (H)
Wave fields
21
Experimental measurements of the energetic ion
tail on C-Mod have been made using a compact
neutral particle analyzer
Good agreement between simulated and measured
tail temperature Tion ? 70 keV
Courtesy of V. Tang, PhD Thesis, MIT (2006)
22
3D (r, V? , V//) distribution function from CQL3D
AORSA reproduces CNPA measurements using a
synthetic code diagnostic
Building this synthetic diagnostic required a
close collaboration between theory and experiment
(V. Tang and R. Harvey)
Courtesy of V. Tang, PhD Thesis, MIT (2006) also
PPCF, 49, 873 (2007).
23
Outstanding challenges for the near
futureFinite ion drift orbit effects
24
Simulations of high harmonic fast wave (HHFW)
fast ion beam interaction in DIII-D are still
unresolved
DIIID high density L-mode
Stronger Beam Interactions at 4?D (60 MHz) Than
at 8?D (116 MHz) Observed in DIII-D CQL3D-AORSA
predicts increased absorption as frequency was
raised in disagreement with expt. Monte Carlo
ORBIT code (ORBIT-RF) combined with an RF
operator (using fields from TORIC solver) does
reproduce the experimental trend.
?1014 (/s)
Neutron reaction rate
Power
Sn neutron enhancement factor
217-05/MC/jy
25
Orbit topology modifies wave-particle resonance

• Shown at right are trajectories for 12 particles
  in the C-Mod case
• 4 equi-spaced velocities
• 3 equi-spaced ? velocities
• 409,600 complete poloidal orbits
• Particle cyclotron resonances and strong
  quasilinear diffusion occur in roughly vertical
  planes in zero-orbit width description.
• But orbit topology can move particles away from
  (or towards) resonances that would be sampled
  (not sampled) in full-wave solver.

26
We are investigating finite ion drift orbit
effects using two approaches

• The diffusion coefficient (D) has been evaluated
  by a direct orbit integration using electric
  fields from AORSA
• The DC code computes averages of the changes in
  velocity, pitch angle, and radial position over a
  complete bounce orbit, to obtain a set of RF
  induced diffusion coefficients.
• Diffusion Coefficient calculations done on CRAY
  XT3 (ORNL) using 256 processors _at_ 10 min.
• The Monte Carlo code ORBIT RF has been combined
  with the TORIC ICRF solver
• Self-consistent iteration not yet carried out.
• ORBIT RF code ported to JAGUAR where good scaling
  to gt 1000 processor cores has been demonstrated.
• We are now examining best way to pass statistical
  distribution from ORBIT RF to TORIC AORSA to do
  self-consistent iteration.

27
Monte Carlo ORBIT Code has been coupled to the
TORIC full-wave solver through an RF
Operator Finite orbit effects can be studied
quantitatively using this approach QL Diffusion
Operator Formulated in terms of Multi-Fourier
Poloidal Modes from the TORIC ICRF Solver and
used to compute increment in magnetic moment due
to the ICRF interaction
Numerical distribution function from ORBIT-RF not
yet coupled back to TORIC or AORSA this is an
important and difficult next step !
28
Beam pressure computed with f(E) from ORBIT-RF
agrees qualitatively with experiment but
iteration is still needed between the ORBIT code
and full-wave solver
Particle distribution f(E)
ORBIT-RF
?104
4th harmonic
NB only
8th harmonic
116MHz (1.7MW)
Beam injection energy (80keV)
Sn(ORBIT-RF)1.2
60 MHz (0.8MW)
Sn(ORBIT-RF)1.9
Beam Ion Pressure (N/m 2 )
r/a
Energy (keV)
217-05/MC/jy
29
Outstanding challenges for the near
futureNonlinear effects at the RF
antenna-tenuous edge plasma
30
ICRF Launchers in Contact with Plasma are Subject
to Nonlinear Effects, Leading to Parasitic Power
Losses.
RF sheaths can form due to mismatch between
equilibrium B and the antenna structure (E//),
resulting in power dissipation
Alcator C-Mod Dipole Antenna
Electrons are preferentially accelerated out of
the sheath region. A DC voltage Vrf is set up
to maintain ambipolarity.
31
Capability to efficiently compute 3D wave fields
will be important for assessing antenna edge
interaction, especially in weak single pass
damping regime
NSTX ITER
Scenario 2
NSTX simulation summed over 81 toroidal
modes.ITER simulation summed over 169 toroidal
modes)AORSA run on JAGUAR using 2048 processors
for 8 hrs
32
Two Approaches are Being Pursued to Study the
Nonlinear ICRF antenna edge Interaction

• Implementation of RF sheath boundary conditions
  in full-wave solver (spectral solution)
• Start with linear field response from a coupled
  full-wave field solver (TORIC) and 3D
  electromagnetic antenna code (TOPICA).
• Modify metal wall BC in field solver to include
  sheath dissipation and then iterate with antenna
  code.
• Approach will quantify how much ICRF power is
  coupled to the plasma.
• Time domain simulations using 3D EM field solver
  - VORPAL
• Fully implicit time domain dielectric response
  module has been implemented for electrons and
  ions.
• Use PIC-treatment in future for ion response
  (fully nonlinear).

33
VORPAL Time Domain Simulation of Antenna LoopD.
Smithe, Poster at this Conference (Tuesday
Evening)
34
VORPAL Time Domain Simulation of Antenna Loop
Surface E-field on Loop antenna
RF B-field
Surface E-field
Wavefronts
Radiation Pattern
Radiation from Behind
35
Summary

• Combined full-wave and Fokker Planck solvers
  (CQL3D-SIGMAD-AORSA) have been used to simulate
  minority ICRF heating in present day tokamaks and
  in ITER
• Full coupling uses self-consistent nonthermal ion
  distributions and quasilinear diffusion
  coefficient in differential form.
• Comparison with synthetic diagnostic (CNPA) and
  experiments validates the use of this simulation
  to predict these experiments in burning plasmas
  such as ITER.
• Both AORSA and TORIC can simulate ICW/IBW mode
  conversion in present day tokamaks and in ITER
• Comparison of synthetic diagnostic (PCI) with
  experiment validates predictive capability of
  simulation
• Can now assess use of mode converted waves for
  pressure and current profile control in present
  day tokamaks and for ITER.

36
Summary

• Fully self-consistent coupling of our full-wave
  solvers to a Monte Carlo orbit code (ORBIT RF) is
  underway
• Full coupling will use statistical nonthermal ion
  distributions and quasilinear diffusion
  coefficient in differential form
• A self-consistent treatment of the RF antenna
  edge plasma is underway
• Linear antenna coupling problem is substantially
  completed using the TORIC TOPICA suite.
• Boundary conditions for sheaths are being
  implemented in our full-wave solver (TORIC)
• Proof of principle time domain simulations of
  sheath formation have been done using VORPAL and
  will now be extended to the nonlinear regime
  using PIC treatment for ions.