Session GI2: Advances in Plasma Simulation I

3D Modeling of the Sawtooth Instability in a Small Tokamak

The sawtooth instability [1] is the most fundamental dynamic of an inductive tokamak discharge such as will occur in ITER. Sawtooth behavior is complex and remains incompletely explained. While the instability is confined to the center of the plasma in low-pressure, low-current, large aspect ratio discharges, under certain conditions it can create magnetic islands at the outer resonant surfaces and may set off a sequence of events that leads to a major disruption. Under some circumstances the reconnection following the sawtooth is observed to be complete; in others, it is incomplete. As part of the CEMM SciDAC project, we have undertaken an ambitious campaign to model this periodic motion as accurately as possible using the most complete fluid- like description of the plasma, the Extended MHD model. Both NIMROD and M3D have been applied to this problem, and we are also using it as a non-trivial test problem to compare these two codes far into the nonlinear regime. Compared to the MHD model, Extended MHD predicts plasma rotation, faster reconnection, and reduced field line stochasticity in the crash aftermath. The multiple time scales associated with the reconnection layer and growth time make this an extremely challenging computational problem. A recent M3D simulation used over 500,000 elements for 400,000 partially implicit time steps, and there still remain some resolution issues. However these calculations are providing insight into the nonlinear mechanisms of surface breakup and healing. We have been able to match many features of a small tokamak and can now project to the computational requirements for simulations of larger, hotter devices such as ITER. These simulations form the basis for studying more complex phenomena such as the effect on these modes of an energetic particle component, or of externally generated electromagnetic waves (RF). [1] R.J. Hastie, Astrophys. Space Sci. {\bf 256} 177 (1997).   

Particle simulation of bursting Alfv\'{e}n modes in JT-60U

Periodic radial redistributions of energetic ions in NNB heated discharges in JT-60U appear to be caused by Abrupt Large-amplitude Events (ALEs): bursting modes, with hundred-microsecond time scales. Between two successive ALEs, weaker fast Frequency Sweeping (FS) modes occur, with millisecond time scales. Particle simulations of a typical NNB-heated JT-60U discharge is presented. Assuming the energetic-ion density profile experimentally observed immediately before an ALE and postulating an anisotropic slowing-down distribution function, we find that a fast-growing mode is destabilized at half radius. Its saturation is accompanied by a macroscopic outward displacement of energetic ions. Frequency range, time scale and broad resonance region support the identification of such a mode with an ALE. The main features of the quiescent phase observed between two ALEs (namely, the destabilization of low-growth-rate, low-amplitude fast FS modes) are then recovered if we initialize energetic ions according to the distribution function modified by the former bursting mode, caring to include both the density profile relaxation and the distortion of the velocity space distribution function. The experiments could then be interpreted as follows: given the energetic-ion source provided by beam injection, the free energy reconstruction rate is set by the rebuilding of the velocity space distribution function. The intermediate fast ion phase-space profiles between two ALEs are characterized by lower mode drive than that of a slowing-down distribution with the same energetic-ion density profile. Weak fast FS modes are then excited, which are unable to contrast the density profile reconstruction. As soon as the combined restoration of the configuration and velocity space distribution provides enough drive for a fast growing Alfv\'{e}n mode, a new ALE occurs.   

Two-Fluid Physics and Field Reversed Configurations

Fluid models of plasmas are a common tool to study fusion devices. In this talk algorithms for the solution of Two-Fluid plasma equations are presented and applied to the study of Field Reversed Configurations (FRCs). The Two-Fluid model is more general than the often used magnetohydrodynamic (MHD) model. The model takes into account electron inertia, charge separation and the full electromagnetic field equations and allows for separate electron and ion motion. Finite Lamor Radii effects are taken into account by self consistently evolving the anisotropic pressure tensor. The algorithm presented is the high resolution wave propagation scheme. The wave propagation method is based on solutions to the Riemann problem at cell interfaces. Operator splitting is used to incorporate the Lorentz and electromagnetic source terms. To preserve the divergence constraints on the electric and magnetic fields the so called perfectly-hyperbolic form of Maxwell equations are used which explicitly incorporate the divergence equations into the time stepping scheme. A detailed study of Field-Reversed Configuration stability and formation is performed. The study is divided into two parts. In the first, FRC stability is studied. The simulation is initialized with various FRC equilibria and perturbed. The growth rates are calculated and compared with MHD results. It is shown that the FRCs are indeed more stable within the Two-Fluid model than the MHD model. In the second part formation of FRCs is studied. In this set of simulations a cylindrical column of plasma is initialized with a uniform axial magnetic field. The field is reversed at the walls. Via the process of magnetic reconnection FRC formation is observed. The effects of Rotating Magnetic Field (RMF) drive on the formation of FRC are also presented. Here, a set of current carrying coils apply a RMF at the plasma boundary, causing a electron flow in the R-Z plane leading to field reversal. The strong azimuthal electron flow causes Lower-Hybrid Drift Instabilities (LHDI), which can be captured if the ion-gyroradius is well resolved. The LHDI is known to be a possible source of anomalous resistivity in many plasma configurations.   

Multidimensional kinetic simulations using dissipative closures and other \textit{reduced} Vlasov methods for differing particle magnetizations

Kinetic plasma simulations in which the phase-space distribution functions are advanced directly via the coupled Vlasov and Poisson (or Maxwell) equations---better known simply as \textit{Vlasov} simulations---provide a valuable low-noise complement to the more commonly employed Particle-in-Cell (PIC) simulations. However, in more than one spatial dimension Vlasov simulations become numerically demanding due to the high dimensionality of $\mathbf{x}$--$\mathbf{v}$ phase-space. Methods that can reduce this computational demand are therefore highly desirable. Several such methods will be presented, which treat the phase-space dynamics along a dominant dimension (e.g., parallel to a beam or current) with the full Vlasov propagator, while employing a \textit{reduced} description, such as moment equations, for the evolution perpendicular to the dominant dimension. A key difference between the moment-based (and other reduced) methods considered here and standard fluid methods is that the moments are now functions of a phase-space coordinate (e.g. moments of $v_y$ in $z$--$v_z$--$y$ phase space, where $z$ is the dominant dimension), rather than functions of spatial coordinates alone. Of course, moment-based methods require closure. For effectively unmagnetized species, new dissipative closure methods inspired by those of Hammett and Perkins [\textit{PRL}, \textbf{64}, 3019 (1990)] have been developed, which exactly reproduce the linear electrostatic response for a broad class of distributions with power-law tails, as are commonly measured in space plasmas. The nonlinear response, which requires more care, will also be discussed. For weakly magnetized species (i.e., $\Omega_s<\omega_s$) an alternative algorithm has been developed in which the distributions are assumed to gyrate about the magnetic field with a fixed nominal perpendicular ``thermal'' velocity, thereby reducing the required phase-space dimension by one. These reduced algorithms have been incorporated into 2-D codes used to study the evolution of nonlinear structures such as double layers and electron holes in Earth's auroral zone.   

A Theory-Based Transport Model With Comprehensive Physics

A new theory based transport model with comprehensive physics (trapping, general toroidal geometry, finite beta, collisions) has been developed. The core of the model is the new trapped-gyro-Landau-fluid (TGLF) equations [1] which provide a fast and accurate approximation to the linear eigenmodes of drift-wave instabilities (trapped ion and electron modes, ion and electron temperature gradient modes and kinetic ballooning modes). This new TGLF transport model employs several new technologies that remove the limitations of its predecessor GLF23. No fitting to experiment is done so applying the model to experiments is a true test of the theory it is approximating. A model for the averaging of the Landau resonance by the trapped particles makes the equations work seamlessly over the whole drift-wave wavenumber range. A fast eigenmode solution method enables unrestricted magnetic geometry. A new model for electron-ion collisions and both parallel and perpendicular electromagnetic fluctuations are included. The linear eigenmodes have been benchmarked against comprehensive physics gyrokinetic calculations over a large range of plasma parameters. Deviation between the gyrokinetic and TGLF linear growth rates averages 11.4\% in shifted circle geometry [1]. The transport model uses the TGLF eigenmodes to compute quasilinear fluxes of energy and particles. A model for the saturated turbulence amplitude is fitted to a large set of non-linear GYRO simulations. The fluxes at each wavenumber are fit by the model. The TGLF model is valid for the low aspect ratio spherical torus which has both a high trapped fraction and strong shaping of magnetic flux surfaces. It is also valid close to the magnetic separatrix so the transport physics of the H-mode pedestal region can be explored.\par \vskip6pt \noindent [1]~G.M. Staebler, et al., Phys.\ Plasmas {\bf 12}, 102508 (2005).   

Coupled ITG/TEM-ETG Gyrokinetic Simulations

Electron temperature gradient (ETG) transport is conventionally defined as the electron energy transport from high-k where ions are adiabatic and there can be no ion energy or plasma transport. Previous simulations have assumed adiabatic ions (ETG-ai). However using the GYRO gyrokinetic code [1], we have found that many simulation cases with trapped electron at moderate shear do not nonlinearly saturate unless fully kinetic ions and some low-k ion scale zonal flow modes are included. We define high-k ETG-ki transport to be that arising from $k_y\rho_{s-i} > 1$ including electron gyroradius scales $k_y\rho_{s-e} \sim 1$, and ion temperature gradient and trapped electron mode (ITG/TEM) transport to be that from $k_y\rho_{s-i} \leq 1$. There has been speculation [2] that ETG transport could be modified by the nonlinear coupling to the ITG/TEM turbulence (or vise-versa). We have done very expensive high Reynold's number $(k_{\perp-max}/k_{\perp-min})^2\propto (\rho_{s-i}/\rho_{s-e})^2=\mu^2$ high-resolution-large-flux-tube simulations with coupled ITG/TEM-ETG-ki turbulence. By comparing expensive simulations with much cheaper uncoupled high-resolution-small-flux-tube ETG-ki and low-resolution-large-flux-tube ITG/TEM simulations, we hope to demonstrate that superposition of the cheaper simulations is sufficiently accurate. Electron energy transport from ETG-ki is 5-10 fold smaller than from ITG/TEM except when the $E\times B$ shear is strong enough to quenched the low-k transport. GYRO compute time for the expensive simulations scales as $\sim\mu^{3-4}$. Reduced mass ratio $\mu=20$ simulations have been done, and $\mu=30$ simulations in progress accurately represent the $\mu=\,$40-60 physical range.\par \vskip6pt \noindent [1] J. Candy and R.E. Waltz, Phys.\ Rev.\ Lett.\ {\bf 91} (2003) 045001.\par \noindent [2] C. Holland and P.H. Diamond, Phys.\ Plasmas {\bf 11} (2004) 1043.