3D nonlinear modeling of the coupling and phase locking of magnetic Islands in tokamaks

Many tokamak discharges develop multiple tearing modes possessing different mode numbers. These modes are observed to phase lock to one another, resulting in a flattening of the core toroidal plasma rotation profile, which can have deleterious effects on transport and MHD stability. In order to study these phenomena with minimum assumptions, we use the M3D-C1 3D nonlinear MHD code to perform initial value simulations of the evolution of equilibria unstable to both the 2/1 and 3/2 modes, but having sheared toroidal rotation. Initial attempts to perform these simulations led to numerical instabilities developing once the islands got to a certain size. In order to study the cause of this instability, we developed a small model code that solves a pure convection equation in 1D. We find that an implicit Crank-Nicholson method in time and Hermite Cubic finite elements (as are used in the toroidal direction in the M3D-C1 code) is not a convergent algorithm. Adding a small second order diffusion term, proportional to the velocity, improves the numerical stability properties but is not convergent in the first-derivative of the solution. Instead, adding a much smaller forth-order spatial derivative term proportional to the velocity leads to an algorithm in which both the solution and the first derivative converge as 1/N$^{\mathrm{2,}}$. Adding similar toroidal forth derivative terms to the M3D-C1 code eliminated the numerical instability. This work was supported by US DOE Contract DE-AC02-09-CH11466.