Simulations of particle dynamics in a global toroidal geometry
The gyrokinetic toroidal code (GTC) has been upgraded for global simulations by coupling the core and scrape-off layer (SOL) regions across the separatrix with field-aligned particle-grid interpolations. A fully kinetic particle pusher for high frequency waves (ion cyclotron frequency and beyond) and a guiding center pusher for low frequency waves have been implemented using cylindrical coordinates in a global toroidal geometry. The two integrators correctly capture the particle orbits and agree well with each other, conserving energy and canonical angular momentum. As a verification and application of this new capability, ion guiding center simulations have been carried out to study ion orbit losses at the edge of the DIII-D tokamak for single null magnetic separatrix discharges. The ion loss conditions are examined as a function of the pitch angle for cases without and with an electric field. The simulations show good agreement with past theoretical results and with many experimentally observed features. A measure of the ion direct orbit loss fraction shows that the loss fraction increases with the ion energy for DIII-D in the initial velocity space.
One of the important challenges in getting to a viable operating regime for ITER and future fusion reactors is associated with the nonlinear turbulent dynamics of the plasma in the scrape-off layer (SOL) loarte2007power (). The plasma characteristics in SOL can greatly affect the overall confinement properties of the device and also regulate the heat load to the tokamak wall. It can also influence the level of fusion ash, the impurity dynamics, sheath physics and plasma shaping effects. Furthermore, the SOL dynamics can degrade the current drive performance of radio frequency (RF) waves through its impact on the density threshold conditions for the onset of parametric decay instabilities cesario2010current (). An indepth undertanding of the mechanisms determining the width of the SOL layer remains an outstanding open problem. A study of the SOL plasma dynamics is challenging due to multiple spatial and temporal scales associated with different energy sources (instabilities) in that region. Fluid simulation transport codes such as UEDGE wising1996simulation (), SOLPS rozhansky2009new () are normally used to simulate the SOL dynamics. These fluid codes use a set of fluid transport equations that are based on the Braginskii equations braginskii1965review (). However, the results show a number of discrepancies between experimental findings and fluid simulations including in the characteristics of the radial electric field, parallel ion flow, impurity radiation, etc. chankin2007discrepancy (); erents2004comparison (); wischmeier2011assessment (). It is believed that kinetic effects could be a significant contributor in the SOL to processes like ion orbit losses stacey2013effect (), X point losses chang2002x (), nonlocal turbulent transport omotani2013non (), plasma sheath dynamicswesson1995effect (), parametric decay instabilities cesario2010current (); kuley2009stabilization (); kuley2010parametric (); kuley2010lower (), etc. To correctly model many of these effects requires a kinetic approach that covers the closed, and open field line regions across the separatrix, and includes realistic SOL physics and tokamak geometry. Due to the difficulty of accessibility of diagnostics in the SOL region such global kinetic simulations can help in developing useful insights for predicting the plasma dynamics in that region for present and future reactors such as ITER and DEMO. A laudable effort in this directions has been the development of the massively parallel kinetic simulation code XGC-1 chang2009compressed () that takes a first-principles approach and has emerged as an efficient method for describing the complex physics of turbulent transport. Another widely used and successful tool, the gyrokinetic toroidal code (GTC) has also undergone continuous development for the last two decades and has been applied in the study of plasma transport in the core region lin1998turbulent (). GTC is a well-benchmarked, first principles code which has been extensively applied to study neoclassical transport lin1997neoclassical (); dong2016effects (), microturbulence xie2017new (); fulton2016gyrokinetic (); xiao2009turbulent (); holod2016effect (); xiao2015gyrokinetic (); zhang2012nonlinear (), mesoscale Alfvn eigenmodes zhang2012nonlinear (); spong2012verification (); wang2013radial () excited by energetic particles, macroscopic MHD modes mcclenaghan2014verification (); liu2014verification (); bao2017conservative () (kink and tearing modes) and radio frequency (RF) waves kuley2013verification (); bao2014particle (); kuley2015verification (); kuley2015nonlinear (); bao2015global (); bao2016electromagnetic (); bao2016nonlinear () in the core region. However, the assumptions used in studying turbulence in the core region may not be valid in the SOL region. GTC normally uses conventional magnetic flux coordinates, in which the equations of motion encounter a mathematical singularity of metric on the magnetic separatrix surface. Recently, the GTC was extended to separately study instabilities in the SOL and core regions of a field reversed configuration (FRC) using Boozer coordinates. However, the code still did not have the capability to couple these two regionsfulton2016gyrokinetic (); fulton2016bgyrokinetic (). The difficulty lay in the discontinuity of the poloidal angle in Boozer coordinates across the separatrix. This limitation restricted the code’s usage to electrostatic simulations either in the core or the SOL region with no cross-separatrix coupling.
In our present work we report a significant enhancement of the GTC through development of a new global nonlinear particle simulation model that couples the core and SOL regions. The code also provides a realistic treatment of the separatrix region through Equilibrium Fitting (EFIT) lao1990equilibrium (); ren2011high () of equilibrium data files
generated by experimental discharges. A particular feature of the upgraded GTC is the use of a cylindrical coordinate system for the advancement of the particle dynamics, which allows particle motion in arbitrary shaped flux surfaces including the magnetic separatrix and the magnetic X-point.
As a first step in developing this nonlinear particle simulation model, we have constructed a field aligned computational mesh to take advantage of the smallest number of grid points in the direction of the magnetic field with high resolution in any given poloidal plane using more fundamental cylindrical coordinates rather than the flux coordinates. Field aligned coordinate approaches have been used in the past in turbulence simulation codes such as GTC xiao2015gyrokinetic (), XGC1 kim2017happens (), FENICA hariri2013flux (), GEM Sturdevant2017 (), etc. The gain in computational efficiency by using appropriate coordinates and computational mesh help to optimize turbulence simulations of large devices like ITER and DEMO. We have also developed a fully kinetic (FK) particle pusher to capture the effect of high frequency waves (ion cyclotron frequency and beyond) and
a guiding center (GC) pusher to describe the particle dynamics associated with low frequency waves (much smaller than the ion cyclotron frequency).
To test the effectiveness of these enhancements and appropriately benchmark the code we have carried out ion guiding center simulations to study ion orbit losses at the edge of the DIII-D tokamak for single null magnetic separatrix discharges. Using model calculations some analytic expressions of such losses have been presented by Miyamoto miyamoto1996direct () and have been used in the past to estimate the loss region in velocity space for JET, JT-60 and ITER magnetic geometries. Stacey has introduced the effect of ion orbit loss and X point loss in his fluid calculations for the interpretation of fluid transport in the edge region stacey2013effect (); stacey2011effect (). Our extended simulation model is meant to benchmark GTC against some of these past results. We have examined the ion orbit losses as a function of the pitch angle both in the presence and absence of an electric field. The simulations show good agreement with past theoretical results and with many experimentally observed features. A measure of the ion direct orbit loss fraction shows that the loss fraction increases with the ion energy for DIII-D in the initial velocity space.
The paper is organized as follows: The next section II contains a detailed description of the coordinate system and the representation of the equilibrium and fluctuating field quantities in the code. Section III describes the computational mesh used in the simulations. The physics model explained the dynamics of the particles is described in section IV followed by our simulation results for the ion orbit losses in section V. Section VI provides a brief summary and some concluding discussions.
Ii Representation of equilibrium and coordinate system
In describing the collective dynamics of a tokamak plasma, it is convenient to represent the physical variables in terms of equilibrium and fluctuating quantities. The equilibrium is described by the Grad-Shafranov equation, while the fluctuations representing wave activity contribute to plasma transport. The input equilibrium magnetic field configuration for GTC can be generated in any of the equilibrium solvers EFIT, VMEC, or M3D-C1 and is normally expressed as a function of the magnetic flux function. In the MHD approximation (ignoring guiding center orbit effects) the plasma profiles of temperature, density, current etc. are also functions of the flux function. In the upgraded GTC version we use a cylindrical set of coordinates where is the distance from the geometry axis, is the toroidal angle and is the direction of the torus symmetry axis. The representation of the magnetic field for an axisymmetric system is then given by,
where is the poloidal flux function that labels the magnetic surfaces for both closed and open field lines [cf. Fig 1(b)] and is the poloidal current function [cf. Fig 1(a)], which provides the components of magnetic field in cylindrical coordinates as:
The Jacobian for this system can be written as
The cylindrical toroidal coordinate system is related to the standard Cartesian system as follows
By defining contravariant basis vectors , , the velocity and the electric field can be written as
The equilibrium inputs from EFIT only provide equilibrium quantities on a coarse mesh, which usually contains a few tens of grid points in the and directions. However, the micro-scale turbulence demands much denser grid points in and directions. Therefore, it is necessary to map the equilibrium mesh to a dense computational mesh to achieve sufficient numerical accuracy. For a 1D function , such as the poloidal current function we can use the following B-spline representation xiao2015gyrokinetic ()
where and are coefficients related to the first and second order differential in direction, which are calculated using finite difference method on a spline mesh. In our simulations we have calculated the grid size in using the following three steps:
where is the poloidal flux function at the limiter point, is the number of grid points in , is the grid point number at the separatix, and is the poloidal flux function at the separatrix. This calculation will provide accurate value of at the separatrix in the simulation, however will give some minor difference in calculating the at the limiter points.
A 2D function can be expressed as
where , , and , . Eq. (LABEL:eq:2Dspline) is derived by using the 1D B-spline functions of and . The spline coefficients are calculated from the spline coefficients and . Fig. 1(a) and Fig. 1(b) represent the equilibrium poloidal current function on uniform flux grid and poloidal flux function on rectangular R-Z grid for DIII-D shot 158103 at 3050 ms, respectively PhysRevLett.114.105002 (). In GTC we use these two functions in Eq.(2) to calculate the magnetic field components for DIII-D [cf. Fig. 2]. In particle pusher we interpolate these field quantities at the particle position using 2D spline function, as described in Eq. (LABEL:eq:2Dspline).
Iii Field aligned mesh for fluctuating quantities
The field aligned mesh has maximum numerical efficiency and accuracy to address the nonlinear physics of drift wave turbulence. Because the particle moves much faster in the direction parallel to the magnetic field than the drifts across the magnetic field, the parallel wave vector is usually much smaller than the perpendicular wave vector for the drift wave instabilities. Thus, it only requires a small number of grid points to resolve the parallel wave length, which greatly saves the computational costs and suppresses the numerical high modes efficiently. Furthermore, it helps to simplify the implementation of the field solver, since the field aligned mesh can exactly decouple the parallel and perpendicular directions.
iii.1 Radial grid:
GTC code has in the past utilized the field aligned mesh in Boozer coordinates and successfully applied it for gyrokinetic simulations of drift waves in the core region of toroidal plasmas lin1998turbulent (). In this work, we extend the original GTC field aligned mesh from Boozer coordinates in the core region to cylindrical coordinates in the whole domain across the separatrix. In order to create the field aligned mesh for the whole tokamak domain, we apply an equally spaced radial grid (spacing of size ) on the outer mid-plane for each field line, which can be calculated as:
where is the radial position as a function of poloidal flux and on the low field side, and is the number of field line in the SOL region. is the maximum value of the poloidal flux on the limiter points (plasma facing components). GTC also has the capability to use a nonuniform grid with grid size in the perpendicular direction correlated with the local gyro-radius. The poloidal flux of each field line for the simulation mesh in the SOL region is:
Then the field line number in the core region can be calculated as:
where is the major radius. The poloidal flux of the innermost grids in the core region can be determined as:
where is not necessarily equal to the poloidal flux value at the magnetic axis. The poloidal flux of each field line for the simulation mesh in the core region is:
Similarly, the field line number in the private region is calculated as:
where is the poloidal flux of the innermost flux surface in the private region, and is the minimum value in direction. The poloidal flux of the innermost grids in the private region can be determined as:
where is not necessarily equal to . The poloidal flux of each field line for the simulation mesh in the private region is:
iii.2 Poloidal grid:
Next we trace each field line along the poloidal direction and calculate the length of the poloidal projection for each field line by using the unit magnetic vector as:
where , and is the poloidal magnetic field. The starting points for tracing the field lines are on the low field side of out-mid plane with for the core region and on the low field side of plane for the SOL and private regions. The step for tracing each field line is at least smaller than the order of magnitude of poloidal projected field line length, and we correct the new position for each tracing step to ensure that it lies on the same flux surface by using Newtonâs method with the error of the order of floating precision. All the advanced positions during tracing poloidal field line are recorded as tracing grids. The outmost simulation boundary is determined by the limiter as shown by the brown line, so the tracing grids outside the limiter in the SOL and private region are removed by using ray casting algorithm wiki.org (). In the closed field line region, we calculate the length of the poloidal projection of the field line for each flux surface. However, in the open field line regions, we set the intersection points formed by field lines and limiter as the simulation boundary points and calculate the length of the poloidal projection of the field line between each two intersection points on one single continuous field line.
The simulation grid number for each continuous poloidal field line is calculated as:
where and are the poloidal length and the poloidal grid number for i-th field line, respectively, and is the given approximate grid size along the poloidal direction. The exact poloidal grid size is then determined by and as:
The position of the field aligned grid along poloidal direction can be derived from poloidal grid size , the much smaller tracing step and the positions of the dense tracing grids.
In the cylindrical coordinate representation, we only apply the field aligned mesh in the poloidal direction, and the meshes on different poloidal plane are identical as shown by Fig. 3 (axisymmetric). In order to decouple the exact parallel and perpendicular directions and use a small number of poloidal planes, the particle-grid gathering and scattering are along the exact parallel direction by interpolating the value on the identical field line.
Iv Physics model for particle dynamics
The efficiency of particle simulation strongly depends on advancement of dynamical quantities. In the upgraded GTC we have developed particle pushers for both fully kinetic particles and guiding center particles using cylindrical coordinate in a global toroidal geometry. The physics model for the fully kinetic dynamics, guiding center dynamics and the numerical methods associated with the time advancement of the physical quantities (particle position and guiding center) are described in the following sections.
iv.1 Fully kinetic particle dynamics
Fully kinetic particle dynamics is described by the six dimensional Vlasov equation
where is the fully kinetic particle distribution function, is the particle charge, and is the particle mass.
The evolution of the particle distribution function can be described by the Newtonian equation of motion in the presence of self-consistent electromagnetic field as follows:
In our simulation we compute the marker particle trajectory Eq. (24) by the time centered Boris push method kuley2013verification (); kuley2015verification (); Boris (); Birdsall (); Tajima (); wei2015method () as discussed in section 4.2. The Lagrangian for the single particle motion in the cylindrical coordinates is written as
where is the magnetic vector potential. Now the components of generalized momenta are
In our simulation we use the poloidal flux function , rather than the vector potential to represent the equilibrium magnetic field components for calculating the toroidal canonical angular momentum. Two constants of motion for fully kinetic dynamics in cylindrical coordinates are defined by Sturdevant2017 () :
iv.2 Boris push for fully kinetic particle dynamics
Boris scheme is the most widely used orbit integrator in explicit particle-in-cell (PIC) simulation of plasmas. In this paper we have extended our Boris push scheme in GTC from Boozer coordinates kuley2015verification () to cylindrical coordinates. This scheme offers second order accuracy while requiring only one force (or field) evaluation per step. The interplay between the PIC cycle and the Boris scheme is schematically represented in Fig.2 of Ref.kuley2015verification (). At the beginning of each cycle the position of the particles and their time centered velocity as well as the grid based electromagnetic fields are given.
At the first step, we add the first half of the electric field to update the velocity from to as follows:
One may write the components of velocity at particle position at as
where . For an orthogonal cylindrical system the above Eq. (28) can be rewritten as as follows:
where , , , and .
In the second step we consider the rotation of the velocity at time . Rotated vector can be written as
where and . The components of rotated vector become
where , , , and . In the third step, we add the other half electric acceleration to the rotated vectors to obtain the velocity at time
To update the particle position we need to recover , which can be done through the following transformation (cf. Fig.2 of Ref.kuley2015verification () dark purple arrow)
where . However, the basis vector is still unknown, since does not exist in standard leap-frog scheme. Here we use a prediction for as
After we find the velocity at time , we can update the particle position using the leap-frog scheme as
In this section, we have described the time advancement of the dynamical quantities such as velocity and position of the particle using time centered approach. However, the self-consistent simulation requires update of particle weight (representing perturbed distribution function), guiding center and electric field. In our future simulations, we will use the second order Runge-Kutta (RK) method, to advance these quantities kuley2015verification (); kuley2015nonlinear ().
iv.3 Guiding center dynamics
Guiding center particle dynamics is described by the five-dimensional phase space
where is the guiding center distribution function, is the guiding center position and is the parallel velocity. The evolution of the guiding center distribution function can be described by the following equations of guiding center motionbrizard2007foundations ():
where , and . The drift velocity , the grad-B drift velocity and curvature drift velocity are given by
In GC description following order is adopted:
where is the frequency of the mode of interest, and are the wave vectors in the parallel and perpendicular direction, respectively.Two constants of motion for guiding center dynamics are defined by:
Toroidal angular momentum
For an axisymmetric system, Eq. (39) can be rewritten in cylindrical coordinate as follows:
For guiding center particle dynamics GTC normally uses second order Runge Kutta (RK) method (cf. Fig.2 of Ref. kuley2015verification ()).
|is the at magnetic axis|
|is the at magnetic axis|
To test our integration schemes (FK and GC), we consider the collisionless single particle motion of the trapped one in the core and cross-separatrix regions using the initial conditions from Table 1 and Table 2 for DIII-D geometry, respectively. The projection of the FK and GC trapped particle orbit on R-Z plane in the core and cross-separatrix regions are shown in Fig. 3. Both integrators correctly capture the trapped particle orbit and agree well with each other. The conservation properties of our integrators are tested with two exact constant of motion, viz., kinetic energy E and toroidal angular momentum .
In numerical simulation the dynamical quantities accumulate error with time. In our integrators the magnitude of the velocity provides better accuracy than the individual components of velocity. As a result the energy conservation is better than the toroidal angular momentum. The toroidal angular momentum is conserved due to the symmetry of the magnetic field along the toroidal direction. Fig. 4 shows the relative variation of E and for FK Boris integrator over a time of 15000 cyclotron period for different time step sizes in core and cross-separatrix regions. Boris algorithm maintains adequate accuracy with 20 time steps per cyclotron period . The relative energy error arises mostly due to floating point cutoff (floating precision). Whereas, the GC second order RK demonstrate that E and can converge with 240 time steps per bounce period in core and cross-separatrix regions, where is the bounce frequency of the particle [cf. Fig.( 5)]. From these convergence studies, it is found that the FK Boris integrator provides better energy convergence than the second order Runge Kutta GC integrator. If there are high frequency electromagnetic perturbations, the condition of will set an upper bound for the time step size. In the present simulations, the particles orbit are studied without electric and perturbed magnetic field. However, previously Wei et al.wei2015method () have demonstrated the effect of electric field on trapped particle orbit (Ware pinch) in the core region of tokamak using Boozer coordinates.
V Ion orbit loss near plasma edge
In a tokamak plasma, ions have drift motions due to the gradient-B and curvature drift. As a result, ion orbits are shifted from a magnetic surface. Due to this shift of ion orbit from magnetic surface, the hot ions that exist close to the separatrix, can pass near the X point region. In this region, the poloidal magnetic field is very weak and the ions have a very small poloidal displacement in time. These ions experience vertical curvature and grad-B drifts and move towards the divertor, resulting in ion orbit loss miyamoto1996direct (); stacey2013effect (); stacey2011effect (). In an axisymmetric configuration, the drift orbit of ion is obtained by solving the following three equations as described in Sec 4.3. These are the conservation of kinetic energy, magnetic moment and canonical angular momentum.
The subscript ‘a’ means the value at the starting point of the ions. The conservation of above quantities with gives the drift orbit surface as follows miyamoto1996direct ():
where is the initial pitch angle. In the above equation, we assume that , and use the approximation . A minus sign must be used instead of a plus sign in Eq. (47) when . We also consider that the poloidal magnetic field is directed clockwise and the ion toroidal drift is away from the X point. The coordinates of the X point are and the poloidal magnetic flux function of X point satisfies and . When the initial position of a test ion is at the outer midplane, that is there are following two cases for orbit ion loss miyamoto1996direct ():
|(1.2748, 0.0)||(1.2748, 0.0)|
|(0.6134, 0.0)||(0.6134, 0.0)|
|(0.8442, -0.6680)||(0.6253, -0.6680)|
|(5.253, 4.688)||(5.138, 4.688)|
Case-I: Trapped ion loss ( changes its sign before escape). The necessary conditions for escape are
where is the radial coordinate of the innermost point along the separatrix, is the null point of the separatrix for the drift surface, and
This case is represented in Fig. 6(a). In this case a test ion starting from the outer midplane initial position (green square) with initial pitch angle is reflected at M, changes the sign of and escapes near the point to the outside divertor plate. Case-I(2) is the condition for the mirror reflection of test ion at and Case-I(3) is the condition for that ion to escape through the region near point to the outer divertor plate.
Case-II: Passing ion loss ( does not change its sign before escape). The necessary conditions for escape are
This case is shown in Fig. 6(b). In this case a test ion starting from the outer midplane initial position (green square) with initial pitch angle escapes without mirror reflection near the point to the inside divertor plate as described by the condition Case-II(2).
When and , the relative position of magnetic surface and drift orbit surface [cf. Eq. (48)], clearly indicates that direct ion orbit loss is not possible. However, when the initial position of the test ion is selected as the following Case-III is possible for ion orbit loss.
Case-III: Passing ion loss. The necessary conditions for escape are
This case is described in Fig. 6(c). In this case a test ion starting from the inner midplane initial position (green square) with initial pitch angle escapes without mirror reflection near the X point to the outside divertor plate.
In the above GTC simulations (cf. Fig. 6) we did not include the radial electric field on ion orbit loss. However, one must mention here that in a tokamak scenario the radial electric field shifts the velocity space boundaries separating trapped and passing orbits. Also the ion orbit loss contributes to the generation of radial electric field, which affects the ion orbit loss. Therefore, the ion orbit loss and the radial electric field should be considered self-consistently for more accurate calculation. Similarly, it is important to study the effect of ion orbit loss and X point loss for understanding the plasma transport stacey2011effect (), particle distribution stacey2013effect (), plasma rotation pan2014co (), etc. Simulation related to these effects will be reported in a future work. To see how the electric field affects the velocity space boundaries, we have calculated the minimum velocity of the ion and particle loss fractions for DIII-D in the following paragraph.
where . For passing ion orbit loss the sign of the second term of the numerator in Eq. (51) is positive and for the trapped ion orbit loss it is negative, respectively. We solve Eq. (51) to find the ion orbit loss region for DIII-D model parameters as described in Table 3. Here, we examine the loss of ion for which . The loss region in the initial velocity space of an ion starting from outer midplane point is described by the solid curves in Fig. 7(a). Case-I and Case-II represent boundary curves CD (CD) and AB (AB) of the loss region without (with ), respectively. Boundary curves BC (BC) i.e, the condition for mirror reflection of the test ion without (with ) is given by Case-I(ii),where and correspond to the outward and inward radial electric field, respectively. The value of minimum energy is obtained from Eq. (51) [cf. Fig. 7(a)] without demonstrate good agreements with the simulation results of Fig. 6(a) and Fig. 6(b).
After the minimum loss energy is determined, particle loss fraction due to direct ion orbit loss is calculated. The ion loss rate is determined by the rate of supply of ions to the loss regions. The corresponding cumulative particle loss fraction of ions, which follow distribution function in velocity space is defined as miyamoto1996direct ():
For a Maxwellian distribution function the above equation turns out to be, , where is the energy corresponding to the minimum velocity for which ion orbit loss is possible [cf. Eq. (51)], is the ion temperature, and is the upper incomplete gamma function of order 3/2. The dependent cumulative particle loss fraction starting from without