Issues with the Boussinesq approximation in the vorticity equation

[yichenfu]

Overview

I'm working on 3-D hermres-3 simulations, and I found that the Boussinesq approximation might lead to some unphysical behavior in the simulations. We may need a non-Boussinesq option to benchmark the results, which is related to #56 .

Setup

Using an input file similar to the tcv-x21 case on dipole geometry. The grid in a poloidal cross section is something like:

The initial density $n_0(\psi)$ and temperature $T_0(\psi)$ are the follows. Initial random purterbations are added using mixmode.

With no source, hermes-3 is run in a profile-driven setting to see how the plasma evolves.

Issue # 1: profile stability depends on Nnorm

I found that with the same physical profile, if I set a different normalization density Nnorm, the stability of the profile could be different, as shown in the figure below. In general, if I use a larger Nnorm, the system tends to be more stable, and vice versa.

I suspect that the stability issue comes from the Boussinesq approximation in the vorticity equation, which is roughly:

$$\frac{\partial \omega}{\partial t}= -\nabla \cdot\left[\frac{m_{\mathrm{i}}}{2 B^2} \nabla_{\perp}\left(\mathbf{v}_E \cdot \nabla p_{\mathrm{i}}\right)+\frac{\omega}{2} \mathbf{v}_E +\frac{m_{\mathrm{i}} n_0}{2 B^2} \tilde{\nabla}_{\perp}^2 \phi\left(\mathbf{v}_E+\frac{\mathbf{b}}{n_0 B} \times \nabla p_{\mathrm{i}}\right)\right]+\nabla \cdot\left(\mathbf{J}_{\|}+\mathbf{J}_{\mathrm{d}}\right).$$

The Boussinesq approximation uses

$$\omega = \nabla \cdot \left[ \dfrac{m_\mathrm{i}}{B^2 q} (n \nabla_\perp \phi + \nabla_\perp p_\mathrm{i}) \right] \approx \nabla \cdot \left[ \dfrac{m_\mathrm{i} {\color{blue}n_0}}{B^2 q} \nabla_\perp \left( \phi + \dfrac{p_\mathrm{i}}{{\color{blue}n_0}} \right) \right],$$

which changes the magnitude the vorticity, therefore breaks the balance of the right-hand side terms in the vorticity equation, and moves a case across the stability boundary.

I have to say that my case is specifially designed near the stability boundary. Therefore, its stability is very susceptible to the Boussinesq approximation.

Issue # 2: background electrostatic potential also depends on Nnorm

I also notice that the background phi I got from the code also depends on Nnorm. For example, at the initial few steps, the vorticity is roughly zero. Inversion gives phi as

$$\omega = 0 \quad \Rightarrow \quad \phi = C - {p_\mathrm{i}}/{n_0}.$$

So, the initial phi is inversely proportional to Nnorm, as the following figure.

So, different reference densities lead to different equilibrium potentials, which may affect the shearing profile and change stability and transport.