Gyrokinetic internal kink simulations

Getting linear kink-type structure is easy using GYGLES, EUTERPE, ORB5
Target: coupling kink with another physics (turbulence, AEs, EPs)
Challenge: getting kink right ( $\delta B_{||}$ etc)
Nonlinear kink evolution

Examples

Internal kink in large aspect ratio tokamak

This case has been considered using GYGLES, EUTERPE, and ORB5 (see Section 4.2).

We consider an “ad-hoc” tokamak with the minor radius $r_{\rm a} = 1~{\rm m}$ , the major radius $R_0 = 10~{\rm m}$ , the magnetic field at the axis $B_0 = 1~{\rm T}$ , flat ion and electron temperatures $T_e = T_i$ defined by $L_{x} = 2 r_{\rm a} / \rho_{\rm s} = 360$ with $\rho_{\rm s} = \sqrt{m_i T_e}/(e B)$ the sound gyroradius. The safety factor is $q(s) = 0.8 (1 + s^2)$ , the flux surface label $s = \sqrt{\psi/\psi_{\rm a}}$ with $\psi_{\rm a}$ the poloidal flux at the plasma edge. The ambient plasma density profile $n_{0{\rm i}}(s) = n_{0{\rm e}}(s)$ is given by

n_{0 e} (s) = n_{0} \exp [- Δ_{n} κ_{n} \tanh (\frac{s - s_{0}}{Δ_{n}})]

$n_{0e}(s) = n_0 \exp\left[\,-\,\Delta_{\rm n} \kappa_{\rm n} \tanh\left(\frac{s-s_0}{\Delta_{\rm n}}\right)\right]$
with

κ_{n} = 3.0

$\kappa_{\rm n} = 3.0$ ,

s_{0} = 0.5

$s_0 = 0.5$ ,

Δ_{n} = 0.2

$\Delta_{\rm n} = 0.2$ , and

n_{0}

$n_0$ corresponding
to

β_{e} = μ_{0} n_{0 e} T_{e} / B^{2} = 0.0052

$\beta_{\rm e} = \mu_0 n_{0\rm e} T_{\rm e}/B^2 = 0.0052$ .

We use $N_{\rm e} = 64 \times 10^6$ electron markers, $N_{\rm i} = 16 \times 10^6$ ion markers, $N_s = 200$ radial grid points, $N_{\theta} = 16$ poloidal grid points, and $N_{\varphi} = 8$ toroidal grid points. The time step is $\omega_{c\rm i} \Delta t = 10$ . The Fourier filter includes the poloidal modes $m \in [-2, 2]$ and the toroidal mode $n = 1$ .

The equilibrium was not consistent (no Grad-Shafranov solved), $\delta B_{\|}=0$ was used.

Internal kink in a straight tokamak (GYGLES)

Cosidered for “usual” safety factor profiles in PoP2012

For “W7-X like” rotational transform profiles in NF2021

Internal kink (linear) in a tokamak (GYGLES)

Internal kink in nonlinear turbulent environment (tokamak, ORB5)

Temperature profiles

Electrostatic potential evolution including turbulence (poloidal plane)
Electrostatic potential evolution, toroidal mode number n = 1
Electrostatic potential evolution, toroidal spectrum

Reduced model for delta B_||

Recipe: replace on each occurance

Ω_{\nabla B} = \frac{b \times \nabla B}{B} ⟹ Ω_{κ} = (\nabla \times b)_{⊥}

$\begin{equation} \newcommand{\veps}{\varepsilon} \newcommand{\vc}[1]{{\bf#1}} \newcommand{\gav}[1]{\langle#1\rangle} \newcommand{\Bgav}[1]{\Big\langle#1\Big\rangle} \vc{\Omega}_{\nabla B} =\frac{{\bf b}\times\nabla B}{B} \Longrightarrow \vc{\Omega}_{\kappa} = (\nabla\times{\bf b})_{\perp} \end{equation}$
Proof via vorticilty equation (GK momentum equation) where

δ B_{‖}

$\delta B_{\|}$ appears as

K = δ {\dot{B}}_{‖} + Ω_{B} \cdot \nabla δ ψ, δ {\dot{A}}_{‖} = - \nabla_{‖} δ ψ

$\begin{equation} K = \delta \dot{B}_{\|} + \vc{\Omega}_B \cdot \nabla\delta\psi \ , \;\; \delta\dot{A}_{\|} = -\nabla_{\|}\delta\psi \end{equation}$
The following conditions are used:

Perpendicular force balance:

$δ B_{‖} = \frac{ξ_{⊥}}{B} \cdot \nabla P ⟺ δ {\dot{B}}_{‖} = \frac{\nabla ϕ \times B}{B^{3}} \cdot \nabla P$ $\delta B_{\|} = \frac{\vc{\xi}_{\perp}}{B}\cdot\nabla P \;\Longleftrightarrow\; \delta\dot{B}_{\|} = \frac{\nabla\phi\times\vc{B}}{B^3} \cdot\nabla P$
Ideal Ohm’s law

$E_{‖} = b \cdot \nabla ϕ + δ {\dot{A}}_{‖} = 0$ $E_{\|} = \vc{b}\cdot\nabla\phi + \delta\dot{A}_{\|} = 0$
MHD equilibrium

$κ = (b \cdot \nabla) b = \frac{μ_{0} \nabla P}{B^{2}} + \frac{\nabla B}{B}$ $\vc{\kappa} = (\vc{b}\cdot\nabla)\vc{b} = \frac{\mu_0\nabla P}{B^2} + \frac{\nabla B}{B}$

If these condition are satisfied, one can show

K = δ {\dot{B}}_{‖} + Ω_{B} \cdot \nabla ψ = Ω_{κ} \cdot \nabla ϕ

$K = \delta \dot{B}_{\|} + \vc{\Omega}_B \cdot \nabla \psi = \vc{\Omega}_{\kappa}\cdot\nabla\phi$

Unperturbed orbits

Original equation (Frieman-Chen, Littlejohn, Hahm, Brizard, etc):

\begin{array}{rcl} {\dot{R}}^{(0)} = v_{‖} b + \frac{1}{q_{s} B_{‖}^{*}} (m_{s} μ B \frac{b \times \nabla B}{B} + m_{s} v_{‖}^{2} \nabla \times b) = v_{‖} b + v_{d} \\ v_{d} \approx \frac{1}{q_{s} B} (m_{s} μ B \underset{Ω_{\nabla B}}{\underset{⏟}{\frac{b \times \nabla B}{B}}} + m_{s} v_{‖}^{2} \underset{Ω_{κ}}{\underset{⏟}{(\nabla \times b)_{⊥}}}), (\nabla \times b)_{⊥} = b \times (b \cdot \nabla) b \end{array}

$\begin{eqnarray} {\bf \dot{R}}^{(0)} = v_{\|} {\bf b} + \frac{1}{q_s B_{\|}^*} \left(m_s \mu B \frac{{\bf b}\times\nabla B}{B} + m_s v_{\|}^2 \nabla\times{\bf b} \right) = v_{\|} {\bf b} + {\bf v}_d \\ {\bf v}_d \approx \frac{1}{q_s B}\left(m_s\mu B \underbrace{\frac{{\bf b}\times\nabla B}{B}}_{\Omega_{\nabla B}} + m_s v_{\|}^2 \underbrace{(\nabla\times{\bf b})_{\perp}}_{\Omega_{\kappa}}\right) \ , \;\; (\nabla\times{\bf b})_{\perp} = {\bf b}\times ({\bf b}\cdot\nabla){\bf b} \end{eqnarray}$

Compressional force balance: replace $\Omega_{\nabla B}$ with $\Omega_{\kappa}$ , see [Graves2019]

{\dot{R}}^{(0)} = v_{‖} b + \frac{m_{s}}{q_{s} B} (μ B + v_{‖}^{2}) (\nabla \times b)_{⊥}

$\begin{equation} \dot{{\bf R}}^{(0)} = v_{\|} {\bf b} + \frac{m_s}{q_s B} (\mu B + v_{\|}^2) (\nabla\times{\bf b})_{\perp} \end{equation}$

No need to change unperturbed velocity equation

{\dot{v}}_{‖}^{(0)} = - μ (b + \frac{m_{s} v_{‖}}{q_{s} B_{‖}^{*}} \nabla \times b) \cdot \nabla B, μ = \frac{v_{⊥}^{2}}{2 B}

$\dot{v}_{\|}^{(0)} = \,-\, \mu \,\left( \vc{b} + \frac{m_s v_{\|}}{q_s B_{\|}^*} \nabla\times\vc{b} \right) \cdot \nabla B \ , \;\; \mu = \frac{v_{\perp}^2}{2 B}$

Perturbed equations of motion

No need to change perturbed part of particle orbits

{\dot{R}}^{(1)} = \frac{b}{B_{‖}^{*}} \times \nabla ⟨ ϕ - v_{‖} A_{‖}^{(s)} - v_{‖} A_{‖}^{(h)} ⟩ - \frac{q_{s}}{m_{s}} ⟨ A_{‖}^{(h)} ⟩ (b + \frac{m_{s} v_{‖}}{q_{s} B_{‖}^{*}} \nabla \times b)

$\dot{{\bf R}}^{(1)} = \frac{{\bf b}}{B_{\|}^*} \times \nabla \Bgav{\phi - v_{\|} A_{\|}^{\rm(s)} - v_{\|} A_{\|}^{\rm(h)}} - \frac{q_s}{m_s} \,\Bgav{A^{\rm(h)}_{\|}} \, \left( \vc{b} + \frac{m_s v_{\|}}{q_s B_{\|}^*} \nabla\times\vc{b} \right)$

Changes in energy equation needed. Original equation

\begin{array}{rcl} {\dot{ε}}^{(1)} & = & v_{‖} {\dot{v}}_{‖}^{(1)} + μ {\dot{R}}^{(1)} \cdot \nabla B = - \frac{1}{B} [μ B \frac{b \times \nabla B}{B} + v_{‖}^{2} (\nabla \times b)_{⊥}] \cdot \nabla ⟨ ϕ ⟩ + \\ \frac{v_{‖}}{B} [\frac{q_{s} B}{m_{s}} v_{‖} b + μ B \frac{b \times \nabla B}{B} + v_{‖}^{2} (\nabla \times b)_{⊥}] \cdot \nabla ⟨ A_{‖}^{(h)} ⟩ + \\ \frac{q_{s}}{m_{s}} μ B [\nabla \cdot b - \frac{m_{s} v_{‖}}{q_{s} B_{‖}^{*}} \frac{\nabla \times B}{B^{2}} \cdot \nabla B] ⟨ A_{‖}^{(h)} ⟩ \end{array}

$\begin{eqnarray} \dot{\veps}^{(1)} &=& v_{\|} \dot{v}_{\|}^{(1)} + \mu \dot{\vc{R}}^{(1)} \cdot \nabla B = \,-\, \frac{1}{B} \left[ \mu B \frac{\vc{b} \times \nabla B}{B} + v_{\|}^2 \, (\nabla \times \vc{b})_{\perp} \right] \cdot \nabla \gav{\phi} \,+\, \nonumber \\ &&{} \frac{v_{\|}}{B} \left[ \frac{q_s B}{m_s} v_{\|}\vc{b} + \mu B \,\frac{\vc{b} \times \nabla B}{B} + v_{\|}^2 \, (\nabla \times \vc{b})_{\perp} \right] \cdot \nabla \gav{A_{\|}^{\rm(h)}} + \\ &&{} \frac{q_s}{m_s} \mu B \left[ \nabla\cdot\vc{b} - \frac{m_s v_{\|}}{q_s B_{\|}^*} \frac{\nabla \times \vc{B}}{B^2} \cdot \nabla B \right] \Bgav{A_{\|}^{\rm(h)}} \end{eqnarray}$

transforms to

\begin{array}{rcl} {\dot{ε}}^{(1)} & = & - \frac{1}{B} (μ B + v_{‖}^{2}) (\nabla \times b)_{⊥} \cdot \nabla ⟨ ϕ ⟩ + \\ \frac{v_{‖}}{B} [\frac{q_{s} B}{m_{s}} v_{‖} b + (μ B + v_{‖}^{2}) (\nabla \times b)_{⊥}] \cdot \nabla ⟨ A_{‖}^{(h)} ⟩ + \\ \frac{q_{s}}{m_{s}} μ B [\nabla \cdot b - \frac{m_{s} v_{‖}}{q_{s} B_{‖}^{*}} \frac{\nabla \times B}{B^{2}} \cdot \nabla B] ⟨ A_{‖}^{(h)} ⟩ \end{array}

$\begin{eqnarray} \dot{\veps}^{(1)} &=& \,-\, \frac{1}{B} (\mu B + v_{\|}^2) \, (\nabla \times \vc{b})_{\perp} \cdot \nabla \gav{\phi} \,+\, \nonumber \\ &&{} \frac{v_{\|}}{B} \left[ \frac{q_s B}{m_s} v_{\|}\vc{b} + (\mu B + v_{\|}^2) \, (\nabla \times \vc{b})_{\perp} \right] \cdot \nabla \gav{A_{\|}^{\rm(h)}} + \\ &&{} \frac{q_s}{m_s} \mu B \left[ \nabla\cdot\vc{b} - \frac{m_s v_{\|}}{q_s B_{\|}^*} \frac{\nabla \times \vc{B}}{B^2} \cdot \nabla B \right] \Bgav{A_{\|}^{\rm(h)}} \end{eqnarray}$

Solving for delta B_|| (EUTERPE)

Equations of motion solved by EUTERPE (note that $\mu$ is the magnetic moment per mass):

\begin{array}{rcl} \dot{\vec{R}} = (v_{‖} - \frac{q}{m} ⟨ A_{‖} ⟩) {\vec{b}}^{*} + \frac{1}{B_{‖}^{*}} \vec{b} \times (\frac{m}{q} μ \nabla B + \nabla ⟨ Ψ ⟩ + \frac{m}{q} μ \nabla δ B_{‖}) \\ {\dot{v}}_{‖} = - \frac{q}{m} {\vec{b}}^{*} (\frac{m}{q} μ \nabla B + \nabla ⟨ Ψ ⟩ + \frac{m}{q} μ \nabla δ B_{‖}) \end{array}

$\begin{eqnarray} \newcommand{\vpar}{v_{\|}} \newcommand{\apar}{A_{\|}} &&{} \dot{\vec{R}} = (\vpar-\frac{q}{m} \langle \apar \rangle) \vec{b}^{*}+\frac{1}{B_{\|}^{*}} \vec{b} \times \left( \frac{m}{q} \mu \nabla B + \nabla \langle \Psi \rangle + \frac{m}{q} \mu \, \nabla \delta B_{\|} \right)\\ &&{} \dot{v}_{\|} = -\frac{q}{m} \vec{b}^{*} \left( \frac{m}{q} \mu \nabla B + \nabla \langle \Psi \rangle + \frac{m}{q} \mu \nabla \delta B_{\|} \right)\\ \end{eqnarray}$

closed by the pressure balance:

δ B_{‖} = - \frac{μ_{0}}{B} \sum_{s = i, e} m_{s} \int μ B δ f_{s} δ (\vec{R} - \vec{x}) d^{6} z

$\begin{equation*} \delta B_{\|} = -\frac{\mu_{0}}{B} \sum_{s=i,e} m_{s} \int \mu B \delta f_{s} \delta(\vec{R}-\vec{x}) {\rm d}^{6} z \end{equation*}$

Quasineutrality equations and parallel Ampere’s law are not modified. Definitions:

μ = \frac{v_{⊥}^{2}}{2 B}, {\vec{B}}^{*} = \vec{B} + \frac{m}{q} v_{‖} \nabla \times \vec{b}

$\begin{equation*} \mu=\frac{v_{\perp}^{2}}{2 B}, \quad \vec{B}^{*}=\vec{B}+ \frac{m}{q} \vpar \nabla \times \vec{b} \end{equation*}$

d^{6} z = B_{‖}^{*} d^{3} R d v_{‖} d μ d α, Ψ = Φ - v_{‖} A_{‖}, ⟨ Φ ⟩ = \frac{1}{2 π} \int Φ (\vec{R} + \vec{ρ}) d α

$\begin{equation*} {\rm d}^{6} z = B_{\|}^{*} \, {\rm d}^{3}R \, {\rm d} \vpar \, {\rm d} \mu \, {\rm d} \alpha, \quad \Psi=\Phi-\vpar \apar, \quad \langle \Phi \rangle = \frac{1}{2 \pi} \int \Phi(\vec{R}+\vec{\rho}) {\rm d} \alpha \end{equation*}$

Background current (driving instability) in GK equations

Which distribution function is appropriate? (shifted Maxwellian is normally used for electrons)
Shift is given by flux-surface averaged $\Bgav{\vc{b}\cdot\nabla\times\vc{B}}$

Pfirsch-Schlueter current vs. Shafranov shift

Cancellation of destabilization (PS current) and stabilization (Shafranov shift). Poloidal variation of ambient current has to be included.