Testing of the new JOREK stellarator-capable model in the tokamak limit

Abstract

In preparation for extending the JOREK nonlinear magnetohydrodynamics (MHD) code to stellarators, a hierarchy of stellarator-capable reduced and full MHD models has been derived and tested. The derivation was presented at the EFTC 2019 conference. Continuing this line of work, we have implemented the reduced MHD model (Nikulsin et al., Phys. Plasmas, vol. 26, 2019, 102109) as well as an alternative model which was newly derived using a different set of projection operators for obtaining the scalar momentum equations from the full MHD vector momentum equation. With the new operators, the reduced model matches the standard JOREK reduced models for tokamaks in the tokamak limit and conserves energy exactly, while momentum conservation is less accurate than in the original model whenever field-aligned flow is present.

1. Introduction

Reduced magnetohydrodynamics (MHD), a system of approximations introduced in its original form in the 1960s (Greene & Johnson 1961), continues to be used by modern MHD codes, such as JOREK (Franck et al. 2015) and M3D-$C^1$ (Breslau, Ferraro & Jardin 2009). While many versions of reduced MHD consisting of different equations and using different methods to derive them have been published, the main idea is the same: the removal of fast magnetosonic waves while retaining relevant physics (Strauss 1997; Jardin et al. 2012). The removal of these waves eliminates the shortest time scale and allows one to use larger time steps due to the Courant condition. Even when implicit time integration methods are used, and the Courant condition is no longer a hard limit, using time steps that are large compared to the shortest time scale can lead to poor accuracy (Kruger, Hegna & Callen 1998; Jardin et al. 2012). An increased time step is especially advantageous in nonlinear simulations, which otherwise tend to be computationally demanding. In particular, stellarator simulations, which are even more computationally expensive than tokamak simulations due to the complex geometry, can benefit from an increased step size.

Originally, reduced MHD was derived via an ordering in a small parameter, often taken to be the inverse aspect ratio. The ordering itself is a system of several approximations and assumptions involving the ordering parameter that allows one to determine the relative order (in terms of the ordering parameter) of any quantity with respect to any other quantity of the same dimension. In this context, terms corresponding to fast magnetosonic waves have a higher order than the terms that one wants to keep, allowing them to be dropped. Different versions of reduced MHD can be derived using different orderings that keep different physical effects (e.g. low-$\beta$ versus high-$\beta$ orderings) (Kruger et al. 1998; Strauss 1976, 1977, 1980, 1997). An alternative ansatz-based approach was introduced by Park et al. (1980), where an ansatz form that eliminates fast magnetosonic waves is used for the velocity and terms of all orders are kept. Later papers also adopt an ansatz for the magnetic field that eliminates field compression (Huysmans & Czarny 2007; Breslau et al. 2009; Franck et al. 2015). The ansatz-based approach makes fewer assumptions and keeps more physical effects, while generally resulting in more complicated equations.

An ordering-based reduced MHD model for stellarators was recently derived in Zocco, Helander & Weitzner (2021) with $\beta \sim \epsilon \ll 1$, where $\epsilon$ is the inverse aspect ratio. Here, as in our previous paper (Nikulsin et al. 2019), we adopt an ansatz-based approach. Previously, we derived a hierarchy of models, which we intended to use for stellarator simulations. We started by introducing ansatzes for the magnetic field and velocity and proving that any arbitrary magnetic field and velocity can be represented by those ansatzes. We also showed that if the background vacuum field is strong enough, then the three terms of velocity ansatz approximately separate the three MHD waves. Dropping the fast magnetosonic wave term from the velocity ansatz and the field compression term from the magnetic field ansatz, we obtain reduced MHD (Nikulsin et al. 2019). However, we later found that the projection operators used had good momentum conservation properties, but not as good energy conservation properties. We thus adopt a different set of projection operators, which makes our new stellarator models a direct generalization of the currently used JOREK tokamak models to three-dimensional geometries; the stellarator models reduce back to the tokamak models in the tokamak limit.

The rest of this paper is structured as follows. In § 2, we briefly review the original set of models, introduce the changes that we made and derive the new set of models. In § 3, we show how the deficiencies of our original model are fixed in the new model. In § 4, we provide some test runs, which show how the deficiencies of the original model manifest themselves in practice by comparing it with the new stellarator model/standard tokamak model (both are the same in this situation) in the case of a tearing mode in a tokamak. In § 5, we consider analytically the local momentum conservation properties of both the original and new models, and then compare the global momentum conservation of the new stellarator model with and without field-aligned flow in the case of a ballooning mode. Finally, in Appendix A, a technicality in the implementation of the original model is discussed.

2. Changes made to the model

The original hierarchy of models, as presented in Nikulsin et al. (2019), was derived from the viscoresistive MHD equations:

(2.1)\begin{equation} \left.\begin{gathered} \frac{\partial \rho}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot} (\rho\boldsymbol{v}) = P, \\ \rho\frac{\partial \boldsymbol{v}}{\partial t} + \rho(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{v} + \boldsymbol{v} P = \boldsymbol{j}\times\boldsymbol{B} - \boldsymbol{\nabla} p + \mu\Delta\boldsymbol{v}, \quad \frac{\partial \boldsymbol{B}}{\partial t} ={-}\boldsymbol{\nabla}\times\boldsymbol{E},\\ \frac{\partial \mathcal{E}}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}\left[\left(\frac{\rho v^2}{2} + \frac{\gamma p}{\gamma-1}\right)\boldsymbol{v} + \frac{p}{\gamma-1}\frac{D_\perp}{\rho}\boldsymbol{\nabla}_\perp\rho + \frac{\boldsymbol{E}\times\boldsymbol{B}}{\mu_0} - \kappa_\perp \boldsymbol{\nabla}_\perp T - \kappa_\parallel \boldsymbol{\nabla}_\parallel T\right] \\ \quad = S_e - \frac{v^2 P}{2}, \\ \mathcal{E} = \frac{\rho v^2}{2} + \frac{p}{\gamma-1} + \frac{B^2}{2\mu_0}, \\ \boldsymbol{\nabla}\times\boldsymbol{B} = \mu_0 \boldsymbol{j}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{B} = 0, \quad \boldsymbol{E} ={-}\boldsymbol{v}\times\boldsymbol{B} + \eta\boldsymbol{j}, \quad P = \boldsymbol{\nabla}\boldsymbol{\cdot} (D_\perp \boldsymbol{\nabla}_{{\perp}} \rho) + S_\rho. \end{gathered}\right\} \end{equation}

The gradient operators parallel and perpendicular to the total magnetic field $\boldsymbol {B}$ are defined as $\boldsymbol {\nabla }_\parallel = ({\boldsymbol {B}}/{B^2})\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }$ and $\boldsymbol {\nabla }_\perp = \boldsymbol {\nabla } - \boldsymbol {\nabla }_\parallel$. The usual MHD notation is followed with $\rho$, $p$, $\boldsymbol {v}$ and $\boldsymbol {B}$ being density, pressure, velocity and magnetic field, respectively. In addition to that, $\eta$ is the resistivity, $\mu$ is the dynamic viscosity, $D_\perp$ is the mass diffusion coefficient perpendicular to field lines, $\kappa _\perp$ and $\kappa _\parallel$ are the thermal conductivity across and along field lines, respectively, and $S_\rho$ and $S_e$ are source terms in the continuity and energy equations, respectively. The ideal gas law $p = \rho RT$ is assumed to hold.

We also introduced the following ansatzes for the magnetic field and velocity:

(2.2)\begin{equation} \left. \begin{gathered} \boldsymbol{B} = \boldsymbol{\nabla}\chi + \boldsymbol{\nabla}\varPsi\times \boldsymbol{\nabla}\chi + \boldsymbol{\nabla}\varOmega\times\boldsymbol{\nabla}\psi_v, \\ \boldsymbol{v} = \frac{\boldsymbol{\nabla}\varPhi\times\boldsymbol{\nabla}\chi}{B_v^2} + v_{{\parallel}}\boldsymbol{B} + \boldsymbol{\nabla}^{{\perp}}\zeta, \end{gathered}\right\} \end{equation}

and proved that any magnetic field and any velocity can be represented in such a way. Here, $\boldsymbol {\nabla }\chi$ is the curl-free vacuum component of the magnetic field, which is generated by the coils, $B_v = |\boldsymbol {\nabla }\chi |$ and $\boldsymbol {\nabla }^{\perp } = \boldsymbol {\nabla } - B_v^{-2}\boldsymbol {\nabla }\chi (\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla })$. Note that, since the ansatzes do not restrict the magnetic and velocity fields, they can still be used within the context of full MHD. The transition to reduced MHD is done by setting $\varOmega = \zeta = 0$ and dropping the evolution equations for $\varOmega$ and $\zeta$ (derived below). Further, in the tokamak limit, i.e. when $\chi = F_0\phi$, where $\phi$ is the toroidal angle,¹ the ansatzes reduce to those of the standard JOREK reduced MHD tokamak models, $\boldsymbol {B} = F_0\boldsymbol {\nabla }\phi + \boldsymbol {\nabla }\psi \times \boldsymbol {\nabla }\phi$ and $\boldsymbol {v} = R^2\boldsymbol {\nabla } u\times \boldsymbol {\nabla }\phi + v_{\parallel }\boldsymbol {B}$, where $\psi = F_0\varPsi$ and $u = \varPhi /F_0$ (Hoelzl et al. 2021).

The scalar $\psi _v$ is one of the Clebsch potentials for the vacuum field, which can locally be written as

(2.3)\begin{equation} \boldsymbol{\nabla}\chi = \boldsymbol{\nabla}\psi_v\times\boldsymbol{\nabla}\beta_v. \end{equation}

As we discussed in Nikulsin et al. (2019), one can construct a Clebsch-type coordinate system with the coordinates $(\psi _v,\beta _v,\chi )$. Note that the variables $\psi _v$ and $\beta _v$ are not unique, as can be seen by replacing $\psi _v$ with $\psi _v' = \psi _v + f(\beta _v)$, where $f$ is an arbitrary function, which leaves (2.3) unchanged. In general, given some vacuum magnetic field $\boldsymbol {\nabla }\chi$, the scalars $\psi _v$ and $\beta _v$ are completely determined by fixing some parameterization $(\psi _v,\tilde {\beta }_v)$ of a surface intersected by the vacuum field lines, with the values of $\psi _v$ and $\tilde {\beta }_v$ in the rest of space being determined by following the field lines off of the surface while keeping $\psi _v$ and $\tilde {\beta }_v$ constant. In a toroidal system, the surface is usually chosen to be a poloidal plane, and, in many cases, a cut needs to be introduced at that poloidal plane to prevent multivaluedness of $\psi _v$ and $\tilde {\beta }_v$ when the field lines travel once around the torus and encounter previously labelled points. Finally, $\beta _v$ is determined by integrating $\partial \beta _v/\partial \tilde {\beta }_v = B_v/|\boldsymbol {\nabla }\psi _v\times \boldsymbol {\nabla }\tilde {\beta }_v|$ (Stern 1970; D'haeseleer et al. 2012). Usually, it is advantageous to find a $\psi _v$ such that the component of $\boldsymbol {\nabla }\varOmega \times \boldsymbol {\nabla }\psi _v$ perpendicular to $\boldsymbol {\nabla }\chi$, $F_v^2|\partial ^{\parallel }\varOmega |^2$, is minimized at $t = 0$, so that most of field line bending is captured by the $\boldsymbol {\nabla }\varPsi \times \boldsymbol {\nabla }\chi$ term.

The variables $\varPsi$, $\varOmega$, $\varPhi$, $v_{\parallel }$ and $\zeta$ are the unknown magnetic field and velocity variables, which are solved for. In order to obtain evolution equations for these variables, we projected Faraday's law onto $\boldsymbol {\nabla }\psi _v$ and $\boldsymbol {\nabla }\chi$, and applied the following projection operators to the vector momentum equation, which was first divided by $\rho$:

(2.4)\begin{equation} \left. \begin{gathered} \boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{\nabla}\times[\boldsymbol{\nabla} \chi\times(\boldsymbol{e}_\chi\times \\ \boldsymbol{\nabla}\chi\boldsymbol{\cdot} \\ \boldsymbol{\nabla}\boldsymbol{\cdot} [B_v^2 \boldsymbol{\nabla}\chi\times(\boldsymbol{e}_\chi\times \end{gathered}\right\} \end{equation}

Here, $\boldsymbol {e}_\chi = \boldsymbol {B}/B^\chi$ is the covariant basis vector in the Clebsch-type coordinate system $(\alpha ,\beta ,\chi )$ associated with the full magnetic field, $\boldsymbol {B} = \boldsymbol {\nabla }\alpha \times \boldsymbol {\nabla }\beta$. These projection operators have the property that each one of them is orthogonal to all but one term in the velocity ansatz (2.2), i.e.

(2.5)\begin{equation} \left. \begin{gathered} \boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{\nabla}\times[ \boldsymbol{\nabla}\chi\times(\boldsymbol{e}_\chi\times\boldsymbol{v})] = \varDelta^{{\perp}}\varPhi, \\ \boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{v} = B^\chi v_{{\parallel}}, \\ \boldsymbol{\nabla}\boldsymbol{\cdot} [B_v^2\boldsymbol{\nabla}\chi\times(\boldsymbol{e}_\chi \times\boldsymbol{v})] ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}(B_v^2 \boldsymbol{\nabla}^{{\perp}} \zeta), \end{gathered}\right\} \end{equation}

where $\varDelta ^{\perp } = \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\nabla }^{\perp }$. However, for the sake of energy conservation, a better set of projection operators would be

(2.6)\begin{equation} \left. \begin{gathered} \boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{\nabla}\times(B_v^{{-}2} \\ \boldsymbol{B}\,\boldsymbol{\cdot} \\ \boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{\nabla}\times(B_v^{{-}2}\boldsymbol{\nabla}\chi\times \end{gathered}\right\} \end{equation}

as we show in the next section.

To obtain the new evolution equations for the velocity variables, we insert the ansatzes (2.2) into the momentum equation and apply the projection operators (2.6), this time without dividing first by $\rho$. The momentum equation can be written as

(2.7)\begin{equation} \rho\boldsymbol{\nabla}\frac{\partial \varPhi}{\partial t}\times\frac{\boldsymbol{\nabla}\chi}{B_v^2} + \rho\frac{\partial }{\partial t}(v_{{\parallel}}\boldsymbol{B}) + \rho\boldsymbol{\nabla}^{{\perp}}\frac{\partial \zeta}{\partial t} + \frac{\rho}{2}\boldsymbol{\nabla} v^2 + \rho\boldsymbol{\omega}\times\boldsymbol{v} + \boldsymbol{v} P = \boldsymbol{j}\times\boldsymbol{B} - \boldsymbol{\nabla} p, \end{equation}

where we used the identity $(\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v} = \frac {1}{2}\boldsymbol {\nabla } v^2 + \boldsymbol {\omega }\times \boldsymbol {v}$ and $\boldsymbol {\omega }$ is the vorticity:

(2.8)\begin{align} \boldsymbol{\omega} &= \boldsymbol{\nabla}\times\boldsymbol{v} ={-}\boldsymbol{\nabla}\chi\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\boldsymbol{\nabla}\varPhi}{B_v^2}\right) + B_v\partial^{{\parallel}}\frac{\boldsymbol{\nabla}\varPhi}{B_v^2} - \frac{\boldsymbol{\nabla}\varPhi\boldsymbol{\cdot}\boldsymbol{\nabla}}{B_v^2}\boldsymbol{\nabla}\chi + \boldsymbol{\nabla} v_{{\parallel}}\times\boldsymbol{B} + \mu_0 v_{{\parallel}}\,\boldsymbol{j}\nonumber\\ &\quad + \boldsymbol{\nabla}\chi\times\boldsymbol{\nabla}\frac{\partial^{{\parallel}}\zeta}{B_v}. \end{align}

Just as in Nikulsin et al. (2019), we do not carry the viscosity term through this derivation, instead adding a generic viscosity term afterwards. To obtain the evolution equation for $\varPhi$, we apply to (2.7) the first projection operator in (2.6), or its equivalent, $-\boldsymbol {\nabla }\boldsymbol {\cdot }(B_v^{-2}\boldsymbol {\nabla }\chi \times$, when appropriate:

(2.9)\begin{align} &\boldsymbol{\nabla}\boldsymbol{\cdot}\left[\frac{\rho}{B_v^2}\boldsymbol{\nabla}^{{\perp}}\frac{\partial \varPhi}{\partial t} + \rho\frac{\partial }{\partial t}(v_{{\parallel}}\boldsymbol{\nabla}^{{\perp}}\varPsi) - \frac{\rho}{B_v}\frac{\partial }{\partial t}(v_{{\parallel}}\partial^{{\parallel}}\varOmega)\boldsymbol{\nabla}\psi_v\right] - B_v\left[\frac{\rho}{B_v^2},\frac{\partial \zeta}{\partial t}\right] = \frac{B_v}{2}\left[\frac{\rho}{B_v^2},v^2\right] \nonumber\\ &\quad + B_v\left[\frac{\rho\omega^\chi}{B_v^4},\varPhi\right] - \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho v_{{\parallel}}\boldsymbol{\omega}^\perp) + B_v\left[\frac{\rho v_{{\parallel}}\omega^\chi}{B_v^2},\varPsi\right] + B_v\left[\frac{\rho v_{{\parallel}}\omega^{\psi_v}}{B_v^2},\varOmega\right] \nonumber\\ &\quad - B_v\left[\frac{\rho v_{{\parallel}}\boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}\varOmega}{B_v^2},\psi_v\right] + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\rho\omega^\chi}{B_v^2}\boldsymbol{\nabla}^{{\perp}}\zeta - \frac{P}{B_v^2}\boldsymbol{\nabla}^{{\perp}}\varPhi - Pv_{{\parallel}}\boldsymbol{\nabla}^{{\perp}}\varPsi + \frac{Pv_{{\parallel}}\partial^{{\parallel}}\varOmega}{B_v}\boldsymbol{\nabla}\psi_v\right) \nonumber\\ &\quad + B_v\left[\frac{P}{B_v^2},\zeta\right] + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{B^\chi}{B_v^2}\boldsymbol{j} - \frac{j^\chi}{B_v^2}\boldsymbol{B}\right) + B_v\left[\frac{1}{B_v^2},p\right] + \boldsymbol{\nabla}\boldsymbol{\cdot}(\mu_\perp\boldsymbol{\nabla}^{{\perp}}\varDelta^{{\perp}}\varPhi), \end{align}

where $[f,g] = B_v^{-1}\boldsymbol {\nabla }\chi \boldsymbol {\cdot }(\boldsymbol {\nabla } f\times \boldsymbol {\nabla } g)$ is a Poisson bracket, $\partial ^{\parallel } = B_v^{-1}\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla }$, superscripts indicate contravariant vector components and, for any vector $\boldsymbol {U}$, $\boldsymbol {U}^\perp = \boldsymbol {U} - U^\chi \boldsymbol {\nabla }\chi /B_v^2$. We have added the generic viscosity term ${\boldsymbol {\nabla }\boldsymbol {\cdot }(\mu _\perp \boldsymbol {\nabla }^{\perp }\varDelta ^{\perp }\varPhi )}$, which we have chosen to match the viscosity term in Franck et al. (2015) in the tokamak limit. As discussed in previous papers (Franck et al. 2015; Nikulsin et al. 2019), the viscosity term $\mu \Delta \boldsymbol {v}$ in the momentum equation (2.1) is only a rough approximation for the divergence of the viscous stress tensor in a plasma, and one does not lose much by using a generic viscosity term in the final equations. Note, however, that the $\mu _\perp$ in the above equation is not the physical viscosity $\mu$ in (2.1). Indeed, from dimensional analysis it follows that $\mu _\perp$ has units of T$^2$ Pa s. Applying scaling analysis to the term $\mu \Delta \boldsymbol {v}$ after inserting just the first term of the ansatz (2.2) for $\boldsymbol {v}$ and applying the first projection operator in (2.6), we see that $\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla }\times [B_v^{-2}\mu \varDelta (\boldsymbol {\nabla }\varPhi \times \boldsymbol {\nabla }\chi /B_v^2)] \sim \mu B_v^{-2}\varPhi /L_\perp ^4$, where $L_\perp$ is the length scale perpendicular to the magnetic field, and $L_\perp \ll L_\parallel$ for most magnetic fusion configurations. Applying the same analysis to the generic viscosity term, we get $\boldsymbol {\nabla }\boldsymbol {\cdot }(\mu _\perp \boldsymbol {\nabla }^{\perp }\varDelta ^{\perp }\varPhi ) \sim \mu _\perp \varPhi /L_\perp ^4$. Comparing the two scalings, we see that $\mu _\perp$ scales as $\mu B_v^{-2}$; for the purposes of the scaling, one can take the value of $B_v$ at the axis as a typical value and write $\mu _\perp \sim \mu B_{v,\textrm {axis}}^{-2}$.

To get the evolution equation for $v_{\parallel }$, we apply the second operator in (2.6), projecting equation (2.7) on $\boldsymbol {B}$:

(2.10)\begin{align} &\rho\left(\frac{\partial \varPhi}{\partial t},\varPsi\right) \!-\! \rho\frac{F_v}{B_v}\partial^{{\parallel}}\varOmega\partial_{\psi_v}\frac{\partial \varPhi}{\partial t} \!+\! \rho B^2\frac{\partial v_{{\parallel}}}{\partial t} + \frac{\rho v_{{\parallel}}}{2}\frac{\partial B^2}{\partial t} + \rho B_v\left[\frac{\partial \zeta}{\partial t},\varPsi\right] + \rho F_v\left[\frac{\partial \zeta}{\partial t},\varOmega\right]_{\psi_v} \nonumber\\ &\quad ={-}\frac{\rho B_v}{2}\partial^{{\parallel}} v^2 - \frac{\rho B_v}{2}\left[v^2,\varPsi\right] - \frac{\rho F_v}{2}\left[v^2,\varOmega\right]_{\psi_v} - \frac{\rho\omega^\chi}{B_v}\partial^{{\parallel}}\varPhi - \frac{\rho\omega^\chi}{B_v}\left[\varPhi,\varPsi\right] \nonumber \\ &\qquad - \frac{\rho\omega^\chi F_v}{B_v^2}\left[\varPhi,\varOmega\right]_{\psi_v} + \frac{\rho B^\chi}{B_v^2}\boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}\varPhi - \rho\boldsymbol{\omega}\boldsymbol{\cdot}(\boldsymbol{\nabla}\zeta\times\boldsymbol{\nabla}\chi) + \rho\omega^\chi(\varPsi,\zeta) \nonumber\\ &\qquad - \rho\boldsymbol{\omega}\boldsymbol{\cdot}\boldsymbol{\nabla}\varOmega F_v\partial_{\psi_v}\zeta + \rho\omega^{\psi_v}(\varOmega,\zeta) - \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{B}\,P \nonumber\\ &\qquad - B_v\partial^{{\parallel}} p - B_v\left[p,\varPsi\right] - F_v\left[p,\varOmega\right]_{\psi_v} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\mu_\parallel\boldsymbol{\nabla}^{{\perp}} v_{{\parallel}}), \end{align}

where we have again added a generic viscosity term with a form analogous to the viscosity term in the previous equation. Here, $(f,g) = \boldsymbol {\nabla }^{\perp } f\boldsymbol {\cdot }\boldsymbol {\nabla }^{\perp } g$ is an inner product, $F_v = |\boldsymbol {\nabla }\psi _v|$, $\partial _{\psi _v} = F_v^{-1}\boldsymbol {\nabla }\psi _v\boldsymbol {\cdot }\boldsymbol {\nabla }$ and $[f,g]_{\psi _v} = F_v^{-1}\boldsymbol {\nabla }\psi_v \boldsymbol {\cdot }(\boldsymbol {\nabla } f\times \boldsymbol {\nabla } g)$ is a Poisson bracket with respect to $\boldsymbol {\nabla }\psi _v$. Finally, to obtain the evolution equation for $\zeta$, we apply the last projection operator in (2.6) (or its equivalent, ${-\boldsymbol {\nabla }\boldsymbol {\cdot }[B_v^{-2}\boldsymbol {\nabla }\chi \times (\boldsymbol {\nabla }\chi \times }$, to some terms) to equation (2.7):

(2.11)\begin{align} &B_v\left[\frac{\rho}{B_v^2},\frac{\partial \varPhi}{\partial t}\right] + B_v\left[\rho\frac{\partial v_{{\parallel}}}{\partial t},\varPsi\right] + B_v\left[\rho v_{{\parallel}},\frac{\partial \varPsi}{\partial t}\right] - B_v\left[\frac{\rho}{B_v}\frac{\partial }{\partial t}(v_{{\parallel}}\partial^{{\parallel}}\varOmega),\psi_v\right]\nonumber\\ &\quad + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\rho\boldsymbol{\nabla}^{{\perp}}\frac{\partial \zeta}{\partial t}\right) ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\rho}{2}\boldsymbol{\nabla}^{{\perp}} v^2 + \frac{\rho\omega^\chi}{B_v^2}\boldsymbol{\nabla}^{{\perp}}\varPhi - \frac{\rho v_{{\parallel}} B^\chi}{B_v^2}\boldsymbol{\nabla}\chi\times\boldsymbol{\omega}\right.\nonumber\\ &\quad\left. + \frac{\rho v_{{\parallel}}\omega^\chi}{B_v^2}\boldsymbol{\nabla}\chi\times\boldsymbol{B}\right) + B_v\left[\frac{\rho\omega^\chi}{B_v^2},\zeta\right] - B_v\left[\frac{P}{B_v^2},\varPhi\right] - B_v\left[Pv_{{\parallel}},\varPsi\right]\nonumber\\ &\quad + B_v\left[\frac{Pv_{{\parallel}}\partial^{{\parallel}}\varOmega}{B_v},\psi_v\right] - \boldsymbol{\nabla}\boldsymbol{\cdot}\left(P\boldsymbol{\nabla}^{{\perp}}\zeta + \frac{B^\chi}{B_v^2}\boldsymbol{\nabla}\chi\times\boldsymbol{j} - \frac{j^\chi}{B_v^2}\boldsymbol{\nabla}\chi\times\boldsymbol{B}\right)\nonumber\\ &\quad - \varDelta^{{\perp}} p + \boldsymbol{\nabla}\boldsymbol{\cdot}(\mu_\zeta\boldsymbol{\nabla}^{{\perp}}\varDelta^{{\perp}}\zeta). \end{align}

This concludes the derivation of the scalar momentum equations.

The evolution equations for the magnetic field variables of Nikulsin et al. (2019) are valid; therefore we keep them as they are:

(2.12)\begin{equation} \left. \begin{gathered} \left[\psi_v,\frac{\partial \varPsi}{\partial t}\right] = \left[\frac{[\varPsi,\varPhi]\!-\!\partial^{{\parallel}}\varPhi}{B_v},\psi_v\right] \!-\! \frac{F_v}{B_v}\left[\varOmega,\frac{[\psi_v,\varPhi]}{B_v}\right]_{\psi_v} + \partial^{{\parallel}}(F_v\partial_{\psi_v}\zeta) + \left[(\zeta,\varPsi),\psi_v\right] \\ - \frac{F_v}{B_v}\left[\varOmega,(\zeta,\psi_v)\right]_{\psi_v} + \frac{1}{B_v}\boldsymbol{\nabla}\boldsymbol{\cdot}(\eta\boldsymbol{\nabla}\psi_v\times\boldsymbol{j}), \\ \left[\frac{\partial \varOmega}{\partial t},\psi_v\right] ={-}\left[\varOmega,\frac{[\psi_v,\varPhi]}{B_v}\right] + \left[\psi_v,\frac{[\varOmega,\varPhi]}{B_v}\right] - 2(B_v,\zeta) - B_v\varDelta^{{\perp}}\zeta - \left[\varOmega,(\zeta,\psi_v)\right] \\ + \left[\psi_v,(\zeta,\varOmega)\right] + \frac{1}{B_v}\boldsymbol{\nabla} \boldsymbol{\cdot}(\eta\boldsymbol{\nabla}\chi\times\boldsymbol{j}). \end{gathered}\right\} \end{equation}

The continuity equation remains unchanged (Nikulsin et al. 2019):

(2.13)\begin{equation} \frac{\partial \rho}{\partial t} ={-} B_v\left[\frac{\rho}{B_v^2},\varPhi\right] - B_v\partial^{{\parallel}}(\rho v_{{\parallel}}) - B_v[\rho v_{{\parallel}},\varPsi] - F_v[\rho v_{{\parallel}},\varOmega]_{\psi_v} - \boldsymbol{\nabla}\boldsymbol{\cdot}(\rho\boldsymbol{\nabla}^{{\perp}}\zeta) + P, \end{equation}

while for the pressure, we use the standard internal energy evolution equation, instead of the full energy conservation equation from the system (2.1):

(2.14)\begin{align} &\frac{\partial p}{\partial t} + \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} p + \gamma p\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{v} \nonumber\\ &\quad = (\gamma-1)\left[\boldsymbol{\nabla}\boldsymbol{\cdot}\left(\kappa_\perp\boldsymbol{\nabla}_\perp T + \kappa_\parallel\boldsymbol{\nabla}_\parallel T + \frac{p}{\gamma - 1}\frac{D_\perp}{\rho}\boldsymbol{\nabla}_\perp\rho\right) + S_e + \eta j^2\right], \end{align}

which now, after using (2.2), reads

(2.15)\begin{align} \frac{\partial p}{\partial t} &={-} \frac{1}{B_v}\left[p,\varPhi\right] - v_{{\parallel}} B_v\partial^{{\parallel}} p - v_{{\parallel}} B_v\left[p,\varPsi\right] - v_{{\parallel}} F_v \left[p,\varOmega\right]_{\psi_v} - (\zeta,p) - \gamma p B_v\left[\frac{1}{B_v^2},\varPhi\right] \nonumber\\ &\quad- \gamma p B_v\partial^{{\parallel}} v_{{\parallel}} - \gamma p B_v\left[v_{{\parallel}},\varPsi\right] - \gamma p F_v\left[v_{{\parallel}},\varOmega\right]_{\psi_v} - \gamma p\varDelta^{{\perp}}\zeta \nonumber\\ &\quad + \boldsymbol{\nabla}\boldsymbol{\cdot}\left[(\gamma-1)\frac{\kappa_\perp}{R}\boldsymbol{\nabla} \left(\frac{p}{\rho}\right)\right.\nonumber\\ &\quad + (\gamma-1)\frac{\kappa_\parallel{-}\kappa_\perp}{RB^2} \boldsymbol{B}\left(B_v\partial^{{\parallel}}\left(\frac{p}{\rho}\right) + B_v\left[\frac{p}{\rho},\varPsi\right] + F_v\left[\frac{p}{\rho},\varOmega\right]_{\psi_v}\right)\nonumber\\ &\quad \left. + \frac{p D_\perp}{\rho}\left(\boldsymbol{\nabla}\rho-\frac{\boldsymbol{B}}{B^2}(B_v\partial^{{\parallel}}\rho + B_v[\rho,\varPsi] + F_v[\rho,\varOmega]_{\psi_v})\right)\right] + (\gamma-1)(S_e + \eta j^2). \end{align}

We note that since any arbitrary magnetic field and velocity can be represented in the ansatz form (2.2), the equations derived above still correspond to full MHD, albeit in a unconventional form. We also showed that if the background vacuum field is the dominant part of the magnetic field, then the MHD waves are approximately separated by the terms of the velocity ansatz, with the first term containing Alfvén waves, the second term containing slow magnetosonic waves and the last term containing fast magnetosonic waves (Nikulsin et al. 2019). Thus, to get the reduced MHD system, we set $\zeta = \varOmega = 0$, eliminating fast magnetosonic waves and field compression, and drop the corresponding evolution equations, namely (2.11) and the second equation in (2.12), to avoid having an overconstrained system. We also approximate the electric field in Ohm's law as $\boldsymbol {E} = -\boldsymbol {v}\times \boldsymbol {B} + \eta \,\boldsymbol {j}^\parallel$, which corresponds to neglecting the components of $\boldsymbol {j}$ perpendicular to $\boldsymbol {\nabla }\chi$ in the last term of the first equation in (2.12). Alternatively, if we introduce a tensor resistivity, $\boldsymbol {\eta }$, this approximation corresponds to neglecting the perpendicular resistivity, $\eta ^\perp \approx 0$, although parallel and perpendicular resistivities are usually of the same order, so it makes more sense to speak about neglecting the perpendicular current (only in the resistive term). This seems to be a fairly good approximation, at least for tokamaks (Pfirsch 1978), and likely also for any configuration with significant bootstrap current. Indeed, we can write $\boldsymbol {\eta }\boldsymbol {\cdot }\boldsymbol {j} = \eta _\perp \boldsymbol {j}_\perp + \eta _\parallel \,\boldsymbol {j}_\parallel$, where subscripts refer to parallel and perpendicular components relative to the total field $\boldsymbol {B}$, not just $\boldsymbol {\nabla }\chi$. When the bootstrap current is high enough, the $\eta _\parallel \boldsymbol {j}_\parallel$ term is dominant, and, since $\boldsymbol {\nabla }\chi$ is the dominant component of $\boldsymbol {B}$, we have $\boldsymbol {j}_\parallel \approx \boldsymbol {j}^\parallel$ in the lowest order. As we show in the next section, the approximation is necessary in order to have energy conservation in reduced MHD. Accordingly, we also have to replace the last term of equation (2.15) with $\eta (\kern1.5pt j^\parallel )^2$ to avoid creating extra internal energy, as now only the parallel component of current contributes to the dissipation of magnetic energy. A further reduction would be to also eliminate slow magnetosonic waves by setting $v_{\parallel } = 0$ and dropping equation (2.10).

Finally, we point out that in the tokamak limit, when $\chi = F_0\phi$, the new reduced models derived in this section match the currently implemented JOREK tokamak models (Hoelzl et al. 2021). The $\varPsi$ evolution equation, once we set $\zeta = \varOmega = 0$ and drop the perpendicular components of the current, becomes

(2.16)\begin{equation} \left[\psi_v,\frac{\partial \varPsi}{\partial t}\right] = \left[\frac{[\varPsi,\varPhi]-\partial^{{\parallel}}\varPhi}{B_v},\psi_v\right] + \frac{1}{B_v}\boldsymbol{\nabla}\boldsymbol{\cdot}(\eta\boldsymbol{\nabla}\psi_v \times\boldsymbol{j}^\parallel), \end{equation}

which can be rewritten as

(2.17)\begin{equation} \boldsymbol{\nabla}\psi_v\boldsymbol{\cdot}\boldsymbol{\nabla} \times\left(\frac{\partial \varPsi}{\partial t}\boldsymbol{\nabla}\chi\right) = \boldsymbol{\nabla}\psi_v\boldsymbol{\cdot}\boldsymbol{\nabla} \times\left(\frac{\partial^{{\parallel}}\varPhi - \left[\varPsi,\varPhi\right]}{B_v} \boldsymbol{\nabla}\chi - \eta\frac{j^\chi}{B_v^2}\boldsymbol{\nabla}\chi\right). \end{equation}

The above equation will be satisfied as long as

(2.18)\begin{equation} B_v\frac{\partial \varPsi}{\partial t} = \partial^{{\parallel}}\varPhi - \left[\varPsi,\varPhi\right] - \eta\frac{j^\chi}{B_v} \end{equation}

is satisfied, which exactly matches the currently implemented form of the magnetic field equation in JOREK when $\chi = F_0\phi$ (Hoelzl et al. 2021). For all other equations, the fact that they match the corresponding JOREK reduced MHD equations in the tokamak limit is obvious and can be verified by simple substitution of $\chi = F_0\phi$.

3. Energy conservation

We begin by showing that the new models conserve energy. If we apply the first projection operator in (2.6) to some vector field $\boldsymbol {Q}$, multiply the result by a test function $\varPhi ^*$ and integrate over the plasma volume, we can get the following expression using the identity $\boldsymbol {\nabla } f\boldsymbol {\cdot }\boldsymbol {\nabla }\times \boldsymbol {U} = -\boldsymbol {\nabla }\boldsymbol {\cdot }(\boldsymbol {\nabla } f\times \boldsymbol {U})$ and integration by parts:

(3.1)\begin{equation} \int_V \varPhi^*\boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{\nabla}\times\left(\frac{ \boldsymbol{Q}}{B_v^2}\right)\textrm{d}V = \int_V \frac{\boldsymbol{\nabla}\varPhi^*\times\boldsymbol{\nabla}\chi}{B_v^2}\boldsymbol{\cdot}\boldsymbol{Q}\,\textrm{d}V, \end{equation}

where we let $\varPhi ^* = 0$ on $\partial V$, so the surface integral term is zero. Doing the same with the third projection operator and test function $\zeta ^*$, we have

(3.2)\begin{align} \int_V \zeta^*\boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{\nabla}\times\left(\frac{\boldsymbol{\nabla}\chi\times \boldsymbol{Q}}{B_v^2}\right)\textrm{d}V &={-}\int_V \frac{\boldsymbol{\nabla}\zeta^*}{B_v^2}\boldsymbol{\cdot}[\boldsymbol{\nabla}\boldsymbol{\chi}\times( \boldsymbol{Q}\times\boldsymbol{\nabla}\chi)]\,\textrm{d}V \nonumber\\ &={-}\int_V \boldsymbol{\nabla}^{{\perp}}\zeta^*\boldsymbol{\cdot}\boldsymbol{Q}\,\textrm{d}V. \end{align}

We can also apply the second projection operator (2.6) to $\boldsymbol {Q}$, multiply by a test function $v_{\parallel }^*$ and integrate over the plasma volume:

(3.3)\begin{equation} \int_V v_{{\parallel}}^*\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{Q} \,\textrm{d}V = \int_V v_{{\parallel}}^*\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{Q} \,\textrm{d}V. \end{equation}

Now, if we let $\boldsymbol {Q}$ be the vector momentum equation, $\boldsymbol {Q} \equiv \rho \partial \boldsymbol {v}/\partial t + \rho (\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v} + \boldsymbol {v} P - \boldsymbol {j}\times \boldsymbol {B} + \boldsymbol {\nabla } p = 0$, as well as $\varPhi ^* = \varPhi$, $v_{\parallel }^* = v_{\parallel }$ and $\zeta ^* = -\zeta$, and sum up the above three equations, we get

(3.4)\begin{equation} 0 = \int_V \boldsymbol{v}\boldsymbol{\cdot}\left[\rho\frac{\partial \boldsymbol{v}}{\partial t} + \rho(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{v} + \boldsymbol{v} P - \boldsymbol{j}\times\boldsymbol{B} + \boldsymbol{\nabla} p\right]\textrm{d}V, \end{equation}

where the left-hand side is zero due to equations (2.9), (2.10) and (2.11). More importantly, if we set $\zeta = 0$ and drop equation (2.11), (3.4) still holds, since with $\zeta ^* = -\zeta = 0$ equation (3.2) becomes $0 = 0$, and (2.11) no longer has to be satisfied in order for the left-hand side of (3.4) to remain zero. Similarly, if we set $v_{\parallel } = 0$ and drop equation (2.10), (3.4) continues to hold due to equation (3.3) becoming $0 = 0$ and no longer requiring equation (2.10) to be satisfied. This property of equation (3.4) being satisfied even when some terms in the velocity ansatz are missing was not present in the original models, leading to problems which we describe in detail further in this section.

We now proceed to the magnetic field evolution equations. In the full MHD case, both of the equations (2.12) are satisfied. These two equations are simply the vector components of Faraday's law in the $\boldsymbol {\nabla }\psi _v$ and $\boldsymbol {\nabla }\chi$ directions, and as such are equivalent to the original Faraday's law in vector form. Only two equations are needed to evolve the full magnetic field, as one degree of freedom is eliminated by the $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B} = 0$ constraint. In the reduced MHD case, Faraday's law is also satisfied, as one can easily see from multiplying equation (2.18) by $\boldsymbol {\nabla }\chi /B_v$ and taking the curl. Thus, we can take the dot product of Faraday's law with $\boldsymbol {B}$ and follow the usual procedure to obtain the magnetic energy equation:

(3.5)\begin{equation} \frac{1}{2\mu_0}\frac{\partial B^2}{\partial t} + \frac{1}{\mu_0}\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{E}\times\boldsymbol{B}) ={-}\boldsymbol{v}\boldsymbol{\cdot}(\kern1.5pt \boldsymbol{j}\times\boldsymbol{B}) - \eta(\kern1.5pt j^*)^2, \end{equation}

where we have $j^* = j$ in the full MHD case and $j^* = j^\parallel$ in the reduced MHD case. To complete the derivation of the energy conservation law, we rewrite (3.4) as

(3.6)\begin{equation} \int_V \left[\frac{\partial }{\partial t}\left(\frac{\rho v^2}{2}\right) + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\rho v^2}{2}\boldsymbol{v}\right)\right]\textrm{d}V = \int_V \left[\boldsymbol{v}\boldsymbol{\cdot}(\kern1.5pt \boldsymbol{j}\times\boldsymbol{B}) - \boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} p - \frac{v^2}{2}P\right]\textrm{d}V, \end{equation}

where we used the continuity equation. We then divide (2.14) by $\gamma -1$ and integrate it, together with (3.5) over the plasma volume. Adding both of the resulting equations to equation (3.6), we obtain

(3.7)\begin{align} &\int_V \left[\frac{\partial }{\partial t}\left(\frac{\rho v^2}{2} + \frac{p}{\gamma-1} + \frac{B^2}{2\mu_0}\right) + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\rho v^2}{2}\boldsymbol{v} + \frac{\gamma p}{\gamma-1}\boldsymbol{v} + \frac{\boldsymbol{E}\times\boldsymbol{B}}{\mu_0} \right.\right. \nonumber\\ &\left.\left.\quad - \kappa_\perp\boldsymbol{\nabla}_\perp T - \kappa_\parallel\boldsymbol{\nabla}_\parallel T - \frac{p}{\gamma - 1}\frac{D_\perp}{\rho}\boldsymbol{\nabla}_\perp\rho\right)\right]\textrm{d}V = \int_V \left(S_e - \frac{v^2}{2}P\right)\textrm{d}V, \end{align}

which can be rewritten as

(3.8)\begin{equation} \frac{\textrm{d}E}{\textrm{d}t} = \oint_{\partial V} \left(\kappa_\perp\boldsymbol{\nabla}_\perp T + \frac{RT}{\gamma - 1}D_\perp\boldsymbol{\nabla}_\perp\rho\right)\boldsymbol{\cdot} \textrm{d}\boldsymbol{S} + \int_V \left(S_e - \frac{v^2}{2}P\right)\textrm{d}V. \end{equation}

Here, we have assumed that the plasma is surrounded by a perfectly conducting wall, which implies the boundary conditions $\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {n} = 0$ and $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {n} = 0$, where $\boldsymbol {n}$ is a unit normal vector to the boundary. These conditions, along with $\varPhi = \zeta = 0$ on $\partial V$, which we had assumed before to make the boundary integral terms in (3.1) and (3.2) zero, are consistent with the boundary conditions that were used in the simulations presented in the following sections.

It is important to point out that, if we replace $P$ by $\partial \rho /\partial t + \boldsymbol {\nabla }\boldsymbol {\cdot }(\rho \boldsymbol {v})$, exact energy conservation continues to hold even when the exact solutions are replaced by Bubnov–Galerkin finite-element approximations. Indeed, a Bubnov–Galerkin solution will ensure that the expressions (3.1), (3.2) and (3.3) equal zero whenever the test functions $\varPhi ^*$, $\zeta ^*$ and $v_{\parallel }^*$ are finite-element basis functions. Since the finite-element solutions $\varPhi$, $\zeta$ and $v_{\parallel }$ are linear combinations of the basis functions, (3.6) will continue to hold. Similarly, since the $\varPsi$ and $\varOmega$ solutions will be linear combinations of the basis functions and the Bubnov–Galerkin method ensures that the weak form of (2.12) is satisfied whenever the test function is a basis function, the integral of (3.5) over volume will continue to hold. Finally, since the constant $1$ can be represented as a linear combination of basis functions, the integral of (2.14) over volume will continue to hold, and so energy will be exactly conserved. Note, however, that this argument does not take into account temporal discretization and Fourier expansion in the toroidal direction, both of which contribute to energy conservation error in practice. However, these errors can be made small, as we discuss in § 4. The practical accuracy of energy conservation in the standard JOREK tokamak model, which is what the new model presented in § 2 reduces to in the tokamak limit, is also considered in Hoelzl et al. (2021).

The necessity of neglecting perpendicular current in the resistive term of the first equation (2.12) when using the reduced MHD model can be clarified by deriving a magnetic field equation for the case when it is not neglected. Let $\varPsi _1$ and $\varPsi _2$ satisfy the equations

(3.9)\begin{equation} \left.\begin{gathered} \left[\psi_v,\frac{\partial \varPsi_1}{\partial t}\right] = \left[\frac{[\varPsi,\varPhi]-\partial^{{\parallel}}\varPhi}{B_v},\psi_v\right] + \frac{1}{B_v} \boldsymbol{\nabla}\boldsymbol{\cdot}(\eta\boldsymbol{\nabla}\psi_v\times\boldsymbol{j}^\parallel), \\ \left[\psi_v,\frac{\partial \varPsi_2}{\partial t}\right] = \frac{1}{B_v}\boldsymbol{\nabla}\boldsymbol{\cdot}(\eta\boldsymbol{\nabla}\psi_v\times\boldsymbol{j}^\perp), \end{gathered}\right\} \end{equation}

then $\varPsi = \varPsi _1 + \varPsi _2$ satisfies the original first equation in (2.12), without the neglect of perpendicular current. The first of the two equations above is equivalent to ${\partial \boldsymbol {B}_1/\partial t = -\boldsymbol {\nabla }\times \boldsymbol {E}_1}$, where $\boldsymbol {B}_1 = \boldsymbol {\nabla }\chi + \boldsymbol {\nabla }\varPsi _1\times \boldsymbol {\nabla }\chi$ and $\boldsymbol {E}_1 = -\boldsymbol {v}\times \boldsymbol {B} + \eta \boldsymbol {j}^\parallel$. Dotting with $\boldsymbol {B} = \boldsymbol {B}_1 + \boldsymbol {B}_2$ and adding $\boldsymbol {B}\boldsymbol {\cdot }\partial \boldsymbol {B}_2/\partial t + \boldsymbol {\nabla }\boldsymbol {\cdot }(\boldsymbol {E}_2\times \boldsymbol {B})$ to both sides, where $\boldsymbol {B}_2 = \boldsymbol {\nabla }\varPsi _2\times \boldsymbol {\nabla }\chi$ and $\boldsymbol {E}_2 = \eta \boldsymbol {j}^\perp$, we get

(3.10)\begin{equation} \frac{1}{2}\frac{\partial B^2}{\partial t} + \boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{E}\times\boldsymbol{B}) = \mu_0[(\boldsymbol{v}\times\boldsymbol{B})\boldsymbol{\cdot}\boldsymbol{j} - \eta j^2] + \boldsymbol{B}\boldsymbol{\cdot}\frac{\partial \boldsymbol{B}_2}{\partial t} + \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\times\boldsymbol{E}_2. \end{equation}

The evolution equation for $\varPsi _2$ can be written as $(\partial \boldsymbol {B}_2/\partial t)^{\psi _v} = -(\boldsymbol {\nabla }\times \boldsymbol {E}_2)^{\psi _v}$, and so the last two terms on the right-hand side of the above equation will only cancel partially, leaving behind the non-conservative terms $B_{\beta _v}\partial B_2^{\beta _v}/\partial t + B_{\beta _v}(\boldsymbol {\nabla }\times \boldsymbol {E}_2)^{\beta _v} + (\boldsymbol {\nabla }\times \boldsymbol {E}_2)^\chi$.

Finally, we show why the projection operators used in our original model allowed for non-physical generation of kinetic energy. To begin with, the velocity field can be written in a covariant form:

(3.11)\begin{equation} \boldsymbol{v} = B_v^2\boldsymbol{\nabla}^{{\perp}}\tilde{\varPhi}\times\boldsymbol{e}_\chi + \widetilde{v_{{\parallel}}}\boldsymbol{\nabla}\chi + B_v^2\boldsymbol{e}_\chi\times(\boldsymbol{\nabla}\tilde{\zeta}\times\boldsymbol{\nabla}\chi), \end{equation}

as opposed to the contravariant form (2.2). The fact that any vector field can be written in this way is easy to show using the same arguments as in Nikulsin et al. (2019), except using the projection operators (2.6) instead of the operators (2.4). Note that the projection operators (2.6) are orthogonal to all but one term in the covariant velocity ansatz (3.11):

(3.12)\begin{equation} \left. \begin{gathered} \boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{\nabla} \times\left(\frac{1}{B_v^2}\boldsymbol{v}\right) ={-}\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{\nabla}^{{\perp}}\tilde{\varPhi}), \\ \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{v} = B^\chi\widetilde{v_{{\parallel}}}, \\ \boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{\nabla} \times\left(\frac{1}{B_v^2}\boldsymbol{\nabla}\chi\times\boldsymbol{v}\right) = \boldsymbol{\nabla}\boldsymbol{\cdot}(B_v^2\boldsymbol{\nabla}^{{\perp}}\tilde{\zeta}). \end{gathered}\right\} \end{equation}

Now, applying the first projection operator in (2.4) to a vector field $\boldsymbol {Q}$ and then multiplying by a test function $\varPhi ^*$ and integrating, we get

(3.13)\begin{equation} \int_V \varPhi^*\boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{\nabla}\times[\boldsymbol{\nabla} \chi\times(\boldsymbol{e}_\chi\times\boldsymbol{Q})]\,\textrm{d}V ={-}\int_V B_v^2(\boldsymbol{\nabla}^{{\perp}}\varPhi^*\times\boldsymbol{e}_\chi)\boldsymbol{\cdot}\boldsymbol{Q}\,\textrm{d}V, \end{equation}

where we again used integration by parts to get the right-hand side. Similarly, for the third projection operator in (2.4), we can write

(3.14)\begin{equation} \int_V \zeta^*\boldsymbol{\nabla}\boldsymbol{\cdot}[B_v^2\boldsymbol{\nabla}\chi \times(\boldsymbol{e}_\chi\times\boldsymbol{Q})]\,\textrm{d}V = \int_V B_v^2[\boldsymbol{e}_\chi\times(\boldsymbol{\nabla}\zeta^*\times\boldsymbol{\nabla}\chi)] \boldsymbol{\cdot}\boldsymbol{Q}\,\textrm{d}V, \end{equation}

and for the second projection operator, we simply multiply by $v_{\parallel }^*$ and integrate, without doing any transformations:

(3.15)\begin{equation} \int_V v_{{\parallel}}^*\boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{Q}\,\textrm{d}V = \int_V v_{{\parallel}}^*\boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{Q}\,\textrm{d}V. \end{equation}

In order to get $\int _V \rho \boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {Q}\,\textrm {d}V = 0$ when we set $\boldsymbol {Q} \equiv \partial \boldsymbol {v}/\partial t + (\boldsymbol {v}\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {v} + \boldsymbol {v} P/\rho - \boldsymbol {j}\times \boldsymbol {B}/\rho + \boldsymbol {\nabla } p/\rho = 0$, we need $\varPhi ^* = -\tilde {\varPhi }$, $v_{\parallel }^* = \widetilde {v_{\parallel }}$, $\zeta ^* = \tilde {\zeta }$ and $\rho = \mathrm {const}$. The first three conditions pose no problem in the full MHD case, but if we try to set $\zeta = 0$, it does not imply $\tilde {\zeta } = 0$, and thus when we drop the evolution equation for $\zeta$, (3.14) is no longer satisfied, and so the kinetic energy equation (3.6) cannot be satisfied. In addition, because the vector momentum equation was divided by density before the projection operators were applied, unless the density is held constant, the kinetic energy equation (3.6) also cannot be satisfied even when $\zeta \neq 0$. Thus, both the reduction and spatially varying density can lead to non-physical generation of kinetic energy.

4. Numerical examples

In this section, we consider several simulations with the JOREK code for a test case of a tearing mode in a circular-cross-section tokamak with an aspect ratio of 10 (figure 1). We compare the reduced model without parallel flow from the original set of stellarator-capable models, as presented in Nikulsin et al. (2019), with the standard JOREK tokamak model without parallel flow (Hoelzl et al. 2021). Note that in the tokamak limit $\chi = F_0\phi$, the latter is equivalent to the the model without parallel flow from the set of new stellarator-capable models, as introduced in § 2. In addition, we point out that the standard tokamak reduced MHD model in JOREK is known to accurately reproduce tearing modes when compared with full MHD (Haverkort et al. 2016; Pamela et al. 2020), thus comparing with the standard tokamak model is sufficient.

Figure 1. A flux-aligned grid used for simulating the tearing mode (a), and the $n=1$ Fourier mode of $\psi = F_0\varPsi$ (JOREK units) in the standard tokamak model at $t=50\,000$ Alfvén times (b).

In the test cases considered below, $F_0 = 10\ \mathrm {T}\ \textrm {m}$. The initial conditions were set up by solving the Grad–Shafranov equation with $FF'(\psi _n) = 1.173(1-\psi _n)$ in units of $\mathrm {T}^2\ \textrm {m}$, where $\psi _n = (\psi - \psi _{\textrm {axis}})/(\psi _{\textrm {edge}}-\psi _{\textrm {axis}})$, and $p(\psi _n)=0$. In the tokamak limit, the $\psi$ in the Grad–Shafranov equation is related to the $\varPsi$ in the magnetic field ansatz (2.2) by $\psi = F_0\varPsi$. The density was initialized to be constant at $3.346\times 10^{-7}\ \mathrm {kg}\ \mathrm {m}^{-3}$, corresponding to $10^{20}$ deuterium ions per cubic metre, and $\varPhi$ was initialized to zero. A viscosity of $\mu = 5.159\times 10^{-8}\ \mathrm {kg}\ (\mathrm {m}\ \mathrm {s})^{-1}$ was used in all simulations, which corresponds to $10^{-5}$ in JOREK units for $\mu _\perp ^t$ in the standard tokamak model.² When using the original stellarator model, the kinematic viscosity was set to be $\nu = \mu /\rho _0 = 0.1542\ \mathrm {m}^2\ \mathrm {s}^{-1}$ ($10^{-7}$ in JOREK units). Finally, since JOREK discretizes the toroidal direction by approximating all solutions as toroidal Fourier series, in all simulations shown in this section, the Fourier series are truncated after the first mode for the sake of simplicity, keeping only the axisymmetric ($n=0$) and the $n=1$ terms, unless noted otherwise. Finite elements are used for discretization in the poloidal plane. Temporal discretization is done via the Crank–Nicolson scheme, i.e. an average of the forward and backward Euler schemes. Given equations of the type

(4.1)\begin{equation} \sum_i a_i(\boldsymbol{U})\frac{\partial A_i(\boldsymbol{U})}{\partial t} = B(\boldsymbol{U}), \end{equation}

where $\boldsymbol {U}$ is the set of all unknowns in the model and $a_i$, $A_i$ and $B$ are the expressions that comprise the equations and can involve spatial derivatives, the Crank-Nicolson scheme gives

(4.2)\begin{equation} \sum_i \left[\frac{1}{2}a_i(\boldsymbol{U}^{n+1}) + \frac{1}{2}a_i(\boldsymbol{U}~^n)\right][A_i(\boldsymbol{U}^{n+1}) - A_i(\boldsymbol{U}^n)] = \Delta t\left[\frac{1}{2}B(\boldsymbol{U}^{n+1}) + \frac{1}{2}B(\boldsymbol{U}^n)\right], \end{equation}

where $n$ in the superscript refers to the already calculated values of the model variables in the current time step and $n+1$ refers to the next time step, yet to be calculated. In addition, a temporal linearization is done around the current time step, leading to the actual numerical scheme that is implemented in the code:

(4.3)\begin{equation} \sum_i a_i(\boldsymbol{U}^n)\left.\frac{\partial A_i}{\partial\boldsymbol{U}}\right|_{\boldsymbol{U}^n}\boldsymbol{\cdot}\delta\boldsymbol{U}^n = \Delta t\left[B(\boldsymbol{U}^n) + \frac{1}{2}\left.\frac{\partial B}{\partial\boldsymbol{U}}\right|_{\boldsymbol{U}^n}\boldsymbol{\cdot}\delta\boldsymbol{U}^n\right], \end{equation}

where $\delta \boldsymbol {U}^n = \boldsymbol {U}^{n+1} - \boldsymbol {U}^n$. For more details on the numerical implementation, see Czarny & Huysmans (2008), Huysmans & Czarny (2007) and Hoelzl et al. (2021).

To begin with, we tested the original stellarator model (Nikulsin et al. 2019) in the linear regime by calculating the tearing mode growth rates at different resistivities and comparing them with the growth rates obtained from the standard tokamak model (Hoelzl et al. 2021). In both cases, we did a spatial and temporal resolution scan at each resistivity to ensure that the growth rate that we recorded was converged. We then repeated the simulations with Fourier modes up to and including $n=4$ with the same spatial and temporal resolutions at which the growth rates converged for the $n=0,1$ simulations. The results of the $n=0,1$ simulations are shown in figure 2(a). Figure 2(b) shows the growth rates of the nonlinearly driven $n=2$ modes in the simulations with $n=0,\ldots ,4$; these modes are not inherently unstable, but are destabilized by the unstable $n=1$ mode. As can be seen, both models show decent agreement. At each resistivity, the $n=2$ growth rates are roughly twice the $n=1$ growth rates, as expected. In both models at higher resistivities, the $n=3$ and $n=4$ modes are destabilized shortly before the onset of the nonlinear regime, and so the corresponding growth rates (not shown here) do not plateau, but rather peak at values roughly three and four times the value of the $n=1$ growth rate and then decline rapidly as nonlinear saturation is reached. All of the growth rates shown in figure 2 were calculated for $\beta =0$, and the pressure was not evolved. In the case of the the original stellarator model, this amounts to not using the full energy conservation equation (third equation in (2.1)). For clarity, we explicitly list the equations solved in table 1.

Figure 2. Tearing mode growth rates at various plasma resistivities in the original stellarator model and the standard tokamak model. (a) The growth rates for the $n=1$ toroidal mode in simulations with $n=0,1$. (b) The growth rates for the $n=2$ toroidal mode in simulations with $n=0,\ldots ,4$. The growth rates of the $n=2$ modes are double those of the $n=1$ modes, indicating that the $n=2$ modes are not naturally unstable but excited by nonlinear mode coupling.

Table 1. The zero-$\beta$ forms of the two reduced MHD models that are compared in this section, namely the original stellarator model from Nikulsin et al. (2019) and the standard tokamak model without field-aligned flow from Franck et al. (2015). Note that the new set of equations presented in § 2 reduces to the standard tokamak model when $\zeta = v_{\parallel } = \varOmega = 0$ and $\chi = F_0\phi$, except for the term containing $P$ (eighth term on the right-hand side of (2.9)). In the cases considered, $P = \boldsymbol {\nabla }\boldsymbol {\cdot }(D_\perp \nabla _\perp \rho )$. When using the original stellarator model, we set $\chi = F_0\phi$ and $\psi _v = R$, where $R$ is the distance from the central axis of symmetry; in addition a subscript $\chi$ means a dot product of the corresponding vector with $\boldsymbol {e}_\chi = \boldsymbol {B}/B_v^2$, and $\boldsymbol {f}_b = -FF'\boldsymbol {\nabla }\varPsi |_{t=0}/R^2$ is the force balancing term (see Appendix A). In order for the initial condition to be a true equilibrium, we have introduced $\boldsymbol {j}_d = \boldsymbol {j} - \boldsymbol {j}_0$, where $\boldsymbol {j}_0$ is the current at $t=0$.

While the original stellarator model (Nikulsin et al. 2019) performs decently in the linear regime and at $\beta =0$, two problems can arise after nonlinear saturation. First, as discussed in the previous section, non-physical kinetic energy can be generated. The rate at which it is generated is negligible in the linear regime, but can become significant as saturation is reached. This can be seen in figure 3, where we consider the tearing mode simulation with a resistivity of $10^{-5}$ JOREK units ($1.9382\times 10^{-5}\ \Omega \ \textrm {m}$). In figure 3, $-\textrm {d}E/\textrm {d}t$ is plotted and compared with the physical energy loss rate. We define $E = \int _V\mathcal {E}\,\textrm {d}V$, where $\mathcal {E}$ is given in (2.1); this is the integrated total energy of the system at a given time. The physical energy loss rate is defined as the sum of the energy fluxes across the plasma boundary and the volume integral of energy sinks due to resistive and viscous dissipation (conversion to internal energy is not accounted for since pressure is not evolved, so dissipated energy is lost). If the inward energy fluxes are greater than outward fluxes plus sinks, the physical loss rate is negative, i.e. energy is gained. As can be seen in figure 3(a), the difference between $-\textrm {d}E/\textrm {d}t$ and the physical loss rate is negligible in the linear regime (until approximately 20 ms) and grows rapidly with the onset of saturation, due to the uncontrollable increase in kinetic energy (figure 4a). The simulation is stopped at $\sim$27 ms, as it would crash shortly thereafter if allowed to continue. In order to allow the simulation to continue after saturation, we can introduce an artificial dissipation. Figure 3(b) shows $-\textrm {d}E/\textrm {d}t$ and physical energy loss rate for a simulation with artificial dissipation in the form of a hyperviscosity term (with $\nu _h = 1.08\times 10^{-3}\ \mathrm {m}^4\ \mathrm {s}^{-1}$, corresponding to $7\times 10^{-10}$ in JOREK units) in the evolution equation for $\varPhi$. In this case, the physical energy loss rate does not include the hyperviscous dissipation. As can be seen, with artificial dissipation, the energy behaves much more reasonably, though it is still not conserved, hence the mismatch between the total rate of energy change and physical energy loss. The reason for such good behaviour is that the non-conservation comes mostly from the kinetic energy, which is kept from exploding by the hyperviscosity term. In addition, the hyperviscosity prevents the plasma from responding too quickly to numerical errors in the magnetic field, thus stabilizing it against numerical instabilities that arise from the $\varPsi$ evolution equation (discussed below).

Figure 3. The negative rate of change of total energy compared with the physical energy loss rate in (a) the original stellarator model without artificial dissipation, (b) the original stellarator model with artificial dissipation, (c) the standard tokamak model/new stellarator model and (d) the original stellarator model with (2.18) replacing the $\varPsi$ equation in table 1.

Figure 4. The negative rate of change of total energy $-\textrm {d}E/\textrm {d}t$ from figure 3 plotted alongside the negative rates of change of the magnetic and kinetic energies, where $\textrm {d}E/\textrm {d}t = \textrm {d}E_{\textrm {mag}}/\textrm {d}t + \textrm {d}E_{\textrm {kin}}/\textrm {d}t$. (a) The negative rates of change in the original stellarator model without artificial dissipation. (b) The same in the case with artificial dissipation.

One can introduce a finite pressure into the original stellarator model by using (2.15) without affecting the energy conservation by much, since all of the error comes from the kinetic energy, which is usually not dominant in fusion-relevant situations. However, we find that the original stellarator model cannot correctly reproduce the growth rates when $\beta$ becomes non-negligible, likely because the pressure terms in the original stellarator model (Nikulsin et al. 2019) (not shown in table 1) contain either $\boldsymbol {\nabla }\rho$ or $\partial ^{\parallel } p$, both of which are zero initially, and so they make the pressure terms small compared to the Lorentz force terms.

Finally, for the sake of comparison, we show the same kind of graph for a simulation with the standard tokamak model/new stellarator model without parallel flow (figure 3c), both of which are equivalent in the tokamak limit. Note the small discrepancy between $-\textrm {d}E/\textrm {d}t$ and the physical energy loss rate in the last plot; unlike the discrepancy in figure 3(b), this discrepancy is purely numerical. As the arguments of § 3 continue to apply even after the exact solutions are replaced by finite-element approximations, the discrepancy is not due to poloidal discretization, and is independent of the resolution of the finite elements, which we have confirmed by running simulations with higher resolution. Instead, the two sources of error are: (a) the temporal discretization, where a small enough time step is used in practice to diminish the deviation, and (b) the truncation of the toroidal Fourier series, where a sufficiently large toroidal resolution is used in practice to avoid this effect. We have carried out further simulations (not shown here), where we have confirmed that by decreasing the time step size and increasing the number of Fourier modes kept, one can make the discrepancy arbitrarily small.

The second problem is that, while looking benign, the $\varPsi$ evolution equation can often produce ill-conditioned time stepping matrices, resulting in numerical instabilities. In fact, the $\varPsi$ equation is what destabilizes the kinetic energy, causing it to explode in figure 4(a). If we replace the $\varPsi$ evolution equation in the original stellarator model (see table 1) by (2.18), we can run the simulations well into the post-saturation regime without needing to introduce hyperviscosity, as can be seen in figure 3(d), albeit with a higher energy conservation error than when hyperviscosity is present. On the other hand, we can replace the $\psi$ evolution equation in the standard tokamak model by the corresponding equation from the original stellarator model (table 1), thus testing the $\varPsi$ equation from the original stellarator model in a setting where energy is conserved in the continuous limit. In this case (figure 5a), energy conservation issues begin around the same time that the original model without hyperviscosity would have crashed (figure 3a). However, instead of exploding, the energy rapidly decreases, and so this simulation can run slightly longer before crashing. Figure 5(b) shows that the kinetic energy plays no part in this numerical instability; all of the energy conservation error comes from the buildup of numerical errors in the magnetic energy. As can be seen, while the $\varPsi$ evolution equation in the original stellarator model in table 1 is equivalent to equation (2.18) in the analytical sense, numerically it behaves very differently. Finally, we verified that replacing the $\psi$ evolution equation in all models considered in this section does not have any significant impact on the growth rates.

Figure 5. A simulation of a tearing mode with the standard tokamak model when the $\psi$ evolution equation is replaced by the corresponding equation from the original stellarator model. (a) A comparison of the negative rate of change of total energy $-\textrm {d}E/\textrm {d}t$ with the physical loss rate. (b) The negative rates of change of magnetic and kinetic energy alongside $-\textrm {d}E/\textrm {d}t$. The inset in (a) zooms in on times between 25 and 35 ms.

To conclude this section, we note that we had initially intended to use the full energy conservation equation, as shown in (2.1), to evolve pressure and ensure energy conservation. However, using the full energy conservation equation would make the internal energy into a reservoir, drawing energy from it when there is a non-physical gain of kinetic energy, or depositing energy into it if there is a non-physical loss of kinetic energy. This shifts the energy conservation error to pressure while ensuring total energy conservation. While this may be acceptable when the internal energy is much larger than the kinetic energy, using the full energy conservation equation would thus disallow low-$\beta$ simulations. It is much better to use equation (2.15) with the original stellarator model instead of the full energy conservation equation, although, as noted above, the original stellarator model is unable to accurately reproduce growth rates when $\beta$ is not negligible.

5. Approximate conservation of momentum

Reduced MHD models usually do not conserve linear momentum exactly, except for some special cases. Locally, this can be seen from the fact that there are only two velocity variables, $\varPhi$ and $v_{\parallel }$, in reduced MHD, so the three components of the full MHD momentum equation cannot be satisfied simultaneously. However, different reduced MHD models can have different amounts of error in the momentum conservation. In the remainder of this section, we first compare momentum conservation properties of the original model presented in Nikulsin et al. (2019) with those of the new model derived in § 2, and then show how global momentum behaves in the new model using some numerical examples. When discussing momentum conservation, we do not exclude exchange of momentum with the walls, such as in the case of a vertical displacement event. Exchange of momentum is a physical process, and, while the total momentum of the simulated system can change, no momentum is created or destroyed. In this section, we are concerned with the non-physical generation of momentum within the plasma due to approximations made in reduced MHD.

5.1. General local momentum conservation properties of the reduced models

We begin by reviewing the momentum conservation properties of the original reduced model, which were discussed in detail in Nikulsin et al. (2019). The action of the first projection operator (2.4) on a vector $\boldsymbol {Q}$ can be written as $\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla }\times [\boldsymbol {\nabla }\chi \times ( \boldsymbol {e}_\chi \times \boldsymbol {Q})] = \boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla }\times (\boldsymbol {e}_\chi Q^\chi - \boldsymbol {Q})$. Thus, if $\boldsymbol {Q}$ is the full MHD momentum equation (2.1), the action of the projection operator corresponds to dropping the two vector components of a vorticity-like equation, $\boldsymbol {\nabla }\times (\boldsymbol {e}_\chi Q^\chi - \boldsymbol {Q})$, that are perpendicular to $\boldsymbol {\nabla }\chi$. If all three components of this equation were satisfied simultaneously (which, in general, is not possible as noted before), then the original full MHD momentum equation would also be satisfied and momentum would be conserved exactly. We can estimate the magnitude of momentum conservation error by considering the components of the vorticity-like equation perpendicular to $\boldsymbol {\nabla }\chi$. This equation can be written as

(5.1)\begin{align} &\frac{\partial \hbox {$\boldsymbol{\rlap{-} \omega}$}}{\partial t} + \mu_0 v_\parallel\frac{\partial \boldsymbol{j}}{\partial t} + \boldsymbol{\nabla} v_\parallel{\times}\frac{\partial \boldsymbol{B}}{\partial t} - \frac{\mu_0\boldsymbol{j}}{B_v}\partial^{{\parallel}} v^2 - \boldsymbol{\nabla}\left(\frac{\partial^{{\parallel}} v^2}{B_v}\right)\times\boldsymbol{B} + \boldsymbol{\nabla}\times(\boldsymbol{\omega}\times\boldsymbol{v})\nonumber\\ &\quad- \frac{\mu_0\,\boldsymbol{j}}{B_v^2}\boldsymbol{\nabla}\chi\boldsymbol{\cdot}(\boldsymbol{\omega}\times\boldsymbol{v}) -\boldsymbol{\nabla}\left[\frac{\boldsymbol{\nabla}\chi\boldsymbol{\cdot}(\boldsymbol{\omega} \times\boldsymbol{v})}{B_v^2}\right]\times\boldsymbol{B} + \hbox{ $\boldsymbol{\rlap{-} \omega}$}\,\frac{P}{\rho} + \boldsymbol{\nabla}\left(\frac{P}{\rho}\right)\times\hbox{$\boldsymbol{\rlap{-} v}$} = \frac{\boldsymbol{B}}{\rho^2}(\kern1.5pt\boldsymbol{j}\boldsymbol{\cdot}\boldsymbol{\nabla}\rho)\nonumber\\ &\quad- \frac{\boldsymbol{j}}{\rho^2}(\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\rho) \!+\! \frac{1}{\rho}(\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{j} \!-\! \frac{1}{\rho}(\kern1.5pt\boldsymbol{j}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{B} \!-\! \frac{\boldsymbol{j}}{\rho B_v^2}\boldsymbol{\nabla}\chi\boldsymbol{\cdot}(\kern1.5pt\boldsymbol{j}\times\boldsymbol{B}) \!-\! \boldsymbol{\nabla}\left[\frac{\boldsymbol{\nabla}\chi\boldsymbol{\cdot}(\kern1.5pt\boldsymbol{j}\times\boldsymbol{B})}{\rho B_v^2}\right]\times\boldsymbol{B}\nonumber\\ &\quad+ \frac{1}{\rho^2}\boldsymbol{\nabla}\rho\times\boldsymbol{\nabla} p + \frac{\boldsymbol{j}}{\rho B_v}\partial^{{\parallel}} p + \boldsymbol{\nabla}\left(\frac{\partial^{{\parallel}} p}{\rho B_v}\right)\times\boldsymbol{B}, \end{align}

where the reduced velocity $\hbox{$\boldsymbol{\rlap{-} v}$} = -\boldsymbol {\nabla }\chi \times (\boldsymbol {e}_\chi \times \boldsymbol {v}) = \boldsymbol {\nabla }\varPhi \times \boldsymbol {\nabla }\chi /B_v^2$ and the reduced vorticity $\hbox{$\boldsymbol{\rlap{-} \omega}$} = \boldsymbol {\nabla }\times \hbox{$\boldsymbol{\rlap{-} v}$}$ were introduced, and the viscosity term is not considered. If the components of this equation perpendicular to $\boldsymbol {\nabla }\chi$ are identically zero, then there is no approximation in the reduction as nothing is being neglected, and momentum is still conserved. The most general case in which these components are zero is the following:

(5.2)\begin{equation} \partial^{{\parallel}} u = 0, \quad u\in\{g^{ik}, \varPhi, \varPsi, v_{{\parallel}}, p, \rho, P\}, \end{equation}

where $g^{ik}$ are the components of the metric tensor of the vacuum field-aligned coordinate system. As can be shown by a simple calculation, in this case both $\hbox{$\boldsymbol{\rlap{-} \omega}$}$ and $\boldsymbol {j}$ will be directed strictly along $\boldsymbol {\nabla }\chi$. If we allow $g^{ik}$, $\varPhi$ or $\varPsi$ to vary along $\boldsymbol {\nabla }\chi$, the same calculation will show that $\hbox{$\boldsymbol{\rlap{-} \omega}$}$ has non-zero perpendicular components. This will cause $\partial \hbox{$\boldsymbol{\rlap{-} \omega}$}/\partial t$ to have non-zero components perpendicular to $\boldsymbol {\nabla }\chi$, which cannot be cancelled by any other terms since there are no more time derivatives involving $\varPhi$ in the equation. Similarly, if we let any of the other quantities vary along $\boldsymbol {\nabla }\chi$, the last term and the seventh term on the right-hand side (pressure), the first term and the seventh term on the right-hand side (density), the last term on the left-hand side ($P$) and the third term on the left-hand side ($v_{\parallel }$) will have non-zero perpendicular components, which will not be cancelled by any other terms. If the conditions (5.2) are met, then only the sixth and eighth terms on the left-hand side have non-zero perpendicular components. As can be shown by a simple expansion of the eighth term:

(5.3)\begin{equation} [\boldsymbol{\nabla}\times(\boldsymbol{\omega}\times\boldsymbol{v})]^\perp{-} \left[\boldsymbol{\nabla}\left[\frac{\boldsymbol{\nabla}\chi\boldsymbol{\cdot}(\boldsymbol{\omega} \times\boldsymbol{v})}{B_v^2}\right]\times\boldsymbol{B}\right]^\perp \equiv 0. \end{equation}

Now we consider the new reduced model. For the first projection operator (2.6), the vorticity-like equation analogous to (5.1) can be written as

(5.4)\begin{align} &\boldsymbol{\nabla}\left(\frac{\rho}{B_v^2}\right)\!\times\!\frac{\partial \hbox{$\boldsymbol{\rlap{-} v}$}}{\partial t} \!+\! \frac{\rho}{B_v^2}\frac{\partial \hbox{$\boldsymbol{\rlap{-} \omega}$}}{\partial t} \!+\! \boldsymbol{\nabla}\left(\frac{\rho}{B_v^2}\frac{\partial v_{{\parallel}}}{\partial t}\right)\times\boldsymbol{B} \!+\! \frac{\mu_0\rho}{B_v^2}\frac{\partial v_{{\parallel}}}{\partial t}\boldsymbol{j} \!+\! \boldsymbol{\nabla}\left(\frac{\rho v_{{\parallel}}}{B_v^2}\right)\times\frac{\partial \boldsymbol{B}}{\partial t} + \frac{\mu_0\rho v_{{\parallel}}}{B_v^2}\frac{\partial \boldsymbol{j}}{\partial t}\nonumber\\ &\quad + \frac{1}{2}\boldsymbol{\nabla}\left(\frac{\rho}{B_v^2}\right)\times\boldsymbol{\nabla} v^2 + \boldsymbol{\nabla}\times\left(\frac{\rho\boldsymbol{\omega}\times\boldsymbol{v}}{B_v^2}\right) + \boldsymbol{\nabla}\left(\frac{P}{B_v^2}\right)\times\hbox{$\boldsymbol{\rlap{-} v}$} + \frac{P}{B_v^2}\hbox{$\boldsymbol{\rlap{-} \omega}$} + \boldsymbol{\nabla}\left(\frac{Pv_{{\parallel}}}{B_v^2}\right)\times\boldsymbol{B} \nonumber\\ &\quad + \frac{\mu_0 Pv_{{\parallel}}}{B_v^2}\boldsymbol{j} = \frac{2\boldsymbol{B}}{B_v^3}\boldsymbol{j}\boldsymbol{\cdot}\boldsymbol{\nabla} B_v - \frac{2\boldsymbol{j}}{B_v^3}\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B_v + \frac{1}{B_v^2}(\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{j} - \frac{1}{B_v^2}(\kern1.5pt\boldsymbol{j}\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{B} + \frac{2}{B_v^3}\boldsymbol{\nabla} B_v\times\boldsymbol{\nabla} p. \end{align}

Equation (2.9) is just the contravariant $\chi$ component of this equation. In the above equation, the third and eleventh terms on the left-hand side will have perpendicular components even when condition (5.2) is met; namely $\boldsymbol {\nabla }^{\perp }(\rho B_v^{-2}\partial v _{\parallel }/\partial t)\times \boldsymbol {\nabla }\chi$ from the third term and $\boldsymbol {\nabla }^{\perp }(Pv_{\parallel } B_v^{-2})\times \boldsymbol {\nabla }\chi$ from the eleventh term. In addition, while the eighth term on the left-hand side of (5.4) is analogous to the sixth term on the left-hand side of (5.1), there is no term in (5.4) analogous to the eighth term of (5.1) to cancel the non-zero perpendicular components of its eighth term. Thus, a more stringent condition is required to ensure that the perpendicular components of equation (5.4) are identically zero and momentum is exactly conserved in the new reduced MHD model. Namely, in addition to the conditions (5.2), we must also have $v_{\parallel } = 0$. This additional condition causes the third and eleventh terms to vanish, and the eighth term becomes

(5.5)\begin{equation} \boldsymbol{\nabla}\times\left(\frac{\rho\hbox{$\boldsymbol{\rlap{-} \omega}$}\times\hbox{$\boldsymbol{\rlap{-} v}$}}{B_v^2}\right) =\boldsymbol{\nabla}\times\left(\frac{\rho\hbox{${\rlap{-} \omega}$}_\chi\boldsymbol{\nabla}^{{\perp}}\varPhi}{B_v^2}\right) =\boldsymbol{\nabla}\left(\frac{\rho\hbox {${\rlap{-} \omega}$}_\chi}{B_v^2}\right)\times\boldsymbol{\nabla}^{{\perp}}\varPhi - \frac{\rho\hbox{${\rlap{-} \omega}$}_\chi}{B_v^2}\boldsymbol{\nabla}\left(\frac{\partial^{{\parallel}}\varPhi}{B_v}\right)\times\boldsymbol{\nabla}\chi. \end{equation}

Clearly, if condition (5.2) is met, the eighth term, as shown above, will not have any perpendicular components.

As can be seen, the original reduced model (Nikulsin et al. 2019) has better momentum conservation properties than the new reduced model (§ 2). This was the reason why we initially attempted to work with the original model. However, as energy conservation is generally more important than momentum conservation, we will use the new model for future simulations.

5.2. Global momentum conservation error (numerical examples)

As a test case, we use a simple ballooning mode in an X-point plasma in a tokamak with an aspect ratio of 3 (figure 6). In the simulations shown here only the Fourier modes with $0 \leqslant n \leqslant 6$ were included. The simulation parameters were as follows: $F_0 = 3\ \mathrm {T}\ \mathrm {m}$, a time step of 3 Alfvén times, a resistivity of $\eta = 3.8764\times 10^{-6}\ \Omega \ \mathrm {m}$ ($2\times 10^{-6}$ JOREK units) and a viscosity of $\mu = 2.293\times 10^{-7}\ \mathrm {kg}\ (\mathrm {m}\ \mathrm {s})^{-1}$, corresponding to $\mu _\perp ^t = 4\times 10^{-6}$ JOREK units. The initial conditions were set up by solving the Grad–Shafranov equation with

(5.6)\begin{equation} FF'(\psi_n) = \frac{1}{2}\left[1.6(1 - \psi_n) - 0.43\cosh^{{-}2}\frac{\psi_n - 0.9}{0.07}\right]\left(1 - \tanh\frac{\psi_n - 1}{0.03}\right), \end{equation}

$T(\psi _n) = 0.015(1 - 0.66\psi _n)[1 - \tanh ((\psi _n - 0.94)/0.08)]/2 + 3\times 10^{-4}$ and $\rho (\psi _n)=[1-\tanh ((\psi _n - 0.94)/0.08)]/2 + 0.01$, where $FF'$ is in units of $\mathrm {T}^2\ \mathrm {m}$, $T$ is in $(10^{20}k_B\mu _0)$ K and $\rho$ is in $3.346\times 10^{-7}\ \mathrm {kg}\ \mathrm {m}^{-3}$. In addition, a hyperresistivity of $\eta _h = 5.8146\times 10^{-10}\ \Omega \ \mathrm {m}^{-2}$ ($3\times 10^{-10}$ JOREK units) was used in the simulations. We compare the new reduced model with $v_{\parallel } = 0$ with the new reduced model with $v_{\parallel }\neq 0$. Since in the tokamak limit these models match the standard JOREK tokamak models, we used the tokamak models for the actual simulations.

Figure 6. A flux-aligned grid used for simulating the ballooning mode, shown here with reduced resolution for clarity (a), and the sum of $n>0$ Fourier mode of $F_0\varPsi$ (JOREK units) in the standard tokamak model at 250 time steps, at $t=641.1$ Alfvén times (b).

To quantify the momentum conservation error, we calculate the total linear momentum in the Cartesian $x$ and $y$ directions for the entire plasma. Physically, this momentum should be zero throughout the simulation; however, due to the inaccuracies discussed in the previous subsection, a non-zero momentum tends to appear. Here, we choose the $x$ direction as the $R$ direction when $\phi =0$ and the $y$ direction as the $R$ direction when $\phi ={\rm \pi} /2$, with the $z$ axis being the axis of symmetry of the torus. We do not consider momentum in the $z$ direction since the global $z$ momentum is conserved in the tokamak limit, as can be seen by setting $\varPhi ^* = \ln R$ in (3.1) and recalling that $B_v = F_0/R$ in the tokamak case.

As can be seen in figure 7, after about 1.5 ms, the momentum conservation error in the model with $v_{\parallel }\neq 0$ is worse by more than an order of magnitude than that in the model with $v_{\parallel } = 0$. In both cases, the instability has reached nonlinear saturation, and at that point the momenta no longer grow in amplitude, but tend to oscillate around zero with the amplitude of the oscillations remaining at the same order of magnitude.

Figure 7. Total linear momentum in the Cartesian $x$ and $y$ directions as a function of time. Momentum evolution is shown for simulations using the reduced model with $v_{\parallel }=0$ (a) and the reduced model with $v_{\parallel }\neq 0$ (b).

6. Conclusion

In continuation of our previous work (Nikulsin et al. 2019), we have implemented the stellarator-capable reduced MHD model derived in our previous paper and tested it in the tokamak limit. We find that, even in the simple case of a tearing mode in a circular high-aspect-ratio tokamak, energy conservation can only be approximately satisfied, due to the possibility of non-physical gain or loss of kinetic energy. Our original plan, as presented in Nikulsin et al. (2019), was to use the full energy conservation equation when evolving pressure, thus ensuring energy conservation. However, that would require the internal energy to be much larger than the kinetic energy, thus excluding low-$\beta$ situations. Even without exact energy conservation, meaningful linear results, such as growth rates, could be obtained without any further modifications; however, the lack of energy conservation required artificial dissipation to be introduced to prevent a crash after nonlinear saturation was reached due to non-physical kinetic energy buildup. Introducing artificial dissipation in the form of a hyperviscosity term kept the kinetic energies in check and allowed the simulations to proceed into the post-saturation regime. Although exact energy conservation could not be achieved, the error in energy conservation remained small because the kinetic energy is small compared to the total energy, and the magnetic energy evolved correctly. Thus, energy should be approximately conserved whenever the kinetic energy is small, which is usually the case in fusion-relevant situations. We also found that the kinetic energy blowup was triggered by numerical errors from the magnetic field. Replacing the magnetic field equation with a different form that is analytically equivalent to the original, but more numerically stable alleviates the need for hyperviscosity, although the energy conservation error in this case is larger than when hyperviscosity is present.

To remedy the lack of energy conservation in our original model, we used a new set of projection operators to obtain scalar momentum equations from the full MHD vector momentum equation. With these new projection operators, the reduced MHD model matches the standard JOREK reduced MHD models for tokamaks (Franck et al. 2015; Hoelzl et al. 2021) in the tokamak limit, i.e. when $\chi = F_0\phi$, where $\boldsymbol {\nabla }\chi$ is the vacuum magnetic field and $\phi$ is the toroidal angle. When using the new model, energy is conserved almost exactly, with a small amount of error due to temporal discretization and the use of a truncated Fourier decomposition in the toroidal direction. However, the new model also has more error in the momentum conservation than the original model whenever $v_{\parallel }\neq 0$.