M3D Code Research Articles

Saturated helical mode in EAST high β hybrid plasmas

A hybrid scenario with the minimum safety factor above unity and a safety factor profile with a wide range of low magnetic shear in the core is one of the candidates for ITER operations. Theory predicts that such equilibria are prone to quasi-interchange mode, which was confirmed by observations from several devices. The presence of those 3D quasi-interchange helical modes has both positive and negative effects on plasma performance. On the one hand, they can avoid low mode number instabilities such as sawteeth. On the other hand, they would lead to confinement degradation, significant fast ion losses and toroidal rotation damping. Such modes are observed for the first time in EAST high βp hybrid plasmas achieved with off-axis neutral beam injection (NBI) heating. The mode manifested itself as long-lived saturated helical instability. Notably, when the mode was destabilized, the toroidal rotation profile in the core was completely flattened. The extent of the flattened rotation profile reaches up to half of the minor radius. As the minimum safety factor qmin approaches unity, strong n = 2 and n = 3 harmonics are successively driven unstable as nonlinear consequences of the n = 1 mode. A 3D resistive nonlinear simulation with a realistic tokamak configuration has been applied to reproduce the evolution of long-lived helical instability with the M3D code. The simulated magnetic perturbation shows a strong n = 2 component. As the mode grows, the amplitude of the n = 2 harmonic grows immensely with respect to that of the n = 1 harmonic. The n = 2 harmonic gradually becomes dominant as the mode saturates. In saturation, the initial weak shear safety factor profile in the core becomes completely flat. This change of safety factor profile in the core corroborates a self-regulating magnetic flux pumping mechanism which is considered responsible for maintaining stationary non-sawtoothing hybrid discharges.

Linear and nonlinear simulations of the visco-resistive internal kink mode using the M3D code

Numerical simulations of visco-resistive internal kink modes with dominant poloidal and toroidal mode numbers m/n = 1 in tokamak plasma configurations have been carried out using the M3D code. The numerical scaling of the growth rate of linear visco-resistive internal kink is found to agree reasonably well with the analytical theory. The impact of toroidal effects is assessed. Nonlinear simulations show that when viscosity increases, the nonlinear evolution of resistive internal kinks transits from cyclic, Kadomtsev-like sawtooth oscillations at relatively low values of the viscosity parameter, to pulsating m/n = 1 magnetic islands at intermediate viscosity, to steady state islands at high viscosity. For the pulsating and steady state scenarios, the safety factor profile is nearly flat and almost equals to unity in the core region. Regimes with pulsating or steady state islands may be relevant to the interpretation of partial sawteeth and of the so-called snake phenomenon sometimes observed in tokamak experiments.

Mode coupling and aspect ratio effects on low and high-n plasma instabilities

In magnetically confined toroidal plasmas such as tokamaks, magnetohydrodynamic (MHD) instabilities experience strong toroidal and nonlinear mode coupling effects. Resistive MHD simulations with the M3D code show the importance of mode coupling and compressible MHD effects, which contribute to stronger mode coupling. For the m/n = 1/1 internal kink mode and sawtooth crash and for the edge localized mode (ELM) at higher n, MHD reproduces many features of the experimental observations, including the fast sawtooth crash and the moderate n ∼ 10 toroidal harmonics of the ELM. A general property of the perpendicular momentum equation in toroidal fusion plasmas is that the unbalanced radial forces remain relatively small, so that the terms that are lowest order in small inverse aspect ratio mostly cancel. The higher order terms then have significant effects, even at small r/Ro and small amplitude. Effects are strongest for the lowest toroidal harmonics n ≃ 1 and the most strongly driven ones with highest amplitude. Unlike the n = 1 internal kink mode, the small amplitude ELM ballooning/peeling-type mode, and thus ELM MHD marginal stability, may be reasonably described by the lowest order in aspect ratio, for moderate and large n ≳ 10. The ELM crash, however, depends on higher order.

Wall-touching kink mode calculations with the M3D code

This paper seeks to address a controversy regarding the applicability of the 3D nonlinear extended MHD code M3D [W. Park et al., Phys. Plasmas 6, 1796 (1999)] and similar codes to calculations of the electromagnetic interaction of a disrupting tokamak plasma with the surrounding vessel structures. M3D is applied to a simple test problem involving an external kink mode in an ideal cylindrical plasma, used also by the Disruption Simulation Code (DSC) as a model case for illustrating the nature of transient vessel currents during a major disruption. While comparison of the results with those of the DSC is complicated by effects arising from the higher dimensionality and complexity of M3D, we verify that M3D is capable of reproducing both the correct saturation behavior of the free boundary kink and the “Hiro” currents arising when the kink interacts with a conducting tile surface interior to the ideal wall.

Influence of resistive internal kink on runaway current profile

An extended magnetohydrodynamic model including the effect of runaway current is developed and implemented in the M3D code. Based on this model, a simulation has been carried out to investigate the linear stability and the nonlinear evolution of the n = 1 resistive kink mode in the presence of runaway current. It is found that sawteeth oscillation is suppressed in a runaway plasma. The nonlinear evolution of the n = 1 mode only leads to a single sawtooth crash before reaching a new steady state axi-symmetric equilibrium with flattened current profile in the plasma core and q(0)>1.

Tokamak toroidal rotation caused by AVDEs and ELMs

Toroidal rotation can be produced by disruptions, as observed in several experiments. There is a concern that rotating asymmetric forces during an ITER disruption might resonate with the blanket and other structures surrounding the plasma. Here it is shown, both computationally using the M3D code, and analytically, that toroidal rotation is produced by magnetohydrodynamic turbulence. In particular, rotation is produced during an asymmetric vertical displacement event (AVDE) disruption. Toroidal and poloidal rotation are also produced during edge localized modes (ELMs), and may be consistent with a scaling law found for intrinsic toroidal rotation in H-mode tokamaks.

Stagnation of a gas puff Z pinch

Simulations of a gas puff Z pinch were performed, using an appropriately modified version of the M3D code [Park et al., Phys. Plasmas 6, 1796 (1999)]. The simulations investigated the stagnation process, including the effects of the shock driven by the magnetic piston and the influence of the Rayleigh Taylor instability. The results compare favorably with recent experimental measurements. The stagnated plasma reaches a quasi equilibrium with approximate balance of plasma pressure and magnetic pressure, measured by β≤1. The dependence of the stagnation radius on β and radiative energy loss are calculated, using a simple radiation model.

Free-boundary magnetohydrodynamic equilibria with flow

The finite-element M3D code [W. Park et al., Phys. Plasmas 6, 1796 (1999)] has been modified to include a free-boundary equilibrium solver with arbitrary toroidal and poloidal flows. With this modification, the M3D code now has the capability to self-consistently model two essential ingredients necessary for equilibrium calculations in the edge region, namely, free-boundary and arbitrary flow. As a free-boundary code, M3D includes the separatrix and scrape-off layer regions in the equilibrium calculation. Poloidal flows in the subsonic, supersonic, and transonic regimes can be calculated with the new version of the M3D code. Calculation results show that the presence of equilibrium flows, in particular those next to the plasma boundary, can considerably influence the position of the X-point and magnetic separatrix shape/location and hence the position of the strike point on the divertor plates. Moreover, it is shown that poloidal flow is not a rigid-body rotation, with the fastest flows occurring on the inboard side of the plasma. A numerical confirmation of the “de Laval nozzle” model of Betti and Freidberg [R. Betti and J. P. Freidberg, Phys. Plasmas 7, 2439 (2000)] for free-boundary equilibrium calculations is obtained, with the formation of the predicted discontinuities between regions of subsonic and supersonic flows (with respect to the poloidal sound speed). Finally, a detailed comparison between isentropic and isothermal equilibria is presented, showing qualitative analogies and quantitative differences.

Wall forces produced during ITER disruptions

Nonlinear simulations with the M3D code [W. Park et al., Phys. Plasmas 6, 1796 (1999)] are performed of disruptions produced by large scale magnetohydrodynamic instabilities. The toroidally asymmetric wall forces produced during a disruption are calculated in an ITER [T. Hender et al., Nucl. Fusion 47, S128 (2007)] model. The disruption is produced by a vertical displacement event and a kink mode. Expressions are derived for the wall force, including the sideways force, using a thin conducting wall model. The scaling of wall force with γτw is obtained, where γ is the kink growth rate and τw is the wall penetration time. The largest force occurs with γτw≈1. A theory is developed of the wall force produced by kink modes. The theory is in qualitative agreement with the simulations and Joint European Torus [V. Riccardo et al., Nucl. Fusion 49, 055012 (2009)] experiments. In particular, the theory and simulations give dependence of the sideways on γτw, correlation of sideways force with sideways plasma displacement, and correlation of toroidally varying plasma current with toroidally varying vertical displacement.

Implementation of an implicit shear Alfvén operator in the M3D code

A new pair of linear operators is introduced into the partially implicit but largely explicit time advance algorithm for the MHD equations in the M3D code. These operators jointly render the algorithm implicit with respect to propagation of the shear Alfvén wave to leading order. As a result, the maximum stable time step is increased by up to a factor of three for many practical applications. The additional computational cost is trivial for the linear mode of operation and adds only approximately 30% additional execution time per step for nonlinear calculations.

Massively parallel modeling of the sawtooth instability in tokamaks

The sawtooth instability is one of the most fundamental dynamics of an inductive tokamak discharge such as will occur in ITER. The sawtooth occurs when the current peaks in a tokamak, creating a region in the core where the safety factor is less than unity, q<1. While this 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 or set off a sequence of events that leads to a major disruption. Sawtooth behavior is complex and remains incompletely explained. The SciDAC Center for Extended MHD Modeling (CEMM) has 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. The multiple time and space scales associated with the reconnection layer and growth time make this an extremely challenging computational problem. The most recent simulation by the M3D code used over 500,000 elements for 400,000 partially implicit time steps for a total of 2×1011 space-time points, 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.