M. V. Umansky

BOUT++: A framework for parallel plasma fluid simulations

Abstract A new modular code called BOUT++ is presented, which simulates 3D fluid equations in curvilinear coordinates. Although aimed at simulating Edge Localised Modes (ELMs) in tokamak x -point geometry, the code is able to simulate a wide range of fluid models (magnetised and unmagnetised) involving an arbitrary number of scalar and vector fields, in a wide range of geometries. Time evolution is fully implicit, and 3rd-order WENO schemes are implemented. Benchmarks are presented for linear and non-linear problems (the Orszag–Tang vortex) showing good agreement. Performance of the code is tested by scaling with problem size and processor number, showing efficient scaling to thousands of processors. Linear initial-value simulations of ELMs using reduced ideal MHD are presented, and the results compared to the ELITE linear MHD eigenvalue code. The resulting mode-structures and growth-rate are found to be in good agreement ( γ BOUT++ = 0.245 ω A , γ ELITE = 0.239 ω A , with Alfvenic timescale 1 / ω A = R / V A ). To our knowledge, this is the first time dissipationless, initial-value simulations of ELMs have been successfully demonstrated.

The magnetic field structure of a snowflake divertor

The snowflake divertor exploits a tokamak geometry in which the poloidal magnetic field null approaches second order; the name stems from the characteristic hexagonal, snowflakelike shape of the separatrix for an exact second-order null. The proximity of the poloidal field structure to that of a second-order null substantially modifies edge magnetic properties compared to the standard X-point geometry (with a first-order null); this, in turn, affects the edge plasma behavior. Modifications include: (1) The flux expansion near the null-point becomes 2–3 times larger. (2) The connection length between the equatorial plane and divertor plate increases. (3) Magnetic shear just inside the separatrix becomes much larger. (4) In the open-field-line region, the squeezing of the flux-tubes near the null-point increases, thereby causing stronger decoupling of the plasma turbulence in the divertor legs and in the main scrape-off layer. These effects can be used to reduce the power load on the divertor plates and/or ...

Nonlinear ELM simulations based on a nonideal peeling–ballooning model using the BOUT++ code

A minimum set of equations based on the peeling–ballooning (P–B) model with nonideal physics effects (diamagnetic drift, E × B drift, resistivity and anomalous electron viscosity) is found to simulate pedestal collapse when using the BOUT++ simulation code, developed in part from the original fluid edge code BOUT. Linear simulations of P–B modes find good agreement in growth rate and mode structure with ELITE calculations. The influence of the E × B drift, diamagnetic drift, resistivity, anomalous electron viscosity, ion viscosity and parallel thermal diffusivity on P–B modes is being studied; we find that (1) the diamagnetic drift and E × B drift stabilize the P–B mode in a manner consistent with theoretical expectations; (2) resistivity destabilizes the P–B mode, leading to resistive P–B mode; (3) anomalous electron and parallel ion viscosities destabilize the P–B mode, leading to a viscous P–B mode; (4) perpendicular ion viscosity and parallel thermal diffusivity stabilize the P–B mode. With addition of the anomalous electron viscosity under the assumption that the anomalous kinematic electron viscosity is comparable to the anomalous electron perpendicular thermal diffusivity, or the Prandtl number is close to unity, it is found from nonlinear simulations using a realistic high Lundquist number that the pedestal collapse is limited to the edge region and the ELM size is about 5–10% of the pedestal stored energy. This is consistent with many observations of large ELMs. The estimated island size is consistent with the size of fast pedestal pressure collapse. In the stable α-zones of ideal P–B modes, nonlinear simulations of viscous ballooning modes or current-diffusive ballooning mode (CDBM) for ITER H-mode scenarios are presented.

Dependence of the L- to H-mode Power Threshold on Toroidal Rotation and the Link to Edge Turbulence Dynamics

The injected power required to induce a transition from L-mode to H-mode plasmas is found to depend strongly on the injected neutral beam torque and consequent plasma toroidal rotation. Edge turbulence and flows, measured near the outboard midplane of the plasma (0.85 < r/a < 1.0) on DIII-D with the high-sensitivity 2D beam emission spectroscopy (BES) system, likewise vary with rotation and suggest a causative connection. The L–H power threshold in plasmas with the ion ∇B drift directed away from the X-point decreases from 4–6 MW with co-current beam injection, to 2–3 MW near zero net injected torque and to <2 MW with counter-injection in the discharges examined. Plasmas with the ion ∇B drift directed towards the X-point exhibit a qualitatively similar though less pronounced power threshold dependence on rotation. 2D edge turbulence measurements with BES show an increasing poloidal flow shear as the L–H transition is approached in all conditions. As toroidal rotation is varied from co-current to balanced in L-mode plasmas, the edge turbulence changes from a uni-modal character to a bi-modal structure, with the appearance of a low-frequency (f = 10–50 kHz) mode propagating in the electron diamagnetic direction, similar to what is observed as the ion ∇B drift is directed towards the X-point in co-rotating plasmas. At low rotation, the poloidal turbulence flow near the edge reverses prior to the L–H transition, generating a significant poloidal flow shear that exceeds the measured turbulence decorrelation rate. This increased poloidal turbulence velocity shear appears to facilitate the L–H transition. No such reversal is observed in high rotation plasmas. The high-frequency poloidal turbulence velocity spectrum exhibits a transition from a geodesic acoustic mode zonal flow to a higher-power, lower frequency zero-mean-frequency zonal flow as rotation varies from co-current to balanced during a torque scan at constant injected neutral beam power, perhaps also facilitating the L–H transition. This reduced power threshold at lower toroidal rotation may benefit inherently low-rotation plasmas such as ITER.

Simulation of edge localized modes using BOUT

The BOUT++ code is used to simulate edge localized modes (ELMs) in a shifted circle equilibrium. Reduced ideal MHD simulations are first benchmarked against the linear ideal MHD code ELITE, showing good agreement. Diamagnetic drift effects are included finding the expected suppression of high toroidal mode-number modes. Nonlinear simulations are performed, making the assumption that the anomalous kinematic electron viscosity is comparable to the anomalous electron thermal diffusivity. This allows simulations with realistically high Lundquist numbers (S = 108), finding ELM sizes of 5–10% of the pedestal stored thermal energy. Scans show a strong dependence of the ELM size on resistivity at low Lundquist numbers, with higher resistivity leading to more violent eruptions. At high Lundquist numbers relevant to high-performance discharges, ELM size is independent of resistivity as hyper-resistivity becomes the dominant dissipative effect.

Edge Gyrokinetic Theory and Continuum Simulations

The following results are presented from the development and application of TEMPEST, a fully nonlinear (full-f) five-dimensional (3d2v) gyrokinetic continuum edge-plasma code. (1) As a test of the interaction of collisions and parallel streaming, TEMPEST is compared with published analytic and numerical results for endloss of particles confined by combined electrostatic and magnetic wells. Good agreement is found over a wide range of collisionality, confining potential and mirror ratio, and the required velocity space resolution is modest. (2) In a large-aspect-ratio circular geometry, excellent agreement is found for a neoclassical equilibrium with parallel ion flow in the banana regime with zero temperature gradient and radial electric field. (3) The four-dimensional (2d2v) version of the code produces the first self-consistent simulation results of collisionless damping of geodesic acoustic modes and zonal flow (Rosenbluth–Hinton residual) with Boltzmann electrons using a full-f code. The electric field is also found to agree with the standard neoclassical expression for steep density and ion temperature gradients in the plateau regime. In divertor geometry, it is found that the endloss of particles and energy induces parallel flow stronger than the core neoclassical predictions in the SOL.

Reduced model simulations of the scrape-off-layer heat-flux width and comparison with experiment

Reduced model simulations of turbulence in the edge and scrape-off-layer (SOL) region of a spherical torus or tokamak plasma are employed to address the physics of the scrape-off-layer heat-flux width. The simulation model is an electrostatic two-dimensional fluid turbulence model, applied in the plane perpendicular to the magnetic field at the outboard midplane of the torus. The model contains curvature-driven-interchange modes, sheath losses, and both perpendicular turbulent diffusive and convective (blob) transport. These transport processes compete with classical parallel transport to set the SOL width. Midplane SOL profiles of density, temperature, and parallel heat flux are obtained from the simulation and compared with experimental results from the National Spherical Torus Experiment [S. M. Kaye et al., Phys. Plasmas 8, 1977 (2001)] to study the scaling of the heat-flux width with power and plasma current. It is concluded that midplane turbulence is the main contributor to the SOL heat-flux width f...

Analysis of geometric variations in high-power tokamak divertors

Quantitative assessment of the performance of high-power tokamak divertors for a range of geometric variations is conducted using the MHD code Corsica (Crotinger et al 1997 Technical Report LLNL) and edge transport code UEDGE (Rognlien et al 1992 J. Nucl. Mater. 196–198 347). In a multi-parametric study the divertor performance is compared for a high-power tokamak with standard and snowflake (Ryutov 2007 Phys. Plasmas 14 064502) configurations for the same core plasma parameters. Divertor and edge quantities that are varied include x-point flux expansion, shape of target plates, and radiating impurity species and concentrations. For a range of studied cases, in the snowflake the peak heat load on the target plates is significantly reduced compared with the standard divertor due to larger plasma-wetted area and a larger fraction of power radiated in the edge.

Effects of parallel electron dynamics on plasma blob transport

The 3D effects on sheath connected plasma blobs that result from parallel electron dynamics are studied by allowing for the variation of blob density and potential along the magnetic field line and using collisional Ohm’s law to model the parallel current density. The parallel current density from linear sheath theory, typically used in the 2D model, is implemented as parallel boundary conditions. This model includes electrostatic 3D effects, such as resistive drift waves and blob spinning, while retaining all of the fundamental 2D physics of sheath connected plasma blobs. If the growth time of unstable drift waves is comparable to the 2D advection time scale of the blob, then the blob’s density gradient will be depleted resulting in a much more diffusive blob with little radial motion. Furthermore, blob profiles that are initially varying along the field line drive the potential to a Boltzmann relation that spins the blob and thereby acts as an addition sink of the 2D potential. Basic dimensionless param...

Local properties of the magnetic field in a snowflake divertor

The power-law series for the poloidal magnetic flux function, up to the third-order terms, is presented for the case where two nulls of the poloidal magnetic field are separated by a small distance, as in a snowflake divertor. Distinct from the earlier results, no assumptions about the field symmetry are made. Conditions for the realization of an exact snowflake are expressed in terms of the coefficients of the power series. It is shown that, by a proper choice of the coordinate frame in the poloidal plane, one can obtain efficient similarity solutions for the separatrices and flux surfaces in the divertor region: the whole variety of flux surface shapes can be characterized by a single dimensionless parameter. Transition from a snowflake-minus to a snowflake-plus configuration in the case of no particular symmetry is described. The effect of the finite toroidal current density in the divertor region is assessed for the case of no particular symmetry.