2.1. Axisymmetric steady-state visco-resistive MHD: self-consistent system of equations

The framework used in this study is MHD. In more precise terms, building on the research initiated by Montgomery and his collaborators (Kamp, Montgomery & Bates Reference Kamp, Montgomery and Bates1998), we assume that the axisymmetric steady states of the plasma are governed by the incompressible visco-resistive MHD. This is consistent with the customary reconstruction of two-dimensional equilibria using the Grad–Shafranov equation, except that we do not assume the velocity field to be zero, and we have a self-consistent model, as we do not have free functions. Then, to describe a tokamak plasma, an essential aspect is to model the external drives involved in the system. One inherent drive in this magnetic confinement fusion device is the external magnetic field. Additionally, the need to wind the magnetic field lines and create a macroscopic poloidal component of the magnetic field requires a second forcing mechanism. Following previous references (Kamp & Montgomery Reference Kamp and Montgomery2003, Reference Kamp, Montgomery and Bates2004; Oueslati et al. Reference Oueslati, Bonnet, Minesi, Firpo and Salhi2019; Oueslati & Firpo Reference Oueslati and Firpo2020; Krupka & Firpo Reference Krupka and Firpo2024), we assume that the poloidal magnetic field component is generated by a toroidal electric field, which drives a toroidal current density.

Denoting by $B_{0}$ the value of the external magnetic field on the magnetic axis, by $\mu _{0}$ the vacuum permeability and by $\rho _{m0}$ the plasma mass density assumed to be constant, the Alfvén velocity is $v_{A0}=B_{0}/(\mu _{0}\rho _{m0})^{1/2}$ . In the remainder of this article, we shall work with dimensionless variables. Specifically, velocities are normalised with respect to the Alfvén speed, $v_{A0}$ , as is the field $\boldsymbol{B}/(\mu _{0}\rho _{m0})^{1/2}$ . Moreover, the space variables are also dimensionless. From the set of cylindrical polar coordinates $(r,\varphi ,z)$ and denoting by $R_{0}$ the tokamak major radius, we define $R=r/R_{0}$ and $Z=z/R_{0}$ . The dimensionless resistivity, $\eta$ , is the resistivity divided by $\mu _{0} R_{0} v_{A0}$ and the dimensionless viscosity, $\nu$ , is the viscosity divided by $R_{0} v_{A0}$ . The computation of visco-resistive axisymmetric steady states involves then solving the steady-state incompressible Navier–Stokes equations (2.1)–(2.2) along with the solenoidal condition (2.3), Ampère law (2.4) and Ohm’s law (2.5) on a tokamak poloidal plasma cross-section $\varOmega$ . The equations are

(2.1) \begin{align} &(\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla })\boldsymbol{v}=\boldsymbol{J}\,{\times}\, \boldsymbol{B}-\boldsymbol{\nabla }p+\nu \boldsymbol{\nabla} ^2\boldsymbol{v},\\[-10pt]\nonumber \end{align}

(2.2) \begin{align} &\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}=0,\\[-10pt]\nonumber \end{align}

(2.3) \begin{align} &\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{B}=0,\\[-10pt]\nonumber \end{align}

(2.4) \begin{align} &\boldsymbol{\nabla }\boldsymbol\,{\times}\,\boldsymbol{B}=\boldsymbol{J},\\[-10pt]\nonumber \end{align}

(2.5) \begin{align} &\boldsymbol{E}+\boldsymbol{v}\boldsymbol\,{\times}\,\boldsymbol{B}=\eta \boldsymbol{J}. \end{align}

With respect to the drives, the externally applied (vacuum) toroidal magnetic field is

(2.6) \begin{equation} \boldsymbol{B}_{\mathrm{ext}}=\frac {1}{R}\boldsymbol{i}_{\varphi }, \end{equation}

with $\boldsymbol{i}_{\varphi }$ a unit vector in the toroidal (azimuthal) direction. As for the electric field, its general expression is

(2.7) \begin{equation} \boldsymbol{E}=-\mathbf{\boldsymbol{\nabla }}\varPhi +V\mathbf{\boldsymbol{\nabla }}\left ( \frac {\varphi }{2\pi }\right ), \end{equation}

where $-\mathbf{\boldsymbol{\nabla }}\varPhi$ is a purely electrostatic contribution and where $V$ is the loop voltage. This second term serves to inductively generate the toroidal current, it is, therefore, non-steady state in nature and amounts to an external drive

(2.8) \begin{equation} \boldsymbol{E}_{\mathrm{ext}}=\frac {E_{0}}{R}\boldsymbol{i}_{\varphi },\end{equation}

with $E_{0}=(2\pi )^{-1}\partial \psi _{pol}/\partial t$ . In the present modelling, we consider the time variation of the poloidal magnetic flux $\psi _{\mathrm{pol}}$ , and thus $E_{0}$ , as constant and non-zero. This field satisfies ${\boldsymbol{\nabla }\,{\times}\, \boldsymbol{E}}_{\mathrm{ext}}=\mathbf{0}$ pointwise for $R\gt 0$ , but its global structure exhibits a non-vanishing toroidal circulation, consistent with a finite $V=\partial \psi _{\mathrm{pol}}/\partial t$ . Consequently, although $\boldsymbol{E}_{\mathrm{ext}}$ is locally curl free in the classical sense, it does not derive from a single-valued scalar potential, and Faraday’s law must be interpreted in the weak (distributional) sense. Thus, the present frame is strictly speaking not steady state as it is required that $\partial \psi _{\mathrm{pol}}/\partial t$ be non-zero, but it is time independent. Nevertheless, once the current drive is implemented (see § 4.1), we can put $E_{0}=0$ in the model equations and just have $\boldsymbol{E}=-\mathbf{\boldsymbol{\nabla }}\varPhi .$ Let us note, indeed, that several tokamak experiments have demonstrated scenarios where the externally applied electric field is curl free while maintaining a zero toroidal component (2.8). This occurs in fully non-inductive discharges driven by mechanisms such as lower hybrid current drive, electron cyclotron current drive and neutral beam current drive, where the electric field is electrostatic or wave driven, not induced by transformer action (see e.g. Sauter et al. (Reference Sauter2000); Litaudon et al. (Reference Litaudon2002)). The magnetic and electric fields in (2.1)–(2.5) are the sum of these external contributions and of the self-consistent plasma fields. This system of equations needs to be solved in the plasma cross-section $\varOmega$ with suitable boundary conditions. From a computational perspective, we solve the system of partial differential equations (PDEs) that we are now presenting.

2.2. Scalar PDE formulation

One can eliminate the unknown pressure term by taking the curl of (2.1). This signifies that the pressure is not an active but a passive variable. Moreover, the single-fluid, mass-averaged plasma velocity $\boldsymbol{v}$ , vorticity $\boldsymbol{\omega }\equiv \boldsymbol{\nabla }\,{\times}\, \boldsymbol{v}$ , magnetic $\boldsymbol{B}$ and current density $\boldsymbol{J}$ vector fields are divergence free, and they admit then the following representations:

(2.9) \begin{align} \boldsymbol{v} &= \frac {1}{R}\boldsymbol{\nabla }\chi \,{\times}\, \boldsymbol{i}_{\varphi }+v_{\varphi }\boldsymbol{i}_{\varphi } ,\\[-10pt]\nonumber \end{align}

(2.10) \begin{align} \boldsymbol{\omega }&=\frac {1}{R}\boldsymbol{\nabla }\left ( Rv_{\varphi }\right ) \,{\times}\, \boldsymbol{i}_{\varphi }-\frac {1}{R}\left ( \triangle ^{\ast }\chi \right ) \boldsymbol{i}_{\varphi },\\[-10pt]\nonumber \end{align}

(2.11) \begin{align} \boldsymbol{B} &= \frac {1}{R} \boldsymbol{\nabla }\psi \,{\times}\, \boldsymbol{i}_{\varphi } + B_{\varphi } \boldsymbol{i}_{\varphi },\\[-10pt]\nonumber \end{align}

(2.12) \begin{align} \boldsymbol{J} &=\frac {1}{R}\boldsymbol{\nabla }( RB_{\varphi }) \,{\times}\, \boldsymbol{i}_{\varphi }-\frac {1}{R}\left ( \triangle ^{\ast }\psi \right ) \boldsymbol{i}_{\varphi } , \end{align}

where $\chi$ is the velocity streamfunction, $\psi$ is the magnetic flux function, $B_{\varphi }$ is the toroidal component of the magnetic field vector and $v_{\varphi }$ is the toroidal component of the velocity field vector. The above system of (2.1)–(2.5) with the external drives (2.6)–(2.8) can be expressed Kamp et al. (Reference Kamp, Montgomery and Bates1998), Kamp & Montgomery (Reference Kamp and Montgomery2003), Kamp & Montgomery (Reference Kamp and Montgomery2004), Krupka & Firpo (Reference Krupka and Firpo2024); Roverc’h et al. (Reference Roverc’h, Oueslati and Firpo2021) as the following set of five scalar elliptic PDEs:

(2.13) \begin{align} \triangle ^{\ast }\chi &=-R\omega _{\varphi },\\[-10pt]\nonumber \end{align}

(2.14) \begin{align} \nonumber \nu \triangle ^{\ast }(R\omega _{\varphi })&=\frac {1}{R^{2}}\frac {\partial }{\partial Z}((RB_{\varphi })^{2}-(Rv_{\varphi })^{2})+\frac {1}{R}\{ Rj_{\varphi },\psi \}\\ &\quad +\frac {1}{R}\big\{ \chi ,R\omega _{\varphi }\big\}+ \frac {2\omega _{\varphi }}{R}\frac {\partial \chi }{\partial Z}-\frac {2j_{\varphi }}{R}\frac {\partial \psi }{\partial Z},\\[-10pt]\nonumber \end{align}

(2.15) \begin{align} \eta \triangle ^{\ast }(RB_{\varphi })&=\frac {1}{R}\big\{ \chi ,RB_{\varphi }\big\} +\frac {1}{R}\big\{ Rv_{\varphi },\psi \big\} +\frac {2RB_{\varphi }}{R^{2}}\frac {\partial \chi }{\partial Z}-\frac {2v_{\varphi }}{R}\frac {\partial \psi }{\partial Z},\\[-10pt]\nonumber \end{align}

(2.16) \begin{align} \nu \triangle ^{\ast }(Rv_{\varphi }) &=\frac {1}{R} \big\{ \chi ,Rv_{\varphi }\big\} +\frac {1}{R}\big\{ RB_{\varphi },\psi \big\},\\[-10pt]\nonumber \end{align}

(2.17) \begin{align} \triangle ^{\ast }\psi &=-Rj_{\varphi }, \end{align}

with the toroidal projection of Ohm’s law giving the constraint

(2.18) \begin{equation} \eta Rj_{\varphi }=E_{0}+\frac {1}{R}\left \{ \psi ,\chi \right \}. \end{equation}

Here, the Poisson bracket $\{u,v\}$ for any spatial functions $u$ and $v$ is defined as

(2.19) \begin{equation} \left \{u,v\right \} \equiv \frac {\partial u}{\partial R}\frac {\partial v}{\partial Z}-\frac {\partial v}{\partial R}\frac {\partial u}{\partial Z}, \end{equation}

and the operator $\triangle ^{\ast }$ is defined by

(2.20) \begin{equation} \triangle ^{\ast } A \equiv \frac {\partial ^{2}A}{\partial R^{2}} - \frac {1}{R} \frac {\partial A}{\partial R} +\frac {\partial ^{2}A}{\partial Z^{2}}. \end{equation}

A relevant dimensionless control parameter has been identified by Montgomery (Reference Montgomery1993), Cappello & Escande (Reference Cappello and Escande2000) and Krupka & Firpo (Reference Krupka and Firpo2024) as the Hartmann number, $H = (\eta \nu )^{-1/2}$ . In fusion-relevant conditions, this is expected to be a large parameter, ranging from $10^6$ to $10^8$ . Simulations for realistic parameter values already exist for this system under various boundary conditions, typically using the JET geometry. The elliptic system (2.13)–(2.17) requires five boundary conditions. The four conditions associated with the divergence-free properties of the magnetic field ( $\boldsymbol{B}$ ), current density ( $\boldsymbol{J}$ ), velocity ( $\boldsymbol{V}$ ) and vorticity ( $\boldsymbol{\omega }$ ) vector fields are determined by ensuring the continuity of their normal components across the plasma boundary. The following boundary conditions are selected in the numerical simulations: $\chi = \psi = 0$ and $B_{\varphi } =1/R$ on $\partial \varOmega$ . We enforce Neumann boundary conditions on both the toroidal velocity $v_{\varphi }$ and toroidal vorticity $\omega _{\varphi }$ through $\partial _n v_{\varphi } = \partial _n \omega _{\varphi } = 0$ on $\partial \varOmega$ . We used the open-source PDE solver FreeFem++, employing the finite element method of Hecht (Reference Hecht2012) to solve the above steady-state axisymmetric system of equations in a weak form on the plasma cross-section domain $\varOmega$ with the specified boundary conditions. For our calculations, we set a tolerance parameter $\epsilon = 10^{-10}$ , allowing the Newton–Raphson scheme to converge in typically 4–5 iterations.

2.3. Examination of the pressure field

Let us now examine the pressure profile in the visco-resistive model (2.1)–(2.5) with the external drives (2.6)–(2.8). Let us assume for now that the steady-state plasma speed is negligible. Then, in the ideal limit, $\eta \to 0$ and $\nu \to 0$ , the steady-state Navier–Stokes equation (2.1) takes the form

(2.21) \begin{equation} \mathbf{\boldsymbol{\nabla }} p =\boldsymbol{J}\,{\times}\,\boldsymbol{B}. \end{equation}

Restricting ourselves to axisymmetric solutions, the projection of this equation on $R$ and $Z$ gives, respectively,

(2.22) \begin{eqnarray} \frac {\partial p}{\partial R}=R^{-1} \left ( -B_{\varphi } \frac {\partial (RB_{\varphi })}{\partial R} + j_{\varphi } \frac {\partial \psi }{\partial R}\right ), \end{eqnarray}

(2.23) \begin{eqnarray} \frac {\partial p}{\partial Z}=R^{-1} \left ( -B_{\varphi } \frac {\partial (RB_{\varphi })}{\partial Z} + j_{\varphi } \frac {\partial \psi }{\partial Z}\right ). \end{eqnarray}

In the toroidal direction, we get

(2.24) \begin{equation} 0=R^{-2} \{\psi ,RB_{\varphi }\},\end{equation}

which amounts to the well-known property of Grad–Shafranov’s theory that the diamagnetic function, $R B_{\varphi }$ , is a function of the magnetic flux $\psi$ only. Moreover, writing that $\boldsymbol{J}\boldsymbol{\cdot} \boldsymbol{B}$ is curl free, which follows from the force balance (2.21), and projecting this on the toroidal direction yields

(2.25) \begin{equation} -2 RB_{\varphi }\frac {\partial B_{\varphi }}{\partial Z} + \{\psi ,Rj_{\varphi }\} +2 Rj_{\varphi } \frac {\partial \psi }{\partial Z}=0. \end{equation}

Then, combining (2.23) and (2.25), the pressure gradient along the $z$ -axis with a zero-flow hypothesis is given by

(2.26) \begin{equation} \frac {\partial p}{\partial Z}=\frac {1}{2 R} \left \{Rj_{\varphi },\psi \right \}. \end{equation}

Yet, assuming no plasma flow, the toroidal projection of Ohm’s law in (2.18) states that $Rj_{\varphi }$ is a constant, with $Rj_{\varphi } = E_{0}/\eta$ . Equation (2.26) indicates then that the pressure field does not depend on $Z$ . However, from the set of (2.22)–(2.24), we can deduce that the pressure is a function of the magnetic flux $\psi$ such that $\{p, \psi \} = 0$ . Thus, we have $\partial _{Z} p = p'(\psi ) \partial _{Z} \psi = 0$ . This implies that the pressure profile is constant. This aligns with the results obtained by Braams (Reference Braams1991), which indicate that in equilibrium configurations, the current density must be proportional to $1/R$ when the pressure gradient is zero.

2.4. Implications and implementation of non-inductive current drives

Section 2.3 demonstrates that the sole inclusion of Ohm’s law to close the system imposes significant limitations on the model. Specifically, the effective pressure in the model arises only because the toroidal geometry and viscous dissipation prevent the steady-state velocity field from being identically zero. This allows the pressure profile to remain non-zero, as $R j_{\varphi }$ is not exactly constant. However, the model lacks a robust mechanism to provide sufficient heating to achieve fusion conditions. Therefore, incorporating alternative heating methods is essential to reaching higher pressure. We will also see that this tends to induce higher plasma rotation velocities in specific drive configurations.

In our previous analysis, we focused on the behaviour of the system for a specific ratio of $E_{0}/\eta$ , which was the only explicit drive in the dimensionless system of equations (Krupka & Firpo (Reference Krupka and Firpo2024). Now, we aim to introduce an additional drive which will manifest as an extra term in the toroidal component of Ohm’s law in (2.5) with

(2.27) \begin{equation} \frac {E_{0}}{R}\mathbf{+}\left ( \boldsymbol{v}\,{\times}\, \boldsymbol{B}\right ) \boldsymbol{\cdot }\boldsymbol{i}_{\varphi } + j_D= \eta \boldsymbol{J}\boldsymbol{\cdot }\boldsymbol{i}_{\varphi }, \end{equation}

where $j_D$ represents a current drive. Equation (2.18) now becomes

(2.28) \begin{equation} \eta Rj_{\varphi }=E_{0}+\frac {1}{R}\left \{ \psi ,\chi \right \}+ j_D R. \end{equation}

Our goal is to investigate the influence of the non-inductive current drive $j_D$ within our system. We will begin by evaluating its effect on the pressure profile. For this purpose, it is necessary to determine the pressure field, which we will now establish.

4.1. Pressure field behaviour without and with non-inductive current drives

Let us now solve Poisson’s equation for pressure, assuming a zero pressure condition at the boundary $\partial \varOmega$ . It is important to note that, in the absence of an additional toroidal current drive and assuming no plasma flow, we previously inferred a zero pressure gradient, resulting in a constant zero pressure in the limits $\eta \to 0$ and $\nu \to 0$ , as discussed in § 2.3. To return to dimensional pressure and compare the results with those obtained from the JET tokamak, we recall that $p^{\ast }$ is the dimensionless total pressure Krupka & Firpo (Reference Krupka and Firpo2024), normalised as $p^{\ast } = \hat {p}^{\ast }/v_{A0}^2 \rho _m$ , where $\hat {p}^{\ast }$ is dimensional pressure. By using parameters from a specific JET deuterium–tritium shot on JET (Team 1992), we obtain $\rho _m = 2.09 \times 10^{-7} \, \text{kg m}^{-3}$ and the Alfvén velocity $v_{A0} = 5.46 \times 10^{6} \, \text{m s}^{-1}$ .

Figure 1. Pressure field in Pascal units computed without the application of the drive ( $j_D=0$ ) for a Hartmann number of $H=10^5$ .

To verify the pressure distribution in the absence of the drive with plasma flow, let us examine figure 1. The pressure profile is presented in Pascal units (Pa). The order of magnitude of the pressure field in the absence of the current drive turns out to be unrealistically small, as predicted in § 2.3.

Figure 2. Toroidal current field without the drive ( $j_D=0$ ) (top) and with the drive $j_D$ set to $A=100$ and offset $B=0$ (bottom), for a Hartmann number of $H=10$ in dimensionless units.

To explore the effects of a non-inductive current drive, we considered a family of drives, $j_D$ , which are solutions to Poisson’s equation $\mathbf{{\nabla} }^{2} j_D = -A$ , with the boundary condition $j_D = B$ on $\partial \varOmega$ . Here, $A$ denotes the magnitude of the drive, while $B$ represents the current distribution offset. This approach provides an initial method for simulating current distributions akin to those generated by heating mechanisms in tokamaks, effectively ‘adding a bump’ to the current profile. Let us also note here that, by introducing this form of the drive in (2.28), $ E_0/\eta$ start to plays the same role as an offset $B$ as it is a constant. Therefore, in the presence of a current drive, we can set $E_0 = 0$ without loss of generality, as this simply shifts $B$ to $B - E_0/\eta$ . This amounts to a purely non-inductive tokamak operation.

Let us now choose a drive that produces realistic pressure profiles. To do so, we select $j_D$ with $A = 100$ and $B = 0$ . Figure 2 compares the toroidal current density fields, calculated at Hartmann number $H = 10$ , using (2.13)–(2.17) for $j_D = 0$ (the reference case) and for $j_D$ with $A = 100$ and $B = 0$ . In the original system of Kamp & Montgomery (2003, 2004), Oueslati et al. (Reference Oueslati, Bonnet, Minesi, Firpo and Salhi2019) and Krupka & Firpo (Reference Krupka and Firpo2024) ( $j_D = 0$ ), the model fails to produce realistic toroidal current density profiles despite yielding a realistic total current. This is because a ratio of $E_{0}/\eta$ of approximately one corresponds to a realistic total current for the JET deuterium–tritium shot described in JET (Team 1992).

Figure 3. Pressure field in Pascal units computed with the application of the drive $j_D$ with $A=100$ and $B=0$ on the toroidal current density field. The Hartmann number is $H=10^5$ .

Let us now examine how the application of the drive affects the pressure profiles. Figure 3 shows the computed pressure, $p = p^{\ast } - v^2/2$ , with the application of the drive $j_D$ using $A = 100$ and $B = 0$ for $H = 10^5$ . The current drive not only prevents unrealistically low pressure levels but also establishes a realistic central pressure. In our incompressible framework, what matters physically is the pressure gradient or difference, since one can always add a constant to the pressure field without affecting the dynamics. In our case, the boundary pressure is fixed at zero, which makes the centre-to-edge pressure drop a meaningful quantity. For the presented configuration, the central pressure reaches values around $ 4 \,{\times}\, 10^6$ Pa, in line with experimental observations in the JET tokamak by JET (Team 1992). This is consistent with the observations in Chahine & Bos (Reference Chahine and Bos2018), where the role of absolute pressure values is discussed in the context of incompressible MHD. As pointed out there, pressure differences, and not the absolute value, are what influence the dynamics, especially when interpreting numerical data or comparing with realistic high- $\beta$ experiments. Our results, in this light, confirm that the application of the drive produces physically relevant pressure gradients and magnitudes that fall within realistic tokamak operating regimes.

Let us now examine how the variation of the magnitude of the drive affects the root-mean-square of the pressure as a function of the Hartmann number. Figure 4 illustrates this relationship. It is evident that all the magnitudes $A$ of the drive $j_D$ chosen in this study produce a realistic pressure response.

Figure 4. Root mean square of the pressure field in Pascals as a function of the Hartmann number, with the application of the drive $j_D$ with $B=0$ on the toroidal current field, for various values of $A$ .

Let us note that, with the application of the drives, it is possible to achieve various configurations of magnetic flux surfaces, including non-nested magnetic field lines with several $n=0$ islands present. Such an example is given in figure 5, which shows the magnetic flux surfaces and the pressure profile when the drive $j_D$ with $A=100$ and $B=-5$ is applied. However, we will primarily focus on drives that induce standard nested magnetic flux surfaces.

Figure 5. Magnetic flux surfaces with internal separatrices (on the left) and pressure profiles (on the right) computed with the application of the drive $j_D$ with $A=100$ and $B=-5$ on the toroidal current field for $H=10^5$ .

Figure 6. Magnetic flux surfaces and pressure isolines computed for the drive $j_D$ with $A = 100$ and $B = 0$ with $H = 1$ (non-ideal case) on the left and $H = 10^3$ (approaching the ideal limit) on the right.

Figure 7. (Left) Dimensionless pressure and diamagnetic function plotted as functions of the poloidal magnetic flux in the $H = 1000$ MHD simulation, along with their best cubic polynomial fits. The $(\psi , P)$ and $(\psi , F)$ plots are constructed by evaluating $P$ , $F$ and $\psi$ at the same finite element mesh nodes $(x_i, y_i)$ . Each pair $(\psi _i, P_i)$ and $(\psi _i, F_i)$ is represented as a red point. (Right) Comparison of the same levels of the normalised poloidal magnetic flux, $\psi _N$ , obtained from the $H = 1000$ visco-resistive MHD simulation (blue curves) and from the Grad–Shafranov equilibrium reconstruction (red curves) using the pressure and diamagnetic functions fitted in the left panel.

4.2. Recovery of Grad–Shafranov equilibrium in the ideal MHD limit

We now consider the case of standard nested magnetic flux surfaces, induced by a non-inductive current drive $j_D$ with magnitude $A = 100$ and offset $B = 0$ . Let us first examine the non-ideal regime, where the dimensionless viscosity $\nu$ and resistivity $\eta$ are both set to 1. This corresponds to a Hartmann number of $H = 1$ , and represents a strongly dissipative case with significant viscous and resistive effects. In this regime, shown on the left of figure 6, the pressure isolines do not align with the magnetic flux surfaces. As we approach the ideal MHD limit, i.e. as $\nu \to 0$ and $\eta \to 0$ (or equivalently, $H \to \infty$ ), the dynamics becomes dominated by the ideal force balance. In this case, shown on the right of figure 6 for $H = 10^3$ , the pressure isolines start to align with the magnetic flux surfaces. This result shows that the pressure tends to become a flux function in the ideal limit.

To explore this further, we extract the data for the pressure field $P$ and the diamagnetic function $F = x B_\varphi$ , and plot them as functions of the corresponding values of the poloidal magnetic flux, $\psi$ . In the left-hand plots of figure 7, we show that, for $H = 1000$ , both $P$ and $F$ are very well approximated by cubic polynomial functions of $\psi$ . They behave then as flux functions. Then we use these fits in a code that solves the Grad–Shafranov equation using a Picard fixed-point iteration.

The right-hand plot of figure 7 compares the level curves of the normalised poloidal magnetic flux, $\psi _N$ , obtained for the same $\psi _N$ values from a self-consistent visco-resistive MHD code at $H = 1000$ , with those obtained from the Grad–Shafranov equilibrium reconstruction. A very good agreement is observed, as expected, since the non-ideal and finite-velocity contributions are small. Quantitatively, using the $L^2$ -norm

(4.1) \begin{equation} \left \| \psi _1 - \psi _2 \right \|_{L^2(\varOmega )} = \left ( \int _\varOmega \left ( \psi _1(x, y) - \psi _2(x, y) \right )^2 \, \textrm{d}x\,\textrm{d}y \right )^{1/2}, \end{equation}

we get in this case $\left \| \psi _{N,GS} - \psi _{N,MHD} \right \|_{L^2(\varOmega )} \approx 0.015$ .

This indicates that the Grad–Shafranov equilibrium can be recovered from the self-consistent MHD solutions in the appropriate limits, but this also highlights the broader validity of our model in more dissipative regimes. Most importantly, the pressure field is self-consistently determined within our system, which is not the case for Grad–Shafranov equilibrium reconstructions.

4.3. Impact of the non-inductive current drive on steady-state velocity and scaling

Let us now examine the impact on the velocity distribution of the application of the current drive $j_D$ . Figure 8 presents the root mean square of the toroidal velocity field while applying the drive $j_D$ with $B = 0$ to the toroidal current field across various values of $A$ , as in figure 4. It can be observed that, while varying the magnitude of the drive causes an increase in velocities in the low- $H$ regime, the large- $H$ , boundary layer regime remains almost unchanged, despite the application of current drives with different magnitudes.

In Krupka & Firpo (Reference Krupka and Firpo2024), we showed that, at large Hartmann numbers, the velocity field develops a distinct boundary layer that becomes progressively thinner as $H$ increases. In this regime, the boundary layer thickness and, consequently, the velocity, scale with the Hartmann number, if nonlinear effects remain negligible. Specifically, by analysing the boundary layer equations analytically, we found that both the toroidal and poloidal velocities scale as $H^{1/4}$ for $H \gg 1$ . This scaling was also confirmed numerically through power-law fitting. Increasing the magnitude of the drive raises the total current, but the velocities appear unaffected by this variation. This insensitivity is consistent with the persistence of the boundary layer regime, which imposes a fixed scaling of the velocities with the Hartmann number, independent of the drive amplitude.

Figure 8. Root mean square of the toroidal velocity in Alfvén velocity units as a function of the Hartmann number, considering the application of the drive $j_D$ with $B=0$ on the toroidal current field, for the different values of $A$ .

Figure 9. Root mean square of the toroidal velocity in Alfvén velocity units as a function of the Hartmann number, considering the application of the drive $j_D$ with $A=100$ on the toroidal current field, for various values of $B$ .

Next, let us examine how the velocities change with the variation of the parameter $B$ , which represents the offset of the drive $j_D$ . Figure 9 presents the same information as figure 8, but with a fixed value of $A=100$ while exploring different values of $B$ . Shifting the drive results in the highest velocities at $B=-5$ . The usual large- $H$ velocity behaviour changes at certain parameters of the current drive. Here, we consider several different toroidal current profiles: two that lie entirely in the positive region (e.g. $B=3$ and $B=0$ ), two that are in between positive and negative regions ( $B=-3$ and $B=-5$ ) and one that is fully in the negative region $B=-7$ . In some of our simulations, the self-consistent MHD evolution leads to the formation of reversed current profiles, where the toroidal current density changes sign across the plasma radius. These current configurations give rise to additional magnetic axes and separatrices, modifying the magnetic topology beyond the typical single-axis case. Such phenomena are not only robust numerical features but are also physically relevant, as reversed current profiles and associated topological changes have been experimentally observed in devices like JET and JT-60U. Moreover, we found that these configurations coincide with a breakdown of standard scaling laws for the velocity fields, showing increased root-mean-square values for both toroidal and poloidal velocities. This makes them particularly interesting for understanding fundamental transport and equilibrium behaviour in such regimes.

Figure 10. Toroidal current density field (on the left) and toroidal velocity field (on the right) with the application of the drive $j_D$ with $A=100$ , $B=-5$ for $H=10^5$ as in figure 5.

Indeed, in the high-Hartmann-number regime, two scenarios emerge: either the boundary layer forms Krupka & Firpo (Reference Krupka and Firpo2024), and the root mean square of the toroidal and poloidal velocities exhibit some scaling law with the Hartmann number, as is the case for the drive $j_D$ with $B=0$ and $B=3$ , or the velocities do not develop a linear behaviour on a log–log scale for $B=-3$ , $-5$ and $-7$ , meaning that the velocities do not scale with the Hartmann number. This departure from scaling behaviour is probably linked to the way the toroidal velocity field changes with the drive offset $B$ . For different values of $B$ , the velocity distribution can differ significantly, which affects whether or not the boundary layer forms. In our view, what we observe here is a competition between different scaling regimes, which leads to the non-monotonic behaviour. Let us now take a closer look at this latest phenomenology, which is novel and pertains to a situation involving non-nested magnetic flux surfaces with internal separatrices. To illustrate this, we examine the non-inductive current drive $j_D$ with $A = 100$ and $B = -5$ for $H = 10^5$ , where the magnetic flux surfaces and pressure field are depicted in figure 5. Figure 10 shows on the left the associated toroidal current field.

The velocity distribution on the right of figure 10 closely resembles the toroidal current distribution, with the highest velocities occurring at the transition point between positive and negative current regions. In this case, we observe no formation of a boundary layer, which is advantageous from a numerical perspective. The absence of a boundary layer contributes to increased stability in the code and yields more robust results. The toroidal velocity field in our simulations exhibits a clear up–down antisymmetry, such that $ v_\varphi (x, y) = -v_\varphi (x, -y)$ , leading to zero net toroidal momentum. This symmetry arises naturally from the geometry and boundary conditions imposed in our set-up. While this excludes the possibility of spontaneous momentum generation in the present context, it is worth noting that previous studies have demonstrated mechanisms through which up–down symmetry breaking can lead to the emergence of net toroidal flow, both in gyrokinetic frameworks (Camenen et al. Reference Camenen, Peeters, Angioni, Casson, Hornsby, Snodin and Strintzi2009) and in MHD (Morales et al. Reference Morales, Bos, Schneider and Montgomery2012). For instance, it was shown in Oueslati & Firpo (Reference Oueslati and Firpo2020) that breaking this symmetry by means of external magnetic perturbations can influence flow topology and might be used to induce finite mean flows. There is significant potential to achieve much higher velocities with these drives; however, accurately predicting which current drive would be optimal for maximising velocities remains a challenge and necessitates further investigation.