Axisymmetric Tokamak Research Articles

Tokamak plasma equilibria with n=1 toroidal asymmetry

A general approach of how to construct plasma equilibrium in a tokamak with n=1 violation of toroidal symmetry is proposed. For an arbitrary axisymmetric tokamak plasma equilibrium, there exists the small n=1 deformation of the initial magnetic configuration that keeps the nesting of the magnetic surfaces (as in the initial configuration) and provides plasma equilibrium; such deformation and final equilibrium configuration are calculated analytically. The asymmetric analogue of the Solov'ev's equilibrium with non-degenerated plasma pressure and current density profiles is presented as an example of the application of the developed algorithm.

Control-Oriented Free-Boundary Equilibrium Solver for Tokamaks

A free-boundary equilibrium solver for an axisymmetric tokamak geometry was developed based on the finite difference method and Picard iteration in a rectangular computational area. The solver can run either in forward mode, where external coil currents are prescribed until the converged magnetic flux function ψ(R,Z) map is achieved, or in inverse mode, where the desired plasma boundary, with or without an X-point, is prescribed to determine the required coil currents. The equilibrium solutions are made consistent with prescribed plasma parameters, such as the total plasma current, poloidal beta, or safety factor at a specified flux surface. To verify the mathematical correctness and accuracy of the solver, the solution obtained using this numerical solver was compared with that from an analytic fixed-boundary equilibrium solver based on the EAST geometry. Additionally, the proposed solver was benchmarked against another numerical solver based on the finite-element and Newton-iteration methods in a triangular-based mesh. Finally, the proposed solver was compared with equilibrium reconstruction results in DIII-D experiments.

Divertor heat load estimates on NSTX and DIII-D using new and open-source 2D inversion analysis code

A thermography inversion algorithm has been developed in the open-source Python-based computer code, HYPERION, to calculate the heat flux incident on plasma-facing components (PFCs) in axisymmetric tokamaks. The chosen mesh size at the surface significantly affects the calculated transient heat flux results. The calculated transient heat flux will exceed the real value when the mesh size tends to zero but will underestimate the real value when the mesh size is large. A criterion for determining the appropriate mesh size for the transient heat flux calculation will be discussed. The numerical scheme for HYPERION uses a 2D fully implicit finite-difference approach, allowing temperature-dependent thermal properties of PFC materials. The inversion algorithm is benchmarked against established heat flux calculation codes, TACO and THEODOR, based on thermography data from NSTX and DIII-D respectively. The primary benefits of HYPERION compared to TACO and THEODOR are that it is open-source and it allows for the optimization of mesh thickness along the substrate. The algorithm also accounts for the thermal properties of thin surface layers that characteristically form on PFCs due to plasma-material interactions. The agreement between HYPERION and THEODOR is excellent, as the percent difference between the codes is ∼5% on average in the case of the DIII-D data for moderate to high heat flux. Verification tests with TACO show slightly higher average percent differences of 8% and 12%. In using HYPERION to study filaments in heat flux, the initial results indicate that small ELMs filaments significantly broaden the divertor heat flux, and decrease divertor peak flux. Compared to the inter-ELM, the small ELM filaments decrease the divertor peak surface temperature. With intermittent divertor filaments, the divertor heat flux width is comparable with that found in L-mode.

Iterative addition of parallel non-local effects to full wave ICRF finite element models in axisymmetric tokamak plasmas

The current response of a hot magnetized plasma to a radio-frequency wave is non-local, turning the electromagnetic wave equation into an integro-differential equation. Non-local physics gives rise to wave physics and absorption processes not observed in local media. Furthermore, non-local physics alters wave propagation and absorption properties of the plasma. In this work, an iterative method that accounts for parallel non-local effects in 2D axisymmetric tokamak plasmas is developed, implemented, and verified. The iterative method is based on the finite element method and Fourier decomposition, with the advantage that this numerical scheme can describe non-local effects while using a high-fidelity antenna and wall representation, as well as limiting memory usage. The proposed method is implemented in the existing full wave solver FEMIC and applied to a minority heating scenario in ITER to quantify how parallel non-local physics affect wave propagation and dissipation in the ion cyclotron range of frequencies (ICRF). The effects are then compared to a reduced local plane wave model, both verifying the physics implemented in the model, as well as estimating how well a local plane wave approximation performs in scenarios with high single pass damping. Finally, the new version of FEMIC is benchmarked against the ICRF code TORIC.

Simulations of tokamak edge plasma turbulent fluctuations based on a minimal 3D model

A new simulation model for tokamak boundary plasma, SOLT3D, is implemented in the BOUT++ framework (Dudson et al 2009 Comput. Phys. Commun. 180 1467). The simulation model includes a set of dynamic equations describing collisional boundary plasma and neutral gas in the tokamak scrape-off layer and divertor region. The model is verified against standard linear plasma instabilities and available nonlinear results. For L-mode like conditions, SOLT3D reproduces characteristics of boundary plasma turbulent fluctuations that are similar to published experimental data, in terms of the amplitude and spatial dependence of the fluctuations. It also reproduces realistic plasma fluxes on material surfaces and expected Bohm-like effective radial transport. Plasma fluctuations inferred from the simulations imply inevitably a significant level of intrinsic ‘noise’ for axisymmetric tokamak plasma transport modeling, introducing errors on the order of unity. In particular, the toroidally averaged atomic rates below 5–10 eV are strongly modified by turbulent plasma fluctuations, which should significantly affect the standard axisymmetric modeling of the tokamak edge plasma and divertor.

Diagnostic weight functions in constants-of-motion phase-space

The fast-ion phase-space distribution function in axisymmetric tokamak plasmas is completely described by the three constants of motion: energy, magnetic moment and toroidal canonical angular momentum. In this work, the observable regions of constants-of-motion phase-space, given a diagnostic setup, are identified and explained using projected velocities of the fast ions along the diagnostic lines-of-sight as a proxy for several fast-ion diagnostics, such as fast-ion Dα spectroscopy, collective Thomson scattering, neutron emission spectroscopy and gamma-ray spectroscopy. The observable region in constants-of-motion space is given by a position condition and a velocity condition, and the diagnostic sensitivity is given by a gyro-orbit and a drift-orbit weighting. As a practical example, 3D orbit weight functions quantifying the diagnostic sensitivity to each point in phase-space are computed and investigated for the future COMPASS-Upgrade and MAST-Upgrade tokamaks.

Nonlinear equilibria and transport processes in burning plasmas

In this work, we put forward a general phase space transport theory in axisymmetric tokamak plasmas based upon the concept of zonal state (ZS). Within this theoretical framework, the ZS corresponds to a renormalized plasma nonlinear equilibrium consisting of phase space zonal structures (PSZS) and zonal electromagnetic fields (ZFs) which evolve self-consistently with symmetry breaking fluctuations and sources/collisions. More specifically, our approach involves deriving governing equations for the evolution of particle distribution functions (i.e, PSZS), which can be used to compute the corresponding macro-/meso-scale evolving magnetized plasma equilibrium adopting the Chew Goldberger Low description, separating the spatiotemporal microscale structures. The nonlinear physics of ZFs and of geodesic acoustic modes (GAMs)/energetic particle driven GAMs is then analyzed to illustrate the applications of our theory.

Coupling finite elements of class C1 on composite curved meshes for second order elliptic problems

SummaryFinite elements of class are suitable for the computation of magnetohydrodynamics instabilities in tokamak plasmas. In addition, isoparametric approximations allow for a precise alignment of the mesh with the magnetic field line. Mesh alignment is crucial to achieve axisymmetric equilibria accurately. It is also helpful to deal with the anisotropy nature of magnetized plasma flows. In this numerical framework, several practical simulations are now available. They help to understand better the operation of existing devices and predict the optimal strategies for using the international ITER tokamak under construction. However, a mesh‐aligned isoparametric representation suffers from the presence of critical points of the magnetic field (magnetic axis, X‐point). We here explore a strategy that combines aligned mesh out of the critical points with non‐aligned unstructured mesh in a region containing these points. By this strategy, we can avoid highly stretched elements and the numerical difficulties that come with them. The mesh‐aligned interpolation uses bi‐cubic Hemite‐Bézier polynomials on a structured mesh of curved quadrangular elements. On the other hand, we assume reduced cubic Hsieh‐Clough‐Tocher finite elements on an unstructured triangular mesh. Both meshes overlap, and the resulting formulation is a coupled discrete problem solved iteratively by a suitable one‐level Schwarz algorithm. In this paper, we will focus on the Poisson problem on a two‐dimensional bounded regular domain. This elliptic equation is a simplified version of the axisymmetric tokamak equilibrium one at the asymptotic limit of infinite major radius (large aspect ratio).

Distribution transforms for guiding center orbit coordinates in axisymmetric tokamak equilibria

Accurate distribution transforms between common 4D axisymmetric tokamak guiding center coordinates {energy,pitch,R,Z}, and orbit-space coordinates (made up of a particle's kinetic energy E, maximum radius in the tokamak Rm, pitch at the point of maximum radius pm, and time-normalised phase τ=t/torbit) are vital for the verification and implementation of orbit tomography. These transforms have proven difficult in the past due to discontinuities across topological boundaries in orbit-space. In this work we exhaustively investigate these transforms, comparing existing and novel transform methods in speed, accuracy and specific sources of error. A distribution transform that samples {energy,pitch,R,Z}-space along pre-computed orbit paths, and relies on Jacobians calculated with autodifferentiation, is benchmarked against a Monte Carlo sampling and binning algorithm. Our new transform demonstrates a hundredfold increase in speed, and better preserves natural discontinuities across the orbit trapped-passing boundary. A potentially damaging source of transform error caused by a Jacobian singularity that occurs for vanishingly small orbits is addressed, ensuring repeated transforms are well-behaved. The application of smoothing splines in {energy,pitch,R,Z}-space and orbit-space is also discussed. A Jacobian-based transform utilising thin-plate polyharmonic splines restricted to subdomains of similar orbit class is presented and benchmarked against its equivalent non-splined transform. This new smoothing transform correctly avoids interpolating across the trapped/passing boundary, doing so without the prohibitively slow computation and hyperparameter tuning required by previous orbit-space splines.

Testing the conservative character of particle simulations: II. Spurious heating of guiding centers and full orbits subject to fluctuations expressed in terms of E and B

For an axisymmetric tokamak plasma, Hamiltonian theory predicts that the orbits of charged particles must stay on invariant tori of conserved energy in the moving frame of reference of a wave that propagates along the torus with fixed angular phase velocity, amplitude, and shape. The mode structure in the poloidal plane is arbitrary if the fluctuations are expressed in terms of potentials Φ and A, which satisfy Faraday's law and the solenoidal condition by definition. Consequently, smoothing operations (such as gyroaveraging and noise suppression) do not violate the conservative laws. However, this is not guaranteed for models expressed in terms of the physical fields E and B. Here, we demonstrate that manipulations of E and B in the poloidal (R, z) plane can cause spurious heating that is independent of time steps or numerical methods, but can be sensitive to geometry. In particular, we show that secular acceleration is enhanced when one imposes nonnormal modes that possess strong up–down asymmetry instead of the usual in–out asymmetry of normal toroidal (eigen)modes. We compare full gyro-orbit and guiding center models and find similar behavior. We also examine the effect of ad hoc N-point gyroaveraging in a guiding center model, as is done in some simulation codes. If one uses Faraday's law to (re)compute B(t) after gyroaveraging E, the guiding center motion remains conservative. Otherwise, spurious heating should be expected and monitored, but it may be tolerable when normal modes dominate.

Electron-cyclotron resonance heating and current drive source for flux-driven gyrokinetic simulations of tokamaks

Electromagnetic waves that resonate with the cyclotron motion of electrons in a magnetized plasma can efficiently transfer their momentum and energy to the plasma. This is routinely used to heat or drive current in tokamak plasmas. The impact of this localized source of momentum and energy on turbulence and the retro-action of turbulence on the resonant interaction between the electromagnetic wave and the plasma has been scarcely studied due to the difficulty in self-consistently simulating the two physical mechanisms. In this paper, a realistic source representing electron-cyclotron resonance heating (ECRH) and electron-cyclotron current drive (ECCD) is derived and implemented in a gyrokinetic code. The implementation of this realistic source in any existing global gyrokinetic code would enable the self-consistent study of turbulence in the presence of ECRH/ECCD using this code. The analytical source derived in this paper is valid for a beam propagating in the equatorial plane of an axisymmetric tokamak plasma. The realistic ECRH/ECCD source is implemented in the global gyrokinetic code ORB5 and successfully benchmarked against analytical theory (Albajar et al 2006 Plasma Phys. Control. Fusion 49 15–29) and the C3PO/LUKE suite of codes (Peysson et al 2011 Plasma Phys. Control. Fusion 53 124028), which is routinely used to study ECRH/ECCD deposition.

Representation and modeling of charged particle distributions in tokamaks

Experimental diagnostics, analysis tools and simulations represent particle distributions in various forms and coordinates. Algorithms to manage these data are needed on platforms like the ITER Integrated Modeling & Analysis Suite (IMAS), performing tasks such as archiving, modeling, conversion and visualization. A method that accomplishes some of the required tasks for distributions of charged particles with arbitrarily large magnetic drifts in axisymmetric tokamak geometry is described here. Given a magnetic configuration, we first construct a database of guiding center orbits, which serves as a basis for representing particle distributions. The orbit database contains the geometric information needed to perform conversions between arbitrary coordinates, modeling tasks, and resonance analyses. Using that database, an imported or newly modeled distribution is mapped to an exact equilibrium, where the dimensionality is reduced to three constants of motion (CoM). The orbit weight is uniquely given when the input is a true distribution: one that measures the number of physical particles per unit of phase space volume. Less ideal inputs, such as distributions estimated without drifts, or models of particle sources, can also be processed. As an application example, we reconstruct the drift-induced features of a distribution of fusion-born alpha particles in a large tokamak, given only a birth profile, which is not a function of the alpha's CoM. Repeated back-and-forth transformations between CoM space and energy-pitch-cylinder coordinates are performed for verification and as a proof of principle for IMAS.

An asymptotic-preserving 2D-2P relativistic drift-kinetic-equation solver for runaway electron simulations in axisymmetric tokamaks

We propose an asymptotic-preserving (AP), uniformly convergent numerical scheme for the relativistic collisional Drift-Kinetic Equation (rDKE) to simulate runaway electrons in axisymmetric toroidal magnetic field geometries typical of tokamak devices. The approach is derived from an exact Green's function solution with numerical approximations of quantifiable impact, and results in a simple, two-step operator-split algorithm, consisting of a collisional Eulerian step, and a Lagrangian orbit-integration step with analytically prescribed kernels. The AP character of the approach is demonstrated by analysis of the dominant numerical errors, as well as by numerical experiments. We demonstrate the ability of the algorithm to provide accurate answers regardless of plasma collisionality on a circular axisymmetric tokamak geometry.

Demonstration of reduced neoclassical energy transport in Wendelstein 7-X

Research on magnetic confinement of high-temperature plasmas has the ultimate goal of harnessing nuclear fusion for the production of electricity. Although the tokamak1 is the leading toroidal magnetic-confinement concept, it is not without shortcomings and the fusion community has therefore also pursued alternative concepts such as the stellarator. Unlike axisymmetric tokamaks, stellarators possess a three-dimensional (3D) magnetic field geometry. The availability of this additional dimension opens up an extensive configuration space for computational optimization of both the field geometry itself and the current-carrying coils that produce it. Such an optimization was undertaken in designing Wendelstein 7-X (W7-X)2, a large helical-axis advanced stellarator (HELIAS), which began operation in 2015 at Greifswald, Germany. A major drawback of 3D magnetic field geometry, however, is that it introduces a strong temperature dependence into the stellarator’s non-turbulent ‘neoclassical’ energy transport. Indeed, such energy losses will become prohibitive in high-temperature reactor plasmas unless a strong reduction of the geometrical factor associated with this transport can be achieved; such a reduction was therefore a principal goal of the design of W7-X. In spite of the modest heating power currently available, W7-X has already been able to achieve high-temperature plasma conditions during its 2017 and 2018 experimental campaigns, producing record values of the fusion triple product for such stellarator plasmas3,4. The triple product of plasma density, ion temperature and energy confinement time is used in fusion research as a figure of merit, as it must attain a certain threshold value before net-energy-producing operation of a reactor becomes possible1,5. Here we demonstrate that such record values provide evidence for reduced neoclassical energy transport in W7-X, as the plasma profiles that produced these results could not have been obtained in stellarators lacking a comparably high level of neoclassical optimization.

Evolution of Heavy Ion Beam Probing from the Origins to Study of Symmetric Structures in Fusion Plasmas

The overview discusses development of the unique fusion plasma diagnostics—Heavy Ion Beam Probing (HIBP) in application to toroidal magnetic plasma devices. The basis of the HIBP measurements of the plasma electric potential and processing of experimental data are considered. Diagnostic systems for probing plasma in tokamaks TM-4, TJ-1, TUMAN-3M and T-10, stellarators WEGA, TJ-II and Uragan-2M are presented. Promising results of the HIBP projects for various existing modern machines, such as TCV, TCABR, MAST, COMPASS, GLOBUS-M2, T-15 MD and W7-X and the international fusion tokamak reactor ITER are given. Results from two machines with similar size and plasma parameters, but with different types of the magnetic con-figuration: axisymmetric tokamak T-10 and helically symmetric stellarator TJ-II are compared. The results of studies of stationary potential profiles and oscillations in the form of quasimonochromatic and broadband fluctuations, turbulent particle flux, fluctuations of density and poloidal magnetic field are presented. The properties of symmetric structures—zonal flows and geodesic acoustic modes of plasma oscillations as well as Alfvén Eigenmodes excited by fast particles from neutral beam injection heating are described. General trends in the behavior of electric potential and turbulence in magnetized fusion plasmas are revealed.

Parametric dependencies of locked mode thresholds in KSTAR L-mode plasmas

Resonant error fields with toroidal harmonics n = 1 and n = 2 in otherwise axisymmetric tokamak plasmas can drive disruptive locked modes when their strength exceeds a threshold value. This paper presents the sensitivity of error field threshold in global parameters such as line-averaged density and toroidal magnetic field as well as toroidal harmonics dependence in KSTAR L-mode plasmas. The error field threshold for mode-locking increases along with density up to a roll-over density, which coincides well with the linear ohmic confinement (LOC) to saturated ohmic confinement (SOC) density for the studied KSTAR plasmas. This non-monotonic density correlation of the error field threshold in LOC–SOC can be explained by a modification of existing theoretical scaling in SOC, while plasma rotation can also play an important role in setting the resonant error field threshold.

The development of an implicit full f method for electromagnetic particle simulations of Alfvén waves and energetic particle physics

In this work, an implicit scheme for particle-in-cell/Fourier electromagnetic simulations is developed and applied to studies of Alfvén waves in one dimension and three dimensional tokamak plasmas. An analytical treatment is introduced to achieve efficient convergence of the iterative solution of the implicit field-particle system. First, its application to the one-dimensional uniform plasma demonstrates its applicability in a broad range of β/me values. Second, toroidicity induced Alfvén eigenmodes (TAE) are simulated in a three dimensional axisymmetric tokamak plasma, using the widely studied case defined by the International Tokamak Physics Activity (ITPA) Energetic Particle (EP) Topical Group. The real frequency and the growth (or damping) rate of the TAE with (or without) EPs agree with previous results reasonably well. The full f electromagnetic particle scheme established in this work provides a possible natural choice for EP transport studies where large profile variation and arbitrary particle distribution functions need to be treated in kinetic simulations.

Quasi-geometric integration of guiding-center orbits in piecewise linear toroidal fields

A numerical integration method for guiding-center orbits of charged particles in toroidal fusion devices with three-dimensional field geometry is described. Here, high order interpolation of electromagnetic fields in space is replaced by a special linear interpolation, leading to locally linear Hamiltonian equations of motion with piecewise constant coefficients. This approach reduces computational effort and noise sensitivity, while the conservation of total energy, magnetic moment and phase space volume is retained. The underlying formulation treats motion in piecewise linear fields exactly and, thus, preserves the non-canonical symplectic form. The algorithm itself is only quasi-geometric due to a series expansion in the orbit parameter. For practical purposes, an expansion to the fourth order retains geometric properties down to computer accuracy in typical examples. When applied to collisionless guiding-center orbits in an axisymmetric tokamak and a realistic three-dimensional stellarator configuration, the method demonstrates stable long-term orbit dynamics conserving invariants. In Monte Carlo evaluation of transport coefficients, the computational efficiency of quasi-geometric integration is an order of magnitude higher than with a standard fourth order Runge–Kutta integrator.

Peeling-ballooning stability of tokamak plasmas with applied 3D magnetic fields

The poloidal harmonics of the toroidal normal modes of an unstable axisymmetric tokamak plasma are employed as basis functions for the minimisation of the 3D energy functional. This approach presents a natural extension of the perturbative method considered in Anastopoulos Tzanis et al (2019 Nucl. Fusion 59 126028). This variational formulation is applied to the stability of tokamak plasmas subject to external non-axisymmetric magnetic fields. A comparison of the variational and perturbative methods shows that for D-shaped, high βN plasmas, the coupling of normal modes becomes strong at experimentally relevant applied 3D fields, leading to violation of the assumptions that justify a perturbative analysis. The variational analysis employed here addresses strong coupling, minimising energy with respect to both toroidal and poloidal Fourier coefficients. In general, it is observed that ballooning unstable modes are further destabilised by the applied 3D fields and field-aligned localisation of the perturbation takes place, as local ballooning theory suggests. For D-shaped high βN plasmas, relevant to experimental cases, it is observed that the existence of intermediate n unstable peeling-ballooning modes, where a maximum in the growth rate spectrum typically occurs, leads to a destabilising synergistic coupling that strongly degrades the stability of the 3D system.

Effect of non-axisymmetric perturbations on the ambipolar Er and neoclassical particle flux inside the ITER pedestal region

The transport dynamics of impurities in the pedestal region of ITER plasmas is of crucial interest since this regulates the penetration of impurities from the edge into the core plasma, where an excessive accumulation of impurities can degrade their fusion performance. In the pedestal region of H-mode tokamak plasmas, anomalous transport is highly reduced and impurity transport is found to be well described by the neoclassical theory. Under these conditions, perturbations to the axisymmetric tokamak geometry can strongly affect both the radial electric field and particle transport. In this work, we describe the results of numerical studies performed to quantify the effects on the pedestal ambipolar electric field and radial particle fluxes of the non-axisymmetric fields, associated with both the intrinsic toroidal field ripple and extrinsic fields applied for ELM control, for ITER Q = 10 plasma conditions with emphasis on high Z impurity transport. It is found that the effect of the ITER toroidal field ripple on high Z impurity transport is negligible. On the contrary, extrinsic three-dimensional fields applied for ELM control cause a strong modification of the pedestal ambipolar electric (to less negative values) and the appearance of multi-valued solutions for the pedestal electric field, analogue to core stellarator transport significantly increasing the outward character of neoclassical pedestal transport for both the main plasma ions (D and T in ITER) and high Z impurities, suggesting a strong modification of the background plasma profiles. Finally, it is found that for the Z impurity, its quantitative evaluation has uncertainties (with important implications for the radial flow direction) associated with the high poloidal Mach number ~ 1, due to the high pedestal electric field and large ion mass, which renders the first order neoclassical theory questionable. Detailed ITER MHD 3D equilibria and a benchmark of SFINCS against NEO and NCLASS codes has also been obtained.