Session NO2: Advances in Plasma Simulation II

Recent progress from Simulation of Wave Interactions with MHD (SWIM)

The SWIM project consists of three elements. 1) A computational platform (the Integrated Plasma Simulator or IPS), that allows efficient coupling of the full range of required fusion modules, that is flexible, permits convenient user access and access to experimental data, and that is robust to evolving physics, code development, and computer hardware. 2) A physics campaign addressing long timescale discharge evolution in the presence of sporadic fast MHD events. This involves interfacing the IPS to the 3D non-linear extended MHD codes and carrying out a program of research on use of RF to study and control fast time-scale MHD phenomena. The primary focus being the understanding of how RF can be employed to optimize discharges and control sawtooth events. 3) A physics campaign for modeling the direct interaction of RF and extended MHD for slowly growing modes, the primary physics focus being to understand how RF can be used to control neoclassical tearing modes. The IPS will contain modules to calculate wave propagation and absorption in all relevant frequency regimes, modules to calculate the modification of the plasma velocity distribution from sources (RF, neutral injection and alpha particles), calculate profile and magnetic evolution, linear MHD stability models and reduced models of non-linear MHD events. The IPS design will be described and preliminary results presented.   

Solution of the Helical Grad-Shafranov Equation for Magnetic Fields with Islands

We have developed a new solver for the helical Grad-Shafranov equation that can handle magnetic islands. Three applications are being pursued: 1) to serve as a testbed for new algorithms to be used in a fully 3D equilibrium code; 2) for benchmarking the PIES 3D equilibrium code for helical equilibria that have magnetic islands; 3) for testing and verification of a new capability being incorporated in the PIES code to handle neoclassical effects on magnetic islands. A Jacobian-Free Newton-Krylov method, including a linesearch algorithm and physics-based preconditioning, is currently being tested as a potential method for speeding the calculation of 3D equilibria with magnetic islands.   

Axisymmetric MHD Simulations of Pellet Ablation

A new 2D magnetohydrodynamic simulation of pellet ablation in the electrostatic approximation was developed. The main features of the model are the explicit tracking of solid pellet/ablated gas, and ablated gas/ambient plasma moving interfaces. Atomic processes in the ablation cloud are included. The major conclusion of the study is that in purely hydrodynamic simulations (without JxB force), changing the heat flux deposition from 1D spherically symmetric to 2D axisymmetric leads to a minor reduction in the ablation rate, contrary to the prevailing expectation of a ``factor of 2'' reduction. However, in the magnetohydrodynamic simulations (with JxB force), the magnetic field channels the flow into an extended plasma shield, and significantly reduces the ablation rate by a factor of 2 to 3, depending on the time it takes for the heat flux to ramp up as seen by a moving pellet. Fast pellets crossing pedestal regions in ITER would lead to shorter warm-up times, which in turn lead to narrower ablation channels, stronger shielding, and reduced ablation rates.   

Full Orbit PIC in NIMROD

The primary goal of the Plasma Science and Innovation Center (PSI Center) is to refine and optimize existing MHD codes to achieve improved predictability for emerging concept (EC) experiments. Kinetic effects have been shown to play a dominant role in some EC experiments, particularly in FRC stability[1]. The PSI Center has extended the hybrid kinetic-MHD implementation in NIMROD[2] from the drift kinetic model to the full kinetic model to include sufficient physics to accurately account for these effects, in particular large Larmour radius. The Boris push has been implemented for particles in NIMROD. However, this places a severe time step restriction on the particle time step. Time step restrictions have been decreased by using orbit averaging. Orbit averaging has been implemented for both full and drift kinetic equations and shows no significant impact on linear growth rates for tested regimes. We will show preliminary results from the implementation of full orbit (Lorentz force) particles coupled to the NIMROD code. \newline \newline [1] E. Belova et.al. ``Numerical Study of tilt stability of prolate field-reversed configurations,'' PoP, 7, 4996, 2000 \newline [2] C.C. Kim et.al. ``Hybrid Kinetic-MHD Simulations in General Geometry,'' CPC, 164, 448, 2004.   

Effects of Ion Landau Damping and Vortex Generation on the Kinetic Internal Kink Model

As the MHD phenomena in tokamaks enter into the collisionless regime, needs for the kinetic (or extended) MHD simulation have been increased. The simulation model in this article includes electron inertia, diamagnetic effect, and ion Landau damping. To simulate the m=1 and n=1 internal kink mode, the collisionless skin depth (d$_{e}\sim $1mm) which is three orders of magnitude smaller than the minor radius is resolved. The results are summarized as follows: (1) As the pressure gradient increases the internal kink mode (IK mode) is stabilized by the diamagnetic effects but the new mode appears with the smaller growth rate and with the frequency of drift wave (DIK mode). There is a stability window between unstable regions of both modes if the ion Landau damping is included. (2) DIK mode can also cause full magnetic reconnection. The linear mode pattern of the DIK mode is basically similar to that of the IK mode but has a strong poloidal shear flow to generate small vortices in the nonlinear stage due to the Kelvin-Helmholz like instability. There appears turbulent region due to the nonlinear coupling of vortices. The DIK mode is first stabilized by the excitation of vortices but later it is destabilized with enhanced growth rate.   

New time-adaptive technique for stiff gas dynamics and MHD

We present a new time-accurate algorithm for explicit integration of multi-scale gas dynamics and MHD equations. Unlike conventional time-stepping schemes, this technique is based on the discrete-event simulation (DES) methodology for stiff PDEs, which allows CPU resource adaptation in accordance with local time scales [1]. DES enables asynchronous (free of the global CFL restriction) flux-conservative integration on arbitrary spatial meshes and eliminates updates in inactive regions of the computation domain. We introduce a new preemptive event processing (PEP) approach, which automatically reduces event-driven integration to synchronous time stepping in regions of configuration space where the change of rate of the solution becomes uniform. This allows efficient parallelization of DES codes. We show that DES can be naturally applied to systems of coupled nonlinear equations. Importantly, we increase the temporal order of accuracy of DES by extending the previous algorithm [1] to second order in time. This is achieved by applying local corrections to the solution obtained with the forward Euler scheme. We demonstrate the accuracy, efficiency and robustness of DES-PEP by comparing numerical solutions obtained in event-driven and equivalent time-stepping simulations of several gas dynamics problems. [1] Y.A. Omelchenko and H. Karimabadi, \textit{J. Comp. Phys.} \textbf{216}, 179-194 (2006).   

Grid-Free Particle Simulations for Electrostatic Plasmas

A novel grid-free particle method for electrostatic plasma simulations is presented. The method employs numerical techniques from vortex methods in computational fluid dynamics including: (1) a multipole treecode algorithm to reduce the cost of evaluating the electric field induced by a set of point charges, (2) kernel smoothing to handle the singularity in the Coulomb potential, and (3) adaptive particle insertion to maintain resolution of the charge distribution in phase space using the Lagrangian flow map. The proposed method is intended as an alternative to particle-in-cell (PIC) methods. Simulations are presented for the instability of cold and warm collisionless electron beams in a neutralizing background of fixed ions.   

On a Primal Coarse Projective Integration Method for Multiscale Simulations

A novel simulation framework called Equation-Free Projective Integration (EFPI) was recently applied to nonlinear plasmas by M. Shay [1] to study propagation and steepening of a 1D ion sound (IS) with a PIC code as a microscopic simulator. To initialize, macro plasma variables are ``lifted'' to a fine micro-representation. PIC code is stepped forward for a short time, and the results are ``restricted'' or smoothed back to macro space. By extrapolation, time derivative is estimated and projected with a large step; the process is repeated. As a simple alternative, we propose a sort of a primal EPFI scheme to simulate nonlinear plasmas including kinetic effects. The micro-simulator is a standard 1D ES PIC code. Ions are assumed inherently coarse grained or ``smoothed'' and tracked to extrapolate in time and project. The potential is averaged over the electron plasma period to extrapolate and project. No adiabatic approximation for electrons is used [2], instead, self-consistently find the non-uniform electron distribution from the Poisson equation and ion density. Preliminary results for nonlinear IS as well as for the IS double layer paradigm are presented and some limitations on the EPFI discussed. [1] M. Shay, J. Drake, W. Dorland, J. of Comp. Phys (APS DPP 2005) [2] G. Stanchev, A. Maluckov et al., in EPS Fusion (Rome, 2006).   

Is PIC-MCC the right tool for the job?

Recent work by M. Turner (Phys. of Plasmas Vol. 13, 2006) point out a potential problem when combining PIC with Monte Carlo Collision (MCC) operators for modeling collisional plasmas. The outcome of Turner's work is that reaction rates in plasma calculations may be three orders of magnitude too fast. This means that most PIC-MCC codes may be providing a poor representation of the true plasma chemistry. However, Turner's work was 1D-1V and it is not clear how the additional degrees of freedom may play a role in the relaxation rates. There are several possible hypotheses for the phenomenon observed by Turner. These include; issues of an incomplete phase space, an interplay between local PIC errors with the statistical noise of the MCC operator and a non-local interplay between PIC and MCC errors. In this work we systematically explore these hypotheses using 1D-3V PIC and Boundary Integral Treecodes (BIT) and propose several solutions.   

Particle dynamics and energization at the earth bow shock.

Simulation results are presented for the dynamics and energization of plasma particles as they cross the earth bow shock. The model is based on a kinetic description of particle dynamics, in which single particle trajectories are calculated using a symplectic integration scheme. Particle distribution functions are calculated at different positions in the magnetosheath using Louville's theorem with a Hamiltonian formulation of single particle equations of motion and time reversed particle tracing. As a first approximation, the bow shock is assumed to be represented by a prescribed local plane shock. A more detailed representation of the shock and magnetosheath regions is then considered, as obtained from a global MHD simulation model. Comparisons are made between the two models, for different possible angles between the incident plasma flow velocity and the shock front. Various scenarios of energization are also considered for majority (H) and minority (He and heavier elements) ion species.   

Viscous Heating of Ions through Saturated Fine-Scale MHD Instabilities in a Z-Pinch at 200-300 KeV Temperature

Pulsed power driven Z-pinches yield large X-ray powers at stagnation, the energy of which can exceed by up to factors of 3 or 4, the estimated kinetic energy of the implosion. Furthermore, when electron temperatures are measured at stagnation similar in temperatures would not lead to pressure balance. These problems can be resolved by a theoretical model in which short wavelength (ka $>>$ 1, and viscous Lundquist number $\sim $ 1), fast growing, m=O MHD instabilities reach a saturated amplitude, and the associated viscous dissipation of the vortices leads to ion heating. Equating this heating rate to the equipartition of energy to electrons leads to an estimate of the ion temperature and pinch radius at pressure balance. Extremely high ion temperatures in the range of 200-300 KeV are predicted from this model for stainless steel wire array experiments on Z at Sandia. These have been confirmed from time-resolved Doppler broadening spectroscopic measurements of the optically thin Fe He-$\delta $ line. This conversion of magnetic energy into ion thermal energy occurs on the nanosecond timescale, and can prevent radiative collapse. Any accompanying loss of magnetic flux in this highly conducting plasma can be explained by the occurrence of a large number of hot spots along the axis, with electron density and temperature variating not exactly in phase. This leads to a significant value of the integral of \underline {E.dl}. Dl along the axis due to the grad Pe term in Ohm's law, analogous to the magnetic field generating term found in laser-plasma interactions. Ref 1. M.G. Haines, et al; Phys. Rev. Lett. 96, 075003 (2006) Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy under Contract DE-ACO4-94AL85000.   

Evidence of Non-Classical Transport in Discharge Hollow Cathodes

Measurements, simplified analyses and 2-D numerical simulations with a fluid plasma model show that classical resistivity can not fully account for the elevated electron temperatures and steep plasma potential gradients measured in a discharge hollow cathode used as an electron source in electric propulsion. The numerical simulations show that classical resistivity yields much colder electron temperatures compared to the measured values in the orifice and near-plume regions of the cathode. Classical transport and Ohm's law also predict exceedingly high electron-ion relative drift speeds compared to the electron thermal speed, which would normally excite streaming instabilities in this plasma. Measurements with a high frequency emissive probe also capture 50-100 kHz plasma potential fluctuations in the near-plume region. It is found that the addition of anomalous resistivity based on existing growth rate formulae for such instabilities improves significantly the comparison between fluid theory and experiment. But the extent of possible deviations from the Maxwellian EEDF, possibly as a result of the micro instabilities, has not yet been quantified in this cathode.   

HDF5WS -- Web Service for Remote Access of Simulation Data

Data produced by modern plasma physics and fusion simulations is growing in size and typically resides on a remote site: a supercomputer or a cluster. In order to analyze and visualize the data, one needs to query and extract subsets of interest rather transfer the bulk of it. In this paper, we introduce our solution to this problem. It is a Web Service based on Globus Toolkit. The service client's API is written in C++ and has a set of methods commonly used to query and access HDF5 data, which is the most popular data format used in plasma physics and fusion simulations. Through this API, users can query attributes of remote datasets, extract particular datasets and hyberslabs into the client memory as if HDF5 files were local. The data transfer mechanism used in the service is gridFTP. In addition to describing the service, we provide multiple benchmarking results, comparing various data formats, types of middleware and data transfer mechanisms. These results determined the design of the service.   

Cross Sections for Charge Exchange and Ionization of High Energy Ions in Noble Gases

Interactions of charge exchange and ionization of fast, low-charged heavy ions are very important in heavy ion beam inertial confinement fusion (HIBF). These effects are crucial in determination of the final focusing in the chamber. However, corresponding cross section data is very limited and/or not accurate over the entire range of energies and ions of interest. This paper reports on our recent studies of cross sections for interactions of heavy ions like Cs+ with noble gas. Since a quantum mechanical treatments encounters with the many body problem, classical trajectory Monte Carlo method is employed. The distribution of inner electrons is estimated by modified Thomas-Fermi Model for the purpose of decreasing the number of electrons to calculate their orbits. Introduced micro-canonical ensemble for the initial electron probability distribution describes quantum mechanical uncertainty. Cross sections are evaluated over a limited energy range. Corresponding scaling laws are then developed to reflect the change probability of beam charge state over a larger energy range. Special attention is given to multi-electron processes, especially double-fold interactions.