Auxiliary numerical projections of the divergence of flow velocity and vorticity parallel to magnetic field are developed and tested for the purpose of suppressing unphysical interchange instability in magnetohydrodynamic simulations. The numerical instability arises with equal-order C 0 finite-and spectral-element expansions of the flow velocity, magnetic field, and pressure and is sensitive to behavior at the limit of resolution. The auxiliary projections are motivated by physical field-line bending, and coercive responses to the projections are added to the flow-velocity equation. Their incomplete expansions are limited to the highest-order orthogonal polynomial in at least one coordinate of the spectral elements. Cylindrical eigenmode computations show that the projections induce convergence from the stable side with first-order ideal-MHD equations during h-refinement and p-refinement. Hyperbolic and parabolic projections and responses are compared, together with different methods for avoiding magnetic divergence error. The projections are also shown to be effective in linear and nonlinear time-dependent computations with the NIMROD code [C. R. Sovinec, et al., J. Comput. Phys. 195 (2004) 355-386], provided that the projections introduce numerical dissipation.2 filtering of incompressible fluid simulations [7], to a penalty method for MHD eigenvalue computations [4], and to finite/spectral element stabilization techniques [for example, 8-11].Our interest is nonlinear time-dependent simulation of nonideal MHD for magnetic confinement. This model admits viscosity that can stabilize localized interchange and facilitate computations. However, dissipation of dynamics perpendicular to magnetic field is weak in high-temperature plasma, and adding enough conventional viscosity to stabilize numerical interchange can distort dynamics. In principle, it is possible to address numerical interchange with sufficiently fine meshing for physical levels of viscosity. However, resonances that are susceptible to interchange can exist across the entire region of closed magnetic flux, and computing with globally fine resolution is computationally challenging and inefficient. The many resonances that occur in a region of bad magnetic curvature would also foil attempts to stabilize selected modes with spatially localized viscous dissipation.The numerical representation of interchange was first studied in the context of ideal-MHD eigenvalue computation for static equilibria. The fundamental dependent fields are the components of the displacement vector, and perturbations in pressure and magnetic field are eliminated analytically prior to discretization. The resulting force operator is a second-order differential operator in space, and it is self-adjoint [12]. Finite-element methods (FEM) for this problem benefit from distinct expansions for the flux-normal, parallel, and "cross" components of the displacement vector, and integration by parts leaves only first-order derivatives of basis and test functions. To conform with analytical ...