Boundary conditions

Boundary conditions

Radial boundaries

Boundary conditions are required at the radially limiting flux surfaces $\rho_{min}$ and $\rho_{max}$.

Radial boundaries are somewhat artificial but have to be chosen consistently via an educated guess such that there is no flux of particles and energy through the radial boundary. Our typical set is $$\begin{align*} \partial_\rho n=&0, & \partial_\rho T_{e,i} =& 0, & \partial_\rho u_\parallel =& 0, & \Omega=\nabla\cdot\left(\frac{n}{B^2}\nabla_\perp\phi\right)=&0, &A_\parallel =& 0, & \partial_\rho j_\parallel=&0. \end{align*}$$

A special role is attributed to the potential. At the outer wall, it became standard to set $\left.\phi\right|_{\rho_{max}}=\Lambda T_e$, for consistency with sheath physics. Note however that this allows fluxes through the wall. If the wall is far enough away from the separatrix, and the plasma has enough time to flow out in parallel direction, this can work fine, but it can also become a problem.

In order to avoid $E\times B$ fluxes through the inner limiting flux surface, the potential must be a constant at the core, $\left.\phi\right|_{\rho_{min}}=const.$. Following GDB 2018, also in GRILLIX this constant is chosen as the zonal average $\left.\phi\right|_{\rho_{min}}=\left\langle\phi\right\rangle_{\rho_{min}}$, allowing the potential to float. This is called 'zonal Neumann' boundary conditions.

Additionally buffer zones are present close to the radial boundaries, where perpendicular diffusion is drastically increased, damping out turbulent fluctuations. Mathematically, this transforms hyperbolic (advective) PDE's into parabolic (diffusive) ones, for which boundary conditions are easier to set and are well defined.

Parallel boundary conditions

Theory on sheath boundary conditions can be found e.g. in "The Plasma Boundary of Magnetic Fusion Devices" by Peter Stangeby. If the magnetic field is strictly orthogonal to the boundary, one may use strictly sonic or allow supersonic sheaths, i.e. $$\begin{align*} u_\parallel= & \pm c_s = \pm\sqrt{(T_e + \gamma_B T_i)/M_i} \text{ (sonic)}, & u_\parallel \gtrless & \pm c_s \text{ (supersonic)}, \end{align*}$$ where the $'+'$ sign corresponds to the boundary where the magnetic field is directed towards the plate and the $'-'$ away from the target. According to the above references, the strictly sonic boundary condition only holds in the isothermal case. Usually, $\gamma_B = 1$ is employed.

In the case where the magnetic field hits the boundary under a glancing angle of incidence, which is the standard case in tokamaks and stellarators, the generalised Bohm boundary condition taking into account the $E\times B$ drift (Rozhansky 2001) has to be employed: the total particle velocity along the boundary normal vector $\mathbf{n}$ pointing outside of the domain (towards the boundary) has to be $$\begin{align*} u_\parallel \mathbf{b}\cdot\mathbf{n} + \mathbf{v}_E\cdot\mathbf{n} \geq c_s |\mathbf{b}\cdot\mathbf{n}|, \end{align*}$$ where again the sonic or the supersonic version can be used.

For the parallel current, isolating sheath boundary conditions are implemented, i.e. $$\begin{align*} & j_\parallel= 0 & \rightarrow & & \phi = \Lambda T_e \end{align*}$$ with $$\begin{align*} \Lambda = - \frac{1}{2} \left[ 2\pi\frac{m_e}{M_i} \left(1+\gamma_B\frac{T_i}{T_e}\right)\right] \approx 3. \end{align*}$$ We note that in Oliveira & Body 2022 it was found that $\nabla_\parallel j_\parallel = 0$ works somewhat better in some cases, but most of the time it does not seem to make a difference.

No parallel boundary condition on $A_\parallel$ is needed as no parallel operator acts on it, but it is set consistently with $j_\parallel=-\nabla_\perp^2A_\parallel$.

The electron temperature is set according to $$\begin{align*} \nabla_\parallel\log(T_e)=&-\frac{\tilde{\gamma}_e}{\chi_\parallel}nu_\parallel, \end{align*}$$ with $\tilde{\gamma}_e \in (1,3)$ the effective sheath heat transmission coefficient (only on the conducted heat flux, i.e. subtracting 2.5 from the convective heat flux). Strictly speaking, $\tilde{\gamma}_e$ is only constant if $j_\parallel=0$, or equivalently if $e\phi=\Lambda T_e$.

A similar option is implemented for $T_i$. However, it is not clear which coefficient $\tilde{\gamma}_i$ should be used. Typical choices are (0 - 0.1).

The insulating sheath boundary conditions run stable, but are not perfetly physically valid, because currents can actually flow through the divertor (and other walls). Describing this is challenging, since it might require the consideration of electric circuits outside the plasma volume. Local conducting sheath boundary conditions, e.g. from Loizu 2012, have been implemented and used in many codes. But in our experience, they are ill defined and ill behaved in tokamak geometry, particularly with complex divertors (they do work in slab). At least this is the case for strongly fluctuating turbulent dynamics, where such boundary conditions can even lead to spurious energy inflows.

It remains to discuss density and vorticity boundary conditions. Since they are supposed to flow out, in principle, no boundary condition is needed. However, in our experience (perhaps due to mimetic centered finite differences), in practise something must be set. The options are $$\begin{align*} \nabla_\parallel^K\log(n)=&0, &\nabla_\parallel^K\Omega=&0, \end{align*}$$ where $K=2$ (extrapolation) would be equivalent to no boundary conditions, but does not work in practise. For density, the only practical boundary condition is $K=1$ (homogeneous Neumann). For vorticity, both $K=0$ and $K=1$ work equally well.

We want to mention that boundary conditions are a complicated and rich topic by itself and are subject of current research.