Axisymmetric equilibria with central current reversal caused by non-parallel plasma flow

TL;DR

This paper derives analytic solutions for axisymmetric plasma equilibria with arbitrary flow, showing that non-parallel flow components can induce central current reversal while maintaining nested magnetic surfaces.

Contribution

It introduces a generalized Grad-Shafranov equation solution that accounts for non-parallel plasma flow and demonstrates how such flow can cause central current reversal.

Findings

Non-parallel flow can produce equilibria with central current reversal.

Solutions include nested and non-nested magnetic surfaces.

Pressure profiles remain monotonic despite current reversal.

Abstract

We obtain analytic solutions of a generalised Grad-Shafranov equation describing steady states with incompressible plasma flow of arbitrary direction, toroidal current reversal and either nested or non-nested magnetic surfaces. It turns out that the component of the flow velocity non parallel to the magnetic field can result in normal equilibria with central current-reversal, i.e. equilibria with nested magnetic surfaces and monotonically varying pressure profiles.

Figures

Tables

Table 1. Table 1: Values of the parameters in ansatz ( 9 - 11 ) for the equilibria of Figs. 2 - 4 .

Parameters	${\alpha }_1$	${\beta }_1$	${\gamma }_1$	$\alpha$
Fig. 2	5.3	1.3	1.3	5.0
Fig. 3	1.5	1.5	0	3.5
Fig. 4	7.5	7.0	0	6.75

Equations

$$(1-M_p^2){\Delta }^{\star }\psi -\frac{1}{2}(M_p^2{)}^{\prime }|\nabla \psi {|}^2+\frac{1}{2}{\left(\frac{X^2}{1-M_p^2}\right)}^{\prime }$$

$$+{\mu }_0{R}^2{P}_s^{\prime }+{\mu }_0\frac{R^4}{2}{\left[\frac{\varrho ({\Phi }^{\prime }{)}^2}{1-M_p^2}\right]}^{\prime }=0$$

$$P=P_s(\psi )-\varrho \left[\frac{v^2}{2}-\frac{R^2({\Phi }^{\prime }{)}^2}{1-M_p^2}\right]$$

$$\mathbf{v}=\frac{M_p}{\sqrt{\varrho }}\mathbf{B}-R{\Phi }^{\prime }{\mathbf{e}}_{\varphi }$$

$$u(\psi )={\int }_0^{\psi }{\left[1-M_p^2(f)\right]}^{1/2}df$$

$${\Delta }^{\ast }u+\frac{1}{2}\frac{d{I}^2}{d\psi }+{\mu }_0{R}^2\frac{d{P}_s}{du}+\frac{{\mu }_0}{2}{R}^4\frac{dG}{du}=0$$

$$[{\partial }_{\rho \rho }-(1/\rho ){\partial }_{\rho }+{\partial }_{\zeta \zeta }]\tilde{u}$$

$$P(kPa)=160{\tilde{P}}_s(\tilde{u})-80\left[\frac{{\tilde{X}}^2}{1-M_p^2}\right]\frac{M_p^2}{{\rho }^2(1-M_p^2)}-$$

$$J_{tor}(MA/m^2)=0.057\frac{1}{\sqrt{1-M_p^2}}\frac{{\tilde{\Delta }}^{\ast }\tilde{u}}{\rho }+$$

$$\frac{1}{2}\frac{d}{du}\left[\frac{X^2}{1-M_p^2}\right]={\alpha }_1+{\alpha }^{2}u$$

$$\frac{d{P}_s}{du}=-{\beta }_1+{\alpha }^{4}u$$

$$\frac{dG}{du}=-{\gamma }_1$$

$$M_p^2=M_0^2{u}^n$$

$$M_p^2=M_0^2{u}^m(u_a-u{)}^n$$

$$[\frac{1}{\rho }{\partial }_{\rho }(\rho {\partial }_{\rho })+{\partial }_{\zeta \zeta }]A+\left({\alpha }^2-\frac{1}{{\rho }^2}\right)A+{\alpha }^4{\rho }^{2}A=$$

$$=-\frac{{\alpha }_1}{\rho }+{\beta }_1\rho +\frac{1}{4}{\gamma }_1{\rho }^3$$

$$A{|}_b=0$$

$$[\frac{1}{\rho }{\partial }_{\rho }(\rho {\partial }_{\rho })+{\partial }_{\zeta \zeta }]A+\left({\lambda }^2-\frac{1}{{\rho }^2}\right)A+{\lambda }^4{\rho }^{2}A=0$$

$$A(1-ϵ,\zeta )=A(1+ϵ,\zeta )=0$$

$$A(\rho ,-\kappa ϵ)=A(\rho ,\kappa ϵ)=0$$

$$\frac{d^{2}Z}{d{\zeta }^2}+{\nu }^{2}Z=0$$

$$Z(-\kappa ϵ)=Z(\kappa ϵ)=0$$

$$Z(\zeta )=cos({\nu }_l\zeta ),{\nu }_l=(l+1/2)\pi /\kappa ϵ,(l=0,1,2,...)$$

$${\rho }^2\frac{d^{2}U}{d{\rho }^2}+\rho \frac{dU}{d\rho }+({\mu }_n^2{\rho }^2+{\lambda }_{nl}^4{\rho }^4-1)U=0$$

$$U(1-ϵ)=U(1+ϵ)=0$$

$$U_{n,l}(\rho )=c_{n,l}{J}_1(\rho ;{\mu }_n,{\nu }_l)+Y_1(\rho ;{\mu }_n,{\nu }_l)$$

$$c_{n,l}=-\frac{Y_1((1-ϵ);{\mu }_n,{\nu }_l)}{J_1((1-ϵ);{\mu }_n,{\nu }_l)}$$

$$-J_1((1+ϵ);{\mu }_n,{\nu }_l)\frac{Y_1((1-ϵ);{\mu }_n,{\nu }_l)}{J_1((1-ϵ);{\mu }_n,{\nu }_l)}+Y_1((1+ϵ);{\mu }_n,{\nu }_l)=0$$

$$A(\rho ,\zeta )=\sum _{n=1}^{N_n}\sum _{l=0}^{N_l}{A}_{nl}{U}_{n,l}(\rho )cos({\nu }_l\zeta )$$

$$-\frac{{\alpha }_1}{\rho }+{\beta }_1\rho +\frac{1}{4}{\gamma }_1{\rho }^3=\sum _{n=1}^{10}\sum _{l=0}^4{A}_{nl}\left[({\alpha }^2-{\lambda }_{n,l}^2)+({\alpha }^4-{\lambda }_{n,l}^4){\rho }^2\right]U_{n,l}(\rho )cos({\nu }_l\zeta )$$

$$U_1:=\sum _{n=0}^{\infty }{u}_n{\rho }^n$$

$$u_n=-\frac{1}{(n^2-1)}[{\mu }^2{u}_{n-2}+{\lambda }^4{u}_{n-4}],(n\ge 5)$$

$$U_2=U_1(\rho )ln(\rho )+\frac{b_{-1}}{\rho }+\sum _{n=0}^{\infty }{b}_n{\rho }^n$$

$$(n^2-1)b_n=-2n{u}_n-{\mu }^2{b}_{n-2}-{\lambda }^4{b}_{n-4},(n\ge 5)$$

Full text

1 Introduction

Advanced confinement in tokamaks is related to Internal Transport Barriers (ITBs) of energy and particle transport (e.g. see the review paper [1]). Experimental evidence indicates that reversed magnetic shear and sheared flow play a role in the establishment of ITBs; specifically according to experimental results of JET [2] and DIII-D [3], on the one hand the reversed magnetic helps in triggering the ITBs development while on the other hand the sheared rotation has an impact on the subsequent growth and allows the formation of strong ITBs. In certain tokamak equilibria with ITBs known as current holes, the central current is nearly zero (e.g. see the review paper [4]). According to current hole experiments in JET [5] and JT-60U [6] the core current density is clamped at zero, indicating the existence of a physical mechanism which prevents it from becoming negative. These observations motivated a number of theoretical studies on the existence of equilibria with current reversal [7]-[23]. The conclusion is that central current reversal in normal equilibria, i.e. global equilibria with nested magnetic surfaces and monotonically varying pressure, can not exist [7]-[10]. Potential physical mechanisms for that prevention are the influence of a resistive kink magnetohydrodynamic instability [7] and reconnection [10]. However, current-reversal equilibria with non nested magnetic surfaces, i.e. having magnetic lobes, can exist [11]-[25] in consistence with experimental evidence in tokamaks [26]-[29].

The present study was partly motivated by the paper of Ref. [30], according to which toroidal rotation opens up the possibility of normal equilibria with current reversal in the vicinity of the magnetic axis. The purpose of our study is twofold: first to construct axisymmetric equilibria with current reversal and incompressible flow with either nested or non-nested magnetic surfaces, and second to examine the existence of normal stationary equilibria with current reversal near the magnetic axis. The construction is based on a generalized Grad-Shafranov (GGS) equation with incompressible flow (Eq. (1) below) involving five free surface-quantity terms. Specifically a linearised form of that equation is solved analytically.

The GGS equation is briefly reviewed in section 2. Solutions with current density reversal are derived for a toroidal plasma of arbitrary aspect ratio and rectangular poloidal cross section with ITER-like shaping in section 3. Part of the construction regarding two independent solutions of the respective homogeneous partial differential equation is given in Appendix. Section 4 summarizes the conclusions.

2 Generalized Grad-Shafranov equation

We consider the generalized Grad-Shafranov (GGS) equation with incompressible flow [31]-[33]:

[TABLE]

Here, the poloidal magnetic flux function $\psi(R,z)$ labels the magnetic surfaces, where ($R,\phi,z$) are cylindrical coordinates with $z$ corresponding to the axis of symmetry; $M_{p}(\psi)$ is the Mach function of the poloidal fluid velocity with respect to the poloidal Alfvén velocity; $X(\psi)$ relates to the toroidal magnetic field, $B_{\phi}=I/R$, through $I=X/(1-M^{2})$; $\Phi(\psi)$ is the electrostatic potential; for vanishing flow the surface function $P_{s}(\psi)$ coincides with the pressure; $B$ is the magnetic field modulus which can be expressed in terms of surface functions and $R$; $\Delta^{\star}=R^{2}\nabla\cdot(\nabla/R^{2})$; and the prime denotes derivatives with respect to $\psi$. Because of incompressibility the density $\varrho(\psi)$ is also a surface quantity and the Bernoulli equation for the pressure decouples from (1):

[TABLE]

where $v$ is the velocity modulus. Also, ${\mathbf{v}}$ can be decomposed in a component parallel and another non-parallel to ${\mathbf{B}}$ as

[TABLE]

The quantities $M_{p}(\psi)$, $X(\psi)$, $P_{s}(\psi)$, $\varrho(\psi)$ and $\Phi(\psi)$ are free functions. Derivation of (1) and (2) is provided in [31]-[33].

Eq. (1) can be simplified by the transformation

[TABLE]

under which (1) becomes

[TABLE]

where $G(u):=\varrho(\psi)(d\Phi(u)/du)^{2}$. Note that no quadratic term as $|{\bf\nabla}u|^{2}$ appears any more in (5). Once a solution of (5) is obtained, the equilibrium can be completely constructed with calculations in the $u$-space.

Instead of (5) we will employ the GGS in the following normalized form

[TABLE]

Here, $\rho:=R/R_{0}$, $\zeta:=z/R_{0}$, $\tilde{u}:=u/u_{0}$ (where $R_{0}$ is the major radius of the toroidal configuration), $\tilde{P}_{s}:=P_{s}/P_{0}$, $\tilde{X}:=X/X_{0}$, $\tilde{G}:=G/G_{0}$, with ITER-relevant data $P_{0}=1.6\times 10^{5}Pa\simeq 1.58Atm$, $R_{0}=6.2m$, $a=1.1m$ (where $a$ is the major radius), $u_{0}^{2}:=P_{0}\mu_{0}R_{0}^{4}$, $u_{0}\simeq 17.23Wb$, $X_{0}:=u_{0}/R_{0}=2.78Wb/m$, $I_{0}:=u_{0}/\mu_{0}R_{0}\simeq 2.21MA$, $G_{0}=4.16\times 10^{3}Kg/m^{3}sec^{2}$. Using the above parametric values the pressure and toroidal current density are given by

[TABLE]

[TABLE]

We hereafter drop tildes for simplicity and choose the following forms for the poloidal current density, pressure and non-parallel-flow function $G$:

[TABLE]

Eqs. (9)-(11) are integrated with zero integration constants for the last two equations and with integration constant $\alpha_{0}$ for the first one. Also, following [15] we introduce the new dependent variable $A=u/\rho$ and adopt the boundary condition $u|_{b}=A|_{b}=0$. Additionally, in connection with potential flow-caused central reversal of the toroidal current density the poloidal Mach function is assigned by

[TABLE]

where $M_{0}^{2}$ and $n$ relate to the maximum value of the poloidal Mach function on the magnetic axis and the shape of $M_{p}^{2}$, respectively. This function becomes maximum on axis and vanishes on the boundary; therefore it is peaked on-axis in association with a respective on-axis momentum source. Alternatively, we have employed the off-axis peaked poloidal Mach function

[TABLE]

where $u_{a}$ is the poloidal flux function on the magnetic axis and $M_{0}^{2}$ relates to the maximum of $M_{p}^{2}$; the parameters $m,n$ determine the shape of $M_{p}^{2}$ and the position of its maximum, i.e. this maximum corresponds to $mu_{a}/(m+n)$ with a measure of the profile width being $u_{a}/(m+n)$.

Furthermore, following [15], [34], we consider a boundary of rectangular cross section so that $\rho\in[1-\epsilon,1+\epsilon]$ and $\zeta\in[-\kappa\epsilon,\kappa\epsilon]$ where $\epsilon:=a/R_{0}\simeq 0.177$ is the inverse aspect ratio and $\kappa\simeq 1.86$ is the elongation. Thus the boundary value problem for the GGS equation assumes the form

[TABLE]

3 Solutions and equilibrium properties

We consider first the homogeneous boundary value problem and boundary conditions corresponding to Eqs. (14) and (15):

[TABLE]

and set $A(\rho,\zeta)=U(\rho)Z(\zeta)$. We thus have

[TABLE]

with solution

[TABLE]

Furthermore setting $\lambda_{nl}^{2}:=\nu_{l}^{2}+\mu_{n}^{2}$, $(n=1,2,...)$ we obtain for $U(\rho)$

[TABLE]

The differential equation (19) is solved in the Appendix and its two independent solutions are hereafter called $J_{1}(\rho;\mu_{n},\nu_{l})$ and $Y_{1}(\rho;\mu_{n},\nu_{l})$. The eigenfunctions of Eq. (19) are therefore of the form

[TABLE]

and implementing the boundary conditions (20) results in the eigenvalue equations

[TABLE]

The last of Eqs. (3) is solved numerically for $\mu_{1},...,\mu_{10}$ since for computational purposes we have retained terms of up to $n=10$. The dependence of the parameters $\mu_{n}$ as a function of $n=1,...,10$ (approximated to vary continuously) for $l=0,1,2$ is shown in Fig. 1.

To construct a particular solution of the inhomogeneous Eq. (14) we employ the method of expansion with respect to eigenfunctions [14, 15]; accordingly $A$ is of the form

[TABLE]

with the coefficients $A_{nl}$ to be determined. Here the sums can have an arbitrary number of terms ($N_{n}=\infty$, and $N_{l}=\infty$), but for computational purposes we have retained terms of up to $N_{n}=10$, $N_{l}=4$. Note that the functions $U_{n,l}(\rho)$ are not orthogonal in the interval $\rho\in[1-\epsilon,1+\epsilon]$. Also (23) satisfies the boundary condition automatically. Then Eq. (14) yields

[TABLE]

Thus, due to the non-orthogonality of the functions $U_{n,l}(\rho)$, we multiply Eq. (24) by $U_{n^{{}^{\prime}},l^{{}^{\prime}}}(\rho)cos(\nu_{l^{{}^{\prime}}}\zeta)$ and integrate in the relevant intervals. We obtain a linear$-(50\times 50)$ system for the expansion coefficients $A_{n,l}$ which is solved numerically.

Assigning proper values to the free parameters involved in ansatz (9)-(11), as given in Table 1, we have constructed equilibria with reversed current density. Three examples of them are shown in Figs. 2-4. The first equilibrium (Fig. 2) has nested magnetic surfaces and the others (Figs. 3, 4) non-nested ones; specifically the second equilibrium consists of a couple of magnetic lobes with magnetic axes orientated perpendicular to the axis of symmetry (Fig. 3) and the third one has three lobes with magnetic axes orientated parallel to the axis of symmetry (Fig. 4). In the second and third equilibria the plasma flows parallel to the magnetic field ($\gamma_{1}=0$). The single magnetic axis of the first equilibrium is located at $\rho\approx 1.002$ and its position remains nearly unaffected when the values of the flow parameters $M_{0}^{2}$ and $\gamma_{1}$ change.

We have examined the impact of the flow on the toroidal current density of equilibria with nested magnetic surfaces and found that the non parallel component of the flow can result in current reversal in the region of the magnetic axis. This is illustrated in Fig. 5 corresponding to the peaked on-axis choice (12) with $n=2$. The continuous curve therein represents a static equilibrium ($M_{0}^{2}=\gamma_{1}=0$) for which current reversal is not possible as already mentioned in section 1. Central current reversal occurs at a critical value of $(\gamma_{1})_{c}=3.0$, nearly irrespective of the value of the parallel-flow parameter $M_{0}^{2}$ [see Eqs. (3) and (12)]. This result is consistent with that of [30] where purely toroidal flow was considered. To check that the parallel flow has no appreciable impact on the central current reversal we also employed the peaked off-axis flow profile (13) with $m=n=2$. The results respective to those of Fig. 5 are shown in Fig. 6. Apparently the current density profiles therein are identical with those of Fig. 5 in the close vicinity of the magnetic axis. The off-axis flow, though, makes the region of current reversal narrower.

The pressure of the equilibrium of Fig. 2 monotonically decreases from the magnetic axis to the plasma boundary as can be seen in Fig. 7. In addition, the safety factor monotonically increasing from the magnetic axis to the plasma boundary is shown in Fig. 8.

Finally, it is noted that current reversal can occur in equilibria with non-nested magnetic surfaces even in the absence of flow [13]-[21].

4 Summary

Solving the GGS equation (6) analytically we have constructed equilibria with incompressible flow, reversed toroidal current density and either nested or non-nested magnetic surfaces. The non-parallel flow results in normal equilibria with nested magnetic surfaces with central current-reversal and monotonically varying pressure profiles. This reversal gets stronger as the non-parallel flow parameter $\gamma_{1}$ in Eq. (11) takes larger values and is nearly unaffected by the parallel component of the flow. Finally, it is noted that the stability of equilibria with non parallel flow remains a tough unsolved problem reflecting the absence of stability criteria even in the framework of hydrodynamics.

Aknowledgments

This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the National Programme for the Controlled Thermonuclear Fusion, Hellenic Republic. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix: Solution of the homogeneous Eq. (19)

Since Eq. (19), put in canonical form, has a regular singular point at $\rho=0$ it can be solved by the Frobenius method of infinite convergent series around this point, e.g. [35]. Here we will employ this method slightly modified. Accordingly we try a solution of Eq. (19) of the form

[TABLE]

We obtain $u_{0}=0$, $u_{1}$ is a freely specified constant which is hereafter, according to standard analysis [35], is set equal to $u_{1}=0.001$, $u_{2}=0$, $u_{3}=-\mu^{2}u_{1}/8$, $u_{4}=0$ and

[TABLE]

where $\lambda^{2}=\mu^{2}+\nu^{2}$. We hereafter call this solution $J_{1}(\rho;\mu,\nu)$. For sufficiently small values of the parameter $\mu$, $J_{1}$ has the usual oscillatory behaviour of the standard Bessel function. More specifically we have found, using extensive numerical tests, that for every $0<\rho_{1}<\rho_{2}$ there exist numbers $\epsilon_{1},\;\epsilon_{2}$ with $0<\epsilon_{2}<\epsilon_{1}$ and an integer $n_{0}$ such that $|U_{1}(\rho_{1})|<\epsilon_{1}$ and $|U_{1}(\rho_{2})|<\epsilon_{2}$ provided that the number of terms $n_{u}$ used in the series expansion of $U_{1}$ satisfies $n_{u}\geq n_{0}$. For computational purposes we have retained terms in the series expansion of up to $N_{s}=100$.

The second independent solution of Eq. (19)will be of the form

[TABLE]

and obtain $b_{0}=0$, $b_{-1}=-1/\mu^{2}$. Also, $b_{1}$ is a free coefficient which is set, according to the analysis in [35] equal to $b_{1}=-0.001$, $b_{2}=0$,

$8b_{3}=\left(\frac{3\mu^{2}}{8}+\frac{\lambda^{4}}{\mu^{2}}-b_{1}\mu^{2}\right)$, $b_{4}=0$ and

[TABLE]

We call this solution hereafter $Y_{1}(\rho;\mu,\nu)$.

Since the closest singular point of (19) to $\rho=0$ is infinity, solutions $J_{1}(\rho;\mu,\nu)$ and $Y_{1}(\rho;\mu,\nu)$ converge for any finite value of $\rho$.