MAGNETOHYDRODYNAMICS AND STABILITY RESEARCH

The GA Theory group has an extensive history of pioneering contributions to the study of magnetohydrodynamics (MHD), including significant advances in analytic theory, large-scale numerical simulation, and quantitative study of actual experiments. These studies have covered a broad range of topics, focusing particularly on the onset conditions and dynamics of important tokamak instabilities, including global beta and current limiting instabilities, edge localized modes (ELMs), resistive wall modes (RWMs), and tearing modes. A variety of theoretical frameworks have been developed and employed, from static ideal MHD to more comprehensive resistive and extended formulations incorporating rotation, resistivity, viscosity, fast ions, and other drift and kinetic effects. A major characteristic of the GA MHD theory work is that the basic theory is informed by issues of numerical tractability, and both the basic theory and numerical tool development are well grounded by continuous testing on applications.

Code development is an important implicit element of our research. The numerous MHD codes under development at GA are actively used at institutions across the world, including IFS, U. of Wisconsin, MIT, Columbia U., PPPL, U. of Washington, IPP Garching, Culham/JET, Korea Basic Science Institute, and the National Institute for Fusion Science (NIFS Japan), and the ASIPP (China). Also, it should be noted that much of the research work has implications for non-tokamak magnetically confined plasmas, such as spherical tori, reversed field pinches (RFP), spheromaks and field reversed configurations.

The following provides a summary of the recent research interests in MHD stability of the GA Theory Group. These range from global ideal MHD stability studies, studies of edge stability and its relation to the phenomena of Edge Localized Modes (ELMs), stability in the presence of a resistive wall, and MHD stability in a resistive plasma. The descriptions of work in each of these areas are intended to provide some motivation and background as well as a summary of some of the more important contributions of the group. References are provided at the end.

Global Ideal MHD: Although the theory of ideal MHD equilibrium and stability is mature and well established, significant advances have continued over the past grant period. Notably, detailed studies using the GATO code [Bernard 1981], and employing synthetic diagnostics, have extensively demonstrated the predictive capability of ideal MHD, including accurate predictions of b limits, growth rates, and mode structures [Turnbull 2002]. Significant theoretical issues, such as the apparent phase reversal of magnetic fluctuation signals on the inboard side of the tokamak, have been resolved, and the correspondence between exceeding the no-wall beta limit and observed rotation slowdown has been clearly established. In addition, studies of fundamental limits on MHD equilibrium have led to the identification of a simple formula for the limit on equilibrium bp [Lin-Liu 2003], and a demonstration that equilibria with ebp as high as 5 are possible, a result with clear positive implications for the prospects of high bootstrap fraction tokamak scenarios. Current hole equilibria have also been studied, with a proof that zero current equilibria do exist, while robust, axisymmetric negative current equilibria do not [Chu 2002b].

As ideal MHD studies have become increasingly routine and well understood, the Theory group's focus is continuing to move toward resistive and extended MHD physics, and the onset and dynamics of modes which impact plasma confinement and performance below global ideal wall stability limits.

ELMs and the Pedestal: Gaining an understanding of ELMs, including onset conditions and dynamic evolution, has come to the fore as a critical issue for ITER and other burning plasma experiments, both because of the potential impact of ELM pulses on material surfaces, and because ELM onset places an effective constraint on the pressure at the top of the edge barrier (the "pedestal height"), which strongly impacts core confinement and overall fusion performance. The GA Theory group has made important breakthroughs in physics understanding of ELMs over the past grant period. In particular, the peeling-ballooning model of ELMs, pioneered by GA in collaboration with Culham, has been quantified, elaborated, and extensively tested against experiment in a series of talks and papers that has led to broad community acceptance of the theory. The peeling-ballooning model posits that intermediate wavelength (typically 3<n<40) MHD modes, driven by a combination of the strong pressure gradient and resulting bootstrap current in the edge barrier region, are responsible for ELMs and constraints on the pedestal height [Connor 1998, Snyder 2002]. The bootstrap current plays a complex dual role, on the one hand driving peeling/kink modes, while on the other lowering magnetic shear and opening second stability access to high n modes in shaped plasmas. The collisionality dependence of the bootstrap current leads to separate dependencies of the stability limit on temperature and density [Snyder 2003]. Second stability effects for finite n modes in shaped toroidal equilibria have been demonstrated and quantified, allowing determination of the most unstable wavelength [Snyder 2002].

Due to the complex dependencies of the stability limit on equilibrium characteristics, and the wide range of relevant wavelengths, a highly efficient numerical tool is needed for comprehensive quantification of peeling-ballooning stability bounds and comparison to experiment. To this end, the ELITE code [Wilson 2002, Snyder 2002], which is based upon an extension of ballooning theory to include higher order terms and appropriate external boundary conditions, was developed and successfully benchmarked. Pedestal stability studies on multiple tokamaks using the intermediate to high n ELITE code, together with low n codes such at GATO and DCON, have consistently found that intermediate wavelength modes become unstable shortly before ELMs are observed, and that several ELM characteristics, such as localization in the outer radial region of the low field side, are consistent with predictions [e.g. Snyder 2002, Lao 2001, Snyder 2004a, Mossessian 2003, Petrie 2003, Turnbull 2003]. In addition, a technique has been developed which systematically employs ELITE calculations to characterize peeling-ballooning limits as a function of key equilibrium parameters. This has allowed comparisons of peeling-ballooning stability bounds to the experimental pedestal database, finding good agreement in the pedestal height (at a given pedestal width) just before an ELM, as well as in its trends with current (see Figure below), triangularity, and density [Snyder 2004a]. Projections of ITER pedestal constraints [Snyder 2003], and systematic scaling studies of peeling-ballooning stability bounds [Snyder 2004b] have also been conducted.

Calculated peeling-ballooning stability threshold (line), as a function of total plasma current, agrees with observed pedestal pressure (circles) shortly before ELMs occur (left figure). Characteristic unstable peeling-ballooning mode structure calculated with ELITE just before an ELM in DIII-D, without (central figure), and with strong toroidal rotation (right figure).

Recently, the dominant compression and toroidal rotation terms have been incorporated into ELITE, revealing that differential rotation shears apart and radially narrows the peeling-ballooning mode structure [See figure below], while having only a relatively small effect on growth rates away from marginal stability [Snyder 2003b]. Related studies have resolved the apparent paradoxical discontinuity in growth rates between the zero rotation case [in which one selects the ballooning angle q0 which maximizes g(q0)], and the small rotation limit of the finite rotation case [in which one averages g(q0) over all q0] [Webster 2003].

The successful development of the ELITE code, and substantial progress of the peeling-ballooning model have laid a foundation of basic understanding of several key ELM physics issues. However, a quantitative study of ELM size, as well as the important issue of transport of heat and particle pulses across the open field line region and into material surfaces, requires moving beyond linear studies to 3D nonlinear, dynamical simulation and supporting theory. Dynamical simulation of the pedestal and boundary region is an enormous challenge, due to the wide range of overlapping spatio-temporal scales associated with electron and ion drifts, Alfvn waves, collisions, and neutral physics, as well as the breakdown in scale separation between perturbed and equilibrium timescales, particularly during the ELM crash. The sharp gradients in the pedestal region generally preclude a significant scale separation between the source and transport physics expected to determine profile shapes, and the stability physics associated with peeling-ballooning modes; both are needed to fully determine constraints on the pedestal height.

Nevertheless, significant initial progress in studies of the evolution of edge instabilities in the presence of non-ideal effects has been made. The NIMROD [Sovinec 2003] code has recently incorporated a model of the open field line region as a high resistivity plasma, allowing study of external modes. In initial studies, edge instabilities with n as high as 20 have been identified in the linear phase, with mode structures and growth rates in the range expected from linear ELITE and GATO calculations. Initial nonlinear studies are underway. In collaboration with LLNL, we have, undertaken nonlinear simulations using an enhanced version of the BOUT [Xu 2003] code, modified to include the kink term, allowing treatment of peeling-ballooning modes in the presence of diamagnetic and resistive effects. Several simulations of linearly unstable equilibria have been undertaken, revealing the expected ballooning-like mode structure and growth rates at early times, followed by poloidal narrowing, and bursts of fast radially propagating filaments [Snyder 2003b], similar to those observed during ELM crashes. An example is shown in the figure below.

3D BOUT simulation of the dynamic evolution of an edge instability: On the left is the time evolution of the perturbed density near the outer midplane, showing the mode growing in the sharp gradient pedestal region on the expected peeling-ballooning timescale (early times), and then (t > ~2000) bursting radially both inward and outward into the SOL. The perturbed density along a bundle of field lines at t = 2105, projected onto an RZ plane is in the central figure (note that the field line bundle extends into and out of the page). The rightmost figure shows a closeup view of the outer midplane region shows negative (dark) density perturbations near the top of the pedestal and a burst of positive (light) density propagating into the SOL.

Resistive Wall Modes: In the next generation of steady-state advanced tokamaks, a high-performance plasma must remain stable over time scales long compared to the flux diffusion time tw of the external resistive wall. It is also desirable to have bN = baB/I as high as possible for optimal performance, if possible exceeding the maximum bN NW stable to global modes in the absence of a conducting wall. This can be achieved by placing a conducting wall close to the plasma. However, when bN exceeds bN NW, the perturbation flux due to the unstable kink mode can diffuse through the external conducting wall on a time scale longer than tw, resulting in destabilization of the resistive wall mode (RWM). To fully realize the advantage of the conducting wall, the RWM must be stabilized. The goal of the research on RWM is to find methods for its stabilization and to develop a quantitative model to predict its stability in present and future advanced tokamaks.

Using the MARS code, Bondeson and Ward [Bondeson 1994] predicted that rotating the plasma relative to an external resistive wall can compensate for the inherent resistivity of the external wall. The plasma rotation needed for the stabilization of the RWM depends critically on the dissipation model for the plasma angular momentum. The kinetic damping model presently implemented in MARS is a large aspect ratio approximation to the usual analytic model. This model neglects the diamagnetic and magnetic curvature drifts, and is generally valid for plasma with large rotation, but is violated in experiments with low rotation.

Utilizing plasma rotation to stabilize the RWM was demonstrated in DIII-D [Strait 1995] and modeled [Chu 1995] using a modified version of the MARS code that includes differential plasma rotation. It was concluded that the plasma rotation needed for the stabilization of the RWM depends critically on the angular momentum dissipation model, and that the required rotation may prove difficult for future reactors to achieve. An alternative approach is magnetic feedback, which makes the resistive wall appear ideally conductive to the plasma. Both rotational stabilization with complete stabilization of the RWM [Garofalo 2002, Strait 2003], and direct feedback control [Okayabashi 2001, Strait 2004] have been demonstrated on DIII-D.

The GA Theory group has made substantial progress in RWM studies of both feedback and rotational stabilization. A nearly complete theoretical model for the feedback stabilization of a static plasma with a general equilibrium configuration has been developed, utilizing a formulation of the stability problem in terms of normal modes of the plasma-resistive wall system. The growth or damping rates and the eigenfunctions of the normal modes are determined via an extended energy principle for the plasma during its open (feedback) loop operation, using a set of equations derived for the time evolution of these normal modes with currents in the feedback coils. The feasibility of a given feedback system is evaluated via the Nyquist diagram method or by solving the characteristic equations. The elements of the characteristic equations are formed from the growth and damping rates of the normal modes, the sensor matrix of the perturbation fluxes detected by the external feedback coils, and the feedback logic. This method has been implemented numerically by coupling the DCON [Glasser 1997] code with an extended version of the VACUUM code [Chance 1997] and a feedback code. The resultant code package (called NMA for normal mode approach), has been applied to the DIII-D tokamak, where it is found that feedback with poloidal sensors is much more effective than feedback with radial sensors. [Chu2003b]

The effectiveness of a differentially rotating resistive wall in stabilizing the RWM is also being investigated. For a non-circular tokamak, wide ranges of flow patterns are all effective [Chu 2003a]. A new general equation was derived for the RWM in a rotating plasma in the presence of an external resonant field [Chu 2003a]. Analysis showed that the RWM growth and damping are determined by the balance of energy flux and dissipation torques between the plasma and the resistive wall. In steady state, only the plasma dissipation contributes to the phase shift between the external resonant field and the plasma response. The phase and amplitude of the plasma resonant response are related to the damping rate and mode frequency of the RWM through a complex amplification factor. These conclusions provide an explanation of the observations of recent MHD spectroscopy experiments on DIII-D [Reimerdes 2004]. We note that our approach is complementary to the approach adopted by the VALEN code [Bialek 2001]. The present VALEN approach includes all the 3D external coil geometry, which is important for engineering considerations, but adopts a simple plasma model and does not include plasma rotation. The approach the GA group has taken, treats the detailed plasma dynamics, necessary for a comprehensive understanding of plasma rotation damping, yet adopts two dimensional idealizations of the geometry of the external coils.

Rotational stabilization is pres­ent in real feedback experiments as well, and the effects of rotation and feedback are synergistic. For plasma with rotation, the linear non-ideal MHD code, MARS-F [Liu 2000], has been found to be a suitable tool. New features of the MARS-F code are the addition of feedback and a new kinetic model. This code has been benchmarked in the ideal MHD limit against NMA, and both codes have been tested for their prediction of the feasibility of feed­ feed­back stabilization. Two dissipation models, the ion sound wave damping model and the kinetic damping model based on the kinetic energy principle, have been employed in the MARS-F code to study the predicted critical rotation speed for rotation stabilization of the RWM. Broad qualitative features of the experiment are reproduced. When scaled up to ITER-FEAT parameters, 33MW of 1MeV negative neutral beam injection is projected to be able to sustain plasma rotation to stabilize the RWM without relying on feedback. MARS-F has also been utilized to study the synergistic effect of feedback and rotation on stabilization of the RWM in a rotating plasma with feedback coil located outside of the plasma chamber, as shown in the figure below [Chu 2004].

Synergistic effects of rotation and feedback using MARS-F with poloidal and radial feedback sensors. Increasing feedback gain is effective up to gains ~1-2, after which effectiveness is strongly reduced. Rotation must first reduce the growth rate below a low level before feedback gains in this range are sufficient to stabilize the model.

Resistive and Extended MHD: The field of linear resistive MHD is relatively mature, yet much basic theory development remains to be done both in linear and nonlinear resistive and extended MHD, and the Theory group is continuing to search for new and innovative approaches. One concept being pursued is the development of new stability criteria based on the Poynting vector, which quantifies net energy flows into and out of the system. In parallel, we are also investigating the physics of the inertial inner layer. This is intended to provide a rigorous physical grounding for using the outer region matching data as a linear stability criterion for general finite b 2D systems; the matching data has only been rigorously and unambiguously shown to be identical to the stability criterion in 1D zero b systems. A new innovative theory of Almost Ideal MHD (AIMHD) [Jensen 2001], based on the concept that resistive MHD relaxes certain constraints but maintains others from the ideal case, is also under development. This essentially nonlinear theory is yielding new insights into the nature of the problem by comparing theoretical ideas about appropriate conservation properties with full resistive simulations. Integral equilibrium constraints are imposed on the helicity contained within regions of magnetic flux space, which determine both the equilibrium and the nonlinear saturated state of the unstable system simultaneously, bypassing questions of actual time dynamics. The distribution of these integral constraints is critical to the prediction of the saturated state.

Linear tearing mode theory has a long history at General Atomics from the early work of Chu, Greene, and Jensen [Jensen 1983, Chu 1989, Chu 1993]. The linear asymptotic matching theory developed in [Chu 1993], provided a new fundamental approach to understanding the linear problem in finite b, 2D systems, and its relation to ideal stability. A crucial innovative development in finding a numerically tractable formulation was subsequently provided [Galkin 2000]. This formulation was incorporated into the 2D TWIST-R code. The TWIST-R code is now working and provides physically meaningful solutions over an extended range of parameters. An example solution is shown in the figure below. This shows the renormalized total displacement x = |x|(1/2+m) x. The same algorithm was incorporated into the inner layer code TWIST-IR, and it was demonstrated that the results match old published results [Galkin 2002]. The considerably extended range and flexibility were also demonstrated and the code applied to consider a non-uniform density, yielding new insights into the nature of the matching and the role of inertia [Galkin 2002, Greene 1995].

Linear resistive MHD solution obtained from TWI ST-R code for a Solovev equilibrium with m = 0.82. The left figure shows contours of the renormalized total displacement x* = lxl(1/2+m) x. On the right are profiles of the resonant m = 2 normal and
tangential components of the normalized displacement and the q profile.

A highly successful phenomenological approach based on the Modified Rutherford Equation (MRE) model [e.g. Sauter 1997] is also being pursued. The crucial contribution from the GA Theory group in this area has been to identify the critical role of the time and b dependent linear driving term D(b(t),t) [Brennan 2002, Chu 2002]. The basic model, in which evolving the linear stability in addition to modeling the nonlinear drive in the neoclassical island evolution are both crucial, has been verified in a few significant cases in DIII-D. Experiments on DIII-D have confirmed predictions from the model both at low b [Chu 2002] and higher b, with a rapid increase in the linear drive on approach to an ideal limit [Brennan 2002, Brennan 2003]. For sawtooth seeded NTMs, we have found very good agreement between the D(b(t),t) computed by PEST-III, and the values estimated from the MRE using measured island decay rates in DIII-D [Turnbull 2002, Brennan 2003] as shown in the leftmost figure below. Note that a key ingredient in the MRE modeling is an accurate finite b calculation for the linear drive. This was obtained from the PEST-III code using equilibria reconstructed from the discharge equilibria. This work has provided a fairly comprehensive description of both seeded and spontaneous NTM phenomena in DIII-D, and now appears to explain sawtooth seeding of 2/1 modes in Alcator C-Mod as well. An interesting feature is the prediction of a lull in the growth rate during the NTM growth under certain conditions, and has now been observed in the experiments and also subsequently found to be common but unnoticed, in many previous discharges [Brennan 2003].

On the left is the theoretical and observed 3/2 island decay and eventual onset after successive seeding from sawteeth. The final island growth occurs when D(b(t),t) computed from the reconstructed equilibrium at each time becomes significantly positive, driving the seeded mode unstable, in agreement with the MRE prediction. The right figure shows the perturbed n=1 electron temperature in a nonlinear simulation of a 1/1 mode seeding various seed islands in DIII-D in the presence of toroidal rotation. During the linear phase, the standard linear eigenfunction is hardly affected by the rotation, while during the nonlinear phase shown here, the eigenfunction appears distorted by the rotation, and yet still drives visible 2/1 seed islands (and 3/2 islands, not shown, in n=2). [See HOME PAGE for 3D visualization]

To test and further advance dynamical aspects of the tearing mode theories described above, large scale nonlinear simulations are being undertaken, employing the NIMROD extended MHD code, as part of a collaboration with the NIMROD team. Recent developments have resulted in inclusion of a finite vacuum region and single fluid sheared plasma rotation, both required for the realization of a predictive numerical tool, plus a fluid closure incorporating neoclassical effects. Extensive calculations using NIMROD have been utilized in support of the MRE predictions. Notably, simulations have reproduced the MRE stability predictions in the presence of nonlinear coupling, and have found super-exponential growth of spontaneous tearing modes with increasing linear instability, lending further support to our MRE based onset theory. Several cases of tearing mode onset driven by nonlinear coupling, in the presence of rotational shear, are also being studied, to determine the effects of rotation on the driven flux. Calculations for DIII-D sawtoothing equilibria showed that the m/n = 2/2 component of the sawtooth does in fact provide the seed for a subsequent 3/2 island. The figure on the right above shows a simulation of the nonlinear evolution of a 1/1 mode seeding islands in DIII-D in the presence of toroidal rotational shear. In the near future, the inclusion of a resistive wall boundary condition will allow the study of the nonlinear interaction of unstable modes such as these with resistive wall modes and imposed error fields.