Implemented equations

Implemented equations (normalised)

In GRILLIX, the equations are solved in normalised form. (Hyper)diffusion, particle/energy sources and diffusion in buffer zones close to the radial boundaries are added.

Normalisation

We apply the following normalisation: $$\begin{align*} \mathbf{B}= & B_0\hat{\mathbf{B}}, & n= & n_0\hat{n}, & T_{e}=&T_{e0}\hat{T}_e, & T_i=&T_{i0}\hat{T}_i, & \phi=&\frac{T_{e0}}{e}\hat{\phi}, & u_\parallel=&c_{s0}\hat{u}_\parallel, \\ v_\parallel=&c_{s0}\hat{v}_\parallel,& j_\parallel=&c_{s0}en_0\hat{j}_\parallel, & \mathbf{x}_\perp= &\rho_{s0}\hat{\mathbf{x}}_\perp, & x_\parallel = & R_0\hat{x}_\parallel, & t = & \frac{R_0}{c_{s0}}\hat{t}, & A_\parallel=&\beta_0B_0\rho_{s0}\hat{A}_\parallel, \end{align*}$$ with $c_{s0}:=\sqrt{T_{e0}/M_i}$, $\rho_{s0}:=c\sqrt{T_{e0}M_i}/(eB_0)$ and $\beta_0:=4\pi n_0T_{e0}/B_0^2$ the dynamical plasma beta. (Note that the classical plasma beta is $\beta=\frac{8\pi n_0\left(T_{e0}+T_{i0}\right)}{B_0^2}$.) We define additionaly as dimensionless quantity $\delta=R_0/\rho_{s0}$ and introduce the normalized quantities: $$\begin{align*} \hat{\mathbf{v}}_E = &\delta\frac{\hat{\mathbf{B}}}{\hat{B}^2}\times\hat{\nabla}\hat{\phi}, & \hat{\mathcal{C}}(f)=\delta\left[-\left(\hat{\nabla}\times\frac{\hat{\mathbf{B}}}{\hat{B}^2}\right)\cdot\hat{\nabla}f\right]. \end{align*}$$

In a cylindrical coordinate system $(\mathbf{e}_R,\mathbf{e}_\varphi,\mathbf{e}_Z)$ with $\hat{\mathbf{B}}\approx \hat{B}^\varphi\mathbf{e}_\varphi$, we can define the Poisson bracket (discretised in GRILLIX with the Arakawa scheme) as $$\begin{align} \left[g,f\right]_{\hat{R},\hat{Z}} = \partial_{\hat{R}}g\partial_{\hat{Z}}f - \partial_{\hat{Z}}g\partial_{\hat{R}}f \approx - \mathbf{b}_0 \cdot \left( \nabla g \times \nabla f \right), \end{align}$$ and the normalised operators become $$\begin{align*} &\hat{\mathbf{v}}_E\cdot\hat{\nabla}\hat{f}=-\frac{\delta}{\hat{B}_\mathrm{tor}}\left[ \hat{\phi}, \hat{f} \right]_{\hat{R},\hat{Z}}, \\ &\hat{\mathcal{C}}(f)=-\frac{\delta}{J}\left[\frac{1}{B^\varphi},f\right], \\ &\hat{\nabla}_\parallel f = \mathbf{b}\cdot\hat{\nabla} f = \hat{\nabla}_\parallel^{(0)} f + \delta\beta_0 \left[\frac{\hat{A}_1}{\hat{B}_\mathrm{tor}}, f\right]_{\hat{R},\hat{Z}}, \end{align*}$$ where $J=R$ the Jacobian of the cylindrical coordinate system and $\hat{B}_\mathrm{tor} = \hat{B}^\varphi J = \hat{B}_\varphi / J$ the toroidal magnetic field strength. $\hat{A}_1$ and $\hat{A}_\parallel$ share the same normalisation. $\hat{A}_1$ is computed according to magnetic shift. Under the assumption $\hat{B}_\mathrm{tor}=1/R$, the flutter terms can ensure the zero divergence of magnetic fluctuations according to Magnetic fluctuations, and the curvature operator reduces to $\hat{\mathcal{C}}(f)=-2\partial_Z f$, which is employed artificially for the curvature operator in slab and circular geometry where $\hat{B}_\mathrm{tor}=1$ otherwise.

Finally, for ion viscosity, we obtain $$\begin{align*} G=n_0T_{i0}\hat{G}, \quad \text{with} \quad \hat{G} = -\eta_{i0}\hat{T}_i^{5/2}\left[\frac{2}{B^{3/2}}\hat{\nabla}\cdot\left( \hat{u}_\parallel B^{3/2} \mathbf{b} \right)-\frac{\hat{\mathcal{C}}(\hat{\phi})}{2}-\frac{\hat{\mathcal{C}}(\hat{p_i})}{2\hat{n}}\right], \end{align*}$$ where $\eta_{i 0}=0.96\tau_{i0}c_{s0}/R_0$ a dimensionless parameter.

In the following we only use normalised quantities and drop the hat. In order to ensure its positivity we develop the logarithms of the density, electron and ion temperature, i.e. $\theta_n:=\log(n)$, $\xi_e:=\log(T_e)$ and $\xi_i:=\log(T_i)$.

Equations

The 6 dynamical equations (plus Ampere's law) solved in GRILLIX are $$\begin{align} \left(\frac{\partial}{\partial t} +\mathbf{v}_E\cdot\nabla+u_\parallel\nabla_\parallel\right)\theta_n+\nabla\cdot \left(\mathbf{b}u_\parallel\right)=&\mathcal{C}(\phi)-T_e\left[\mathcal{C}(\xi_e)+\mathcal{C}(\theta_n)\right]+\frac{1}{n}\nabla\cdot\left(\mathbf{b}j_\parallel\right) \notag\\ &+\frac{1}{n}\left[\mathcal{D}_n+S_n\right], \end{align}$$ $$\begin{align} \nabla\cdot\left[\frac{n}{B^2}\left(\frac{\partial}{\partial t}+\mathbf{v}_E\cdot\nabla+u_\parallel\nabla_\parallel\right)(\nabla_\perp\phi +\zeta\frac{\nabla_\perp p_i}{n}) \right]=&-p_e\left[\mathcal{C}(\theta_n)+\mathcal{C}(\xi_e)\right]-\zeta p_i\left[\mathcal{C}(\theta_n)+\mathcal{C}(\xi_i)\right] \notag \\ & +\nabla\cdot\left(\mathbf{b}j_\parallel\right)-\frac{1}{6}\zeta\mathcal{C}(G) \notag \\ & + \mathcal{D}_\Omega, \label{eq_final:vorticity} \end{align}$$

$$\begin{align} \left(\frac{\partial}{\partial t}+\mathbf{v}_E\cdot\nabla +u_\parallel\nabla_\parallel\right)u_\parallel=&-T_e\left[\nabla_\parallel\theta_n+\nabla_\parallel\xi_e\right] -\zeta T_i\left[\nabla_\parallel\theta_n+\nabla_\parallel\xi_i\right] \notag \\ & + \zeta T_i\mathcal{C}(u_\parallel) -\frac{2}{3n}\zeta B^{3/2}\nabla_\parallel \frac{G}{B^{3/2}} \notag \\ &+\mathcal{D}_u, \label{eq_final:upar} \end{align}$$

$$\begin{align} \beta_0\frac{\partial}{\partial t}A_\parallel+\mu\left[\frac{\partial}{\partial t}+\mathbf{v}_E\cdot\nabla+v_\parallel\nabla_\parallel\right]\left(\frac{j_\parallel}{n}\right)= & -\left(\frac{\eta_{\parallel 0}}{T_e^{3/2}}\right)j_\parallel+T_e\left[\nabla_\parallel\theta_n+(1+\color{blue}{0.71})\nabla_\parallel\xi_e\right]-\nabla_\parallel\phi \notag \\ &+\mathcal{D}_\Psi, \label{eq_final:Ohm} \end{align}$$

$$\begin{align} \left(\frac{\partial}{\partial t}+\mathbf{v}_E\cdot\nabla+v_\parallel\nabla_\parallel\right)\xi_e & = \frac{2}{3}\mathcal{C}(\phi)-T_e\left[\frac{2}{3}\mathcal{C}(\theta_n)+\frac{7}{3}\mathcal{C}(\xi_e)\right]-\frac{2}{3}\nabla\cdot\left(\mathbf{b}v_\parallel\right)+\frac{2}{3n}0.71\nabla\cdot\left(\mathbf{b}j_\parallel\right) \notag \\ & + \frac{2}{3n}\frac{\eta_{\parallel0}}{T_e^{5/2}}j_\parallel^2+\frac{2}{3nT_e}\nabla\cdot \left( \mathbf{b}q_\parallel^e \right) \notag \\ & - 2\nu_{e0}\mu \frac{n}{T_e^{3/2}}\left(1-\zeta \frac{T_i}{T_e}\right) \notag \\ & + \frac{1}{T_e}\left(\mathcal{D}_{T_e}+S_{et}\right), \label{eq_final:electrontemperature} \end{align}$$

$$\begin{align} \left(\frac{\partial}{\partial t}+\mathbf{v}_E\cdot\nabla+u_\parallel\nabla_\parallel\right)\xi_i=&\frac{2}{3}\mathcal{C}(\phi)-\frac{2}{3}T_e\left[\mathcal{C}(\theta_n)+\mathcal{C}(\xi_e)\right]+\frac{10}{6}\zeta T_i\mathcal{C}(\xi_i)-\frac{2}{3}\nabla\cdot\left(\mathbf{b}v_\parallel\right) \notag \\ &+\frac{2}{3nT_i}\nabla\cdot \left( \mathbf{b}q_\parallel^i \right) + \frac{2}{3n} j_\parallel \nabla_\parallel \theta_n \notag \\ & +\frac{2}{9nT_i^{7/2}\eta_{i0}}G^2 + 2\nu_{e0}\mu \frac{n}{T_e^{3/2}}\left(\frac{1}{\zeta} (T_e-T_i) - 1\right) \notag \\ &+\frac{1}{T_i}\left(\mathcal{D}_{T_i}+S_{it}\right), \label{eq_final:iontemperature} \end{align}$$

$$\begin{align} \nabla_\perp^2A_\parallel=-j_\parallel. \end{align}$$

For the electron and ion heat fluxes $q_\parallel^{e,i}$, different models are avilable (see Beyond Braginskii for the transcollisional extensions). The Braginskii expressions are $$\begin{align} q_\parallel^{e,i} = \chi_{\parallel 0}^{e,i} T_{e,i}^{7/2}\nabla_\parallel\xi_{e,i}. \end{align}$$

$S_n,S_{et},S_{it}$ are source functions driving the system with particles and energy.

The dimensionless physics parameters of the system are $$\begin{align} \delta = & \frac{R_0}{\rho_{s0}}, & \beta_0 = & \frac{4\pi n_0T_{e0}}{B_0^2}, & \mu = & \frac{m_e}{M_i}, \\ \nu_{e0} = & \left( \frac{c_{s0}}{R_0}\tau_{e0} \right)^{-1}, % \nu_{e0}: normalized collisionality / electron collision frequency & \eta_{\parallel0} = & 0.51\mu\nu_{e0}, & \chi_{\parallel 0}^e = & \color{blue}{3.16}\frac{c_{s0}}{R_0}\tau_{e0}\frac{M_i}{m_e}, \\ \zeta = & \frac{T_{i0}}{T_{e0}}, & \eta_{i0} = & 0.96\frac{c_{s0}}{R_0}\tau_{i0}, % \eta_{i0}: ion viscosity & \chi_{\parallel 0}^i = & 3.9\frac{c_{s0}}{R_0}\tau_{i0}\frac{T_{i0}}{T_{e0}}. \label{dimlessparams} \end{align}$$ $\tau_{e0}$ and $\tau_{i0}$ are the electron and ion collision times evaluated at reference $n_0$ and $T_0$, i.e $$\begin{align} \tau_{e0}= &\frac{3\sqrt{m_e}T_{e0}^{3/2}}{4\sqrt{2\pi}\lambda e^4Zn_0}=\frac{3.5\cdot10^4}{\lambda/10}\frac{T_{e0}^{3/2}[eV]}{Zn_0[cm^{-3}]}[s], \\ \tau_{i0}= &\frac{\tau_{e0}}{Z^2}\left(\frac{T_{i0}}{T_{e0}}\right)^{3/2}\sqrt{\frac{2M_i}{m_e}}. \end{align}$$ Note: not only collision times depend on ion charge state $Z$, or $Z_\mathrm{eff}$, but also blue marked coefficients in the thermal force (Ohm's law & $T_e$) and electron heat conductivity. For $\chi_{\parallel0}^e$ we have $\gamma_\chi \approx 13.58\times\frac{Z+0.21}{Z+4.2}$, and for thermal force $\gamma_{j_\parallel T_e} \approx 1.5\times\frac{Z+0.55}{Z+2.273}$ [see Chankin 2018]. The additional dependence in resistivity $\eta_{\parallel0}$ is weak and can be ignored [compare V M Zhdanov, Transport Processes in Multicomponent Plasma, eq.~(3.4.14)].

There are two auxiliary quantities: the "generalized vorticity" $$\begin{equation} \Omega = \nabla\cdot\ \left[ \frac{n}{B^2}\left(\nabla_\perp\phi+\zeta \frac{\nabla_\perp p_i}{n}\right)\right] \end{equation}$$ for the numerical integration of \eqref{eq_final:vorticity}, and $$\begin{equation} \Psi_m = \beta_0 A_\| + \mu \frac{j_\|}{n}=\beta_0 A_\| - \mu \frac{\nabla_\perp^2A_\|}{n} \end{equation}$$ for the numerical integration of \eqref{eq_final:Ohm}. Note that the dissipation term in Ohm's law \eqref{eq_final:Ohm} is also applied to the dynamical variable $\Psi_m$.

Numerical dissipation

Each equation can include dissipation in three forms. The first is via a perpendicular (in-plane) hyperdissipation operator, defined as: $$\begin{align} \mathcal{D}_\perp(f)=\nabla\cdot\left[\nu_\perp\nabla_\perp^{2N-1}f\right], \end{align}$$ where $N$ is the order of hyperdissipation. The perpendicular dissipation in buffer zones is given by: $$\begin{align} \mathcal{D}_{buff}(f)=\nabla\cdot\left[\nu_{buff}\nabla_\perp f\right], \end{align}$$ Different dissipation coefficients,$\nu_\perp$ and $\nu_{buff}$ can be selected for various quantities. Furthermore, both coefficients can exhibit spatial dependence, allowing for the configuration of mutually exclusive hyperdissipation and buffer dissipation regions.

Finally, parallel dissipation can be applied along the static background magnetic field: $$\begin{align} \mathcal{D}_{\parallel}(f)=\nabla\cdot\left[\nu_{\parallel}\mathbf{b}_0\nabla_\parallel^{(0)} f\right]. \end{align}$$

On density and vorticity all three dissipations are applied: $$\begin{align} \mathcal{D}_{n}=&\mathcal{D}_\perp(n) + \mathcal{D}_{buff}(n) +\mathcal{D}_{\parallel}(n), \\ \mathcal{D}_{\Omega}=&\mathcal{D}_\perp(\Omega) + \mathcal{D}_{buff}(\Omega) +\mathcal{D}_{\parallel}(\Omega). \end{align}$$

On parallel velocity and in Ohm's law, only perpendicular dissipation is applied: $$\begin{align} \mathcal{D}_{u}=&\mathcal{D}_\perp(u) + \mathcal{D}_{buff}(u), \\ \mathcal{D}_{\Psi}=&\mathcal{D}_\perp(\Psi) + \mathcal{D}_{buff}(\Psi). \end{align}$$

Due to the presence of physical parallel diffusion, only perpendicular dissipation is applied to the temperatures. Two distinct forms of dissipation can be chosen. The first form directly applies the dissipation to the temperature: $$\begin{align} \mathcal{D}_{T_{i,e}}=&\mathcal{D}_\perp(T_{i,e}) + \mathcal{D}_{buff}(T_{i,e}). \end{align}$$ In this form, the dissipation term acts as a sink (and may even act locally as a source) for thermal energy. Alternatively, a second form is available, where dissipation is applied to the pressure instead: $$\begin{align} \mathcal{D}_{p_{i,e}}=&\frac{1}{n}\left(\mathcal{D}_\perp(p_{i,e}) + \mathcal{D}_{buff}(p_{i,e})\right) - \frac{T_{e,i}}{n}\left(\mathcal{D}_\perp(n) + \mathcal{D}_{buff}(n)\right). \end{align}$$ This version ensures that the perpendicular dissipation does not act as a sink or source. Instead, the dissipated energy is absorbed by the background, giving rise to a dissipative flux.