J. A. Breslau

A high-order implicit finite element method for integrating the two-fluid magnetohydrodynamic equations in two dimensions

We describe a new method for solving the time-dependent two-fluid magnetohydrodynamic (2F-MHD) equations in two dimensions that has significant advantages over other methods. The stream-function/potential representation of the velocity and magnetic field vectors, while fully general, allows accurate description of nearly incompressible fluid motions and manifestly satisfies the divergence condition on the magnetic field. Through analytic manipulation, the split semi-implicit method breaks the full matrix time advance into four sequential time advances, each involving smaller matrices. The use of a high-order triangular element with continuous first derivatives (C1 continuity) allows the Galerkin method to be applied without introduction of new auxiliary variables (such as the vorticity or the current density). These features, along with the manifestly compact nature of the fully node-based C1 finite elements, lead to minimum size matrices for an unconditionally stable method with order of accuracy h4. The resulting matrices are compatible with direct factorization using SuperLU_dist. We demonstrate the accuracy of the method by presenting examples of two-fluid linear wave propagation, two-fluid linear eigenmodes of a tilting cylinder, and of a challenging nonlinear problem in two-fluid magnetic reconnection.

3D MHD VDE and disruptions simulations of tokamaks plasmas including some ITER scenarios

Tokamaks vertical displacement events (VDEs) and disruptions simulations in toroidal geometry by means of a single fluid visco-resistive magneto-hydro-dynamic (MHD) model are presented in this paper. The plasma model is completed with the presence of a 2D wall with finite resistivity which allows the study of the relatively slowly growing magnetic perturbation, the resistive wall mode (RWM), which is, in this paper, the main drive of the disruption evolution. Amplitudes and asymmetries of the halo currents pattern at the wall are also calculated and comparisons with tokamak experimental databases and predictions for ITER are given.

Linear stability and nonlinear dynamics of the fishbone mode in spherical tokamaks

Extensive linear and nonlinear simulations have been carried out to investigate the energetic particle-driven fishbone instability in spherical tokamak plasmas with weakly reversed q profile and the qmin slightly above unity. The global kinetic-MHD hybrid code M3D-K is used. Numerical results show that a fishbone instability is excited by energetic beam ions preferentially at higher qmin values, consistent with the observed appearance of the fishbone before the “long-lived mode” in MAST and NSTX experiments. In contrast, at lower qmin values, the fishbone tends to be stable. In this case, the beam ion effects are strongly stabilizing for the non-resonant kink mode. Nonlinear simulations show that the fishbone saturates with strong downward frequency chirping as well as radial flattening of the beam ion distribution. An (m, n) = (2, 1) magnetic island is found to be driven nonlinearly by the fishbone instability, which could provide a trigger for the (2, 1) neoclassical tearing mode sometimes observed afte...

Nonlinear simulation studies of tokamaks and STs

The multilevel physics, massively parallel plasma simulation code, M3D has been used to study spherical toris (STs) and tokamaks. The magnitude of outboard shift of density profiles relative to electron temperature profiles seen in NSTX under strong toroidal flow is explained. Internal reconnection events in ST discharges can be classified depending on the crash mechanism, just as in tokamak discharges; a sawtooth crash, disruption due to stochasticity, or high-β disruption. Toroidal shear flow can reduce linear growth of internal kink. It has a strong stabilizing effect nonlinearly and causes mode saturation if its profile is maintained, e.g. through a fast momentum source. Normally, however, the flow profile itself flattens during the reconnection process, allowing a complete reconnection to occur. In some cases, the maximum density and pressure spontaneously occur inside the island and cause mode saturation. Gyrokinetic hot particle/MHD hybrid studies of NSTX show the effects of fluid compression on a fast-ion driven n = 1 mode. MHD studies of recent tokamak experiments with a central current hole indicate that the current clamping is due to sawtooth-like crashes, but with n = 0.

Simulation of non-resonant internal kink mode with toroidal rotation in the National Spherical Torus Experiment

Plasmas in spherical and conventional tokamaks, with weakly reversed shear q profile and minimum q above but close to unity, are susceptible to an non-resonant (m,n) = (1,1) internal kink mode. This mode can saturate and persist and can induce a (2,1) seed island for Neoclassical Tearing Mode. [Breslau et al. Nucl. Fusion 51, 063027 (2011)]. The mode can also lead to large energetic particle transport and significant broadening of beam-driven current. Motivated by these important effects, we have carried out extensive nonlinear simulations of the mode with finite toroidal rotation using parameters and profiles of an NTSX plasma with a weakly reversed shear profile. The numerical results show that, at the experimental level, plasma rotation has little effect on either equilibrium or linear stability. However, rotation can significantly influence the nonlinear dynamics of the (1,1) mode and the induced (2,1) magnetic island. The simulation results show that a rotating helical equilibrium is formed and maint...

Three-dimensional modeling of the sawtooth instability in a small tokamak

The sawtooth instability is one of the most fundamental dynamics of an inductive tokamak discharge such as will occur in ITER [R. Aymar et al., Plasma Phys. Controlled Fusion 44, 519 (2002)]. Sawtooth behavior is complex and remains incompletely explained. The Center for Extended MHD Modeling (CEMM) SciDAC project has undertaken an ambitious campaign to model this periodic motion in a small tokamak as accurately as possible using the extended MHD model. Both M3D [W. Park et al., Phys. Plasmas 6, 1796 (1999)] and NIMROD [C. R. Sovinec et al., Phys. Plasmas 10, 1727 (2003)] have been applied to this problem. Preliminary nonlinear MHD results show pronounced stochasticity in the magnetic field following the sawtooth crash but are not yet fully converged. Compared to the MHD model, extended MHD predicts plasma rotation, faster reconnection, and reduced field line stochasticity in the crash aftermath. The multiple time and space scales associated with the reconnection layer and growth time make this an extreme...

Simulation of two fluid and energetic particle effects in stellarators

MHD and resistive MHD are inadequate to understand the stability of stellarators properly. Ideal MHD ballooning mode theory predicts β limits substantially below the values that can be expected in experiments. Resistive MHD is even more pessimistic, predicting that many stellarators are completely unstable. Including two fluid effects, ideally and resistively stable stellarator equilibria can be obtained. It may be possible to completely stabilize ballooning modes. The two fluid computations are done with a realistic value of the Hall parameter, the ratio of the ion skin depth to the major radius. Hybrid gyrokinetic simulations with energetic particles indicate that global shear Alfven TAE modes can be more stable in stellarators than in tokamaks. Computations in a two-period compact stellarator obtained a predominantly n = 1 toroidal mode with the expected TAE frequency. The TAE modes are more stable in the two-period compact stellarator than in a tokamak with the same q and pressure profiles. The cause for the stabilization is believed to be the increased damping rate due to 3D geometry. Simulations were performed with the M3D extended MHD code.

MHD simulations with resistive wall and magnetic separatrix

Abstract A number of problems in resistive MHD magnetic fusion simulations describe plasmas with three regions: the core, the halo region, and the resistive boundary. Treating these problems requires maintenance of an adequate resistivity contrast between the core and halo. This can be helped by the presence of a magnetic separatrix, which in any case is required for reasons of realistic modeling. An appropriate mesh generation capability is also needed to include the halo region when a separatrix is present. Finally a resistive wall boundary condition is required, to allow both two dimensional and three dimensional magnetic perturbations to penetrate the wall. Preliminary work is presented on halo current simulations in ITER. The first step is the study of VDE (vertical displacement event) instabilities. The growth rate is consistent with scaling inversely proportional to the resistive wall penetration time. The simulations have resistivity proportional to the −3/2 power of the temperature. Simulations have been done with resistivity contrast between the plasma core and wall of 1000 times, to model the vacuum region between the core and resistive shell. Some 3D simulations are shown of disruptions competing with VDEs. Toroidal peaking factors are up to about 3.

Simulation studies of the role of reconnection in the “current hole” experiments in the Joint European Torus

Injection of lower hybrid current drive into the current ramp-up phase of Joint European Torus (JET) discharges has been observed to produce an annular current distribution with a core region of essentially zero current density [Hawkes et al., Phys. Rev. Lett. 87, 115001 (2001)]. Similar “current holes” have been observed in Japan Atomic Energy Research Institute Tokamak 60 Upgrade discharges with off-axis current drive supplied by bootstrap current [T. Fujita et al., Phys. Rev. Lett. 87, 245001 (2001)]. In both cases, the central current does not go negative, although current diffusion calculations indicate that there is sufficient noninductive current drive for this to occur. This is explained by the multilevel 3D code (M3D) nonlinear 2D and 3D resistive magnetohydrodynamic (MHD) simulations in toroidal geometry, which predict that these discharges undergo n=0 reconnection events—“axisymmetric sawteeth”—that redistribute the current to hold its core density near zero. Unlike conventional sawteeth, these...

Symmetric solution in M3D

Abstract The coefficient matrices in M3D are reformed here to have symmetric structures. They are further categorized into 3 types: weak diagonally-dominant matrix, moderate diagonally-dominant matrix, and strong diagonally-dominant matrix. The weak diagonally-dominant matrix corresponds to the solution of auxiliary quantity F of the perturbed toroidal flux Ĩ with Neumann boundary conditions. The moderate diagonally-dominant matrix corresponds to the solution of the toroidal current C and the scalar potential Φ with Dirichlet boundary conditions. The strong diagonally-dominant matrix corresponds to the solution of the perturbed toroidal flux Ĩ , the poloidal flux ψ , the pressure p , and the 3 components of velocity: χ , U , and v ϕ with Dirichlet boundary conditions, respectively. We compare LU, GMRES, and ICCG algorithms for linear equation with these 3 types of matrices. It is observed that ICCG greatly accelerates the solution process as compared to GMRES, especially for the weak diagonally-dominant matrix. In this case we achieved 4 to 44 times speedup when the matrix order ranges from ∼10 1 to ∼10 4 . For the moderate diagonally-dominant case, there is a 4 to 24 times speedup. For the strong diagonally-dominant case, an average 2 times speedup is observed. It is also shown that the GMRES algorithm is much slower for the weak diagonally-dominant type than for the other 2 types: 4.4 times slower than the moderate diagonally-dominant case, 333 times slower than the strong diagonally-dominant case. The LU algorithm is faster than GMRES for the weak diagonally-dominant matrix since the matrix is strong ill-conditioned, but GMRES outperforms LU for the strong diagonally-dominant matrix since the matrix is well-conditioned. However, the ICCG algorithm outperforms both ILU and GMRES for all matrix types.