2. Visco-resistive magnetohydrodynamic description with heat convection-diffusion

We introduce now the modelling frame used to compute axisymmetric visco-resistive MHD steady states allowing for non-vanishing velocity fields. Rather than assuming a constant resistivity as in some previous studies, we shall allow here its space variation by imposing some temperature inhomogeneity. We shall assume that some extra heating is applied on the top of the tokamak so that there is a vertical gradient of the plasma boundary temperature. The temperature field within the tokamak plasma is assumed to be governed by a convection-diffusion (heat) equation (2.3g). We shall use the fact that the electrical resistivity $\eta$ varies with the temperature $T$ according to Spitzer's law

(2.1)\begin{equation} \eta (T)\equiv \alpha T^{{-}3/2} \end{equation}

with

(2.2)\begin{equation} \alpha \equiv \frac{4}{3}\frac{\sqrt{2{\rm \pi} m_{e}}e^{2}\ln \varLambda }{\left( 4{\rm \pi} \varepsilon _{0}\right) ^{2}}, \end{equation}

where $\ln \varLambda$ denotes the Coulomb logarithm.

The system of equations considered in the present study reads then

(2.3a)\begin{gather} (\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} )\boldsymbol{v} = \boldsymbol{J}\times \boldsymbol{B}-\boldsymbol{\nabla} p+\nu \nabla^{2}\boldsymbol{v}, \end{gather}

(2.3b)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{v} = 0, \end{gather}

(2.3c)\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot}\boldsymbol{B} = 0, \end{gather}

(2.3d)\begin{gather}\boldsymbol{\nabla} \times \boldsymbol{E} = \boldsymbol{0}, \end{gather}

(2.3e)\begin{gather}\boldsymbol{\nabla} \times \boldsymbol{B} = \boldsymbol{J}, \end{gather}

(2.3f)\begin{gather}\boldsymbol{E}+\boldsymbol{v}\times \boldsymbol{B} = \eta (T) \boldsymbol{J}, \end{gather}

(2.3g)\begin{gather}(\boldsymbol{v}\boldsymbol{\cdot}\boldsymbol{\nabla} )T = \chi_{{T}} \triangle T. \end{gather}

The physical unknown quantities are the electric field $\boldsymbol {E}$, the magnetic field $\boldsymbol {B}$, the velocity field $\boldsymbol {v}$ and the scalar pressure $p$ and temperature $T$ fields. The conducting plasma fluid is assumed to have a uniform density so that the mass conservation equation turns into the incompressibility condition (2.3b). The steady-state Navier–Stokes equation normalized to the mass density, $\rho$, is written in (2.3a) where the kinematic viscosity $\nu$ has been assumed to be scalar and constant. These fluid equations are coupled to the Maxwell equations (2.3c)–(2.3e), to Ohm's law (2.3f) and to the heat steady-state convection-diffusion equation (2.3g). This system of equations is written for dimensionless variables in the usual Alfvèn units where velocities are normalized by the characteristic Alfvèn velocity $v_{{A}}=(B^{2}/\mu _{0}\rho )^{1/2}$, where $\mu _{0}$ is the permeability of free space. Concerning (2.3g), let us mention that it does not include the heat generation by viscous dissipation. Since we have in mind plasmas for magnetic confinement fusion applications that have very low viscosity, the temperature inhomogeneity induced by this term is expected to be negligible in front of the one induced by an external inhomogeneous heating. Besides, since our focus is on toroidally invariant steady states, we retain in the anisotropic heat flux only the transverse part and identify this to $- \chi _{{T}} \boldsymbol {\nabla } T$. The diffusion coefficient $\chi _{{T}}$ is typically of the order of $1\ \mathrm {m}^2\,\mathrm {s}^{-1}$ in tokamaks (Freidberg Reference Freidberg2007) and this numerical value is taken in our simulations.

The system (2.3) reduces to the system of coupled parabolic differential equations (A15) presented in Appendix A. It needs to be completed by the specification of the problem geometry and boundary conditions. As for the problem geometry, we consider the JET parameters with major radius $r_{0}=3\ \mathrm {m}$, semi-minor axis $r_{1}=1.25\ \mathrm {m}$, semi-major axis $r_{2}=\kappa r_{1}$ with plasma elongation $\kappa = 1.55$ and triangularity, δ, given by $\arcsin \delta = 0.5$. Then, six boundary conditions have to be given. Four of them are associated with the divergence-free nature of the magnetic $\boldsymbol {B}$, current-density $\boldsymbol {J}$, velocity $\boldsymbol {v}$ and vorticity $\boldsymbol {\omega }$ vector fields and may be fixed by the continuity of their normal component on the plasma boundary. There remains then two conditions on the toroidal vorticity and temperature. As detailed below, one chooses some vertically inhomogeneous heating that fixes the boundary condition for $T$; whereas one assumes, with some arbitrariness, that the toroidal vorticity vanishes on the plasma border. This is written explicitly in Appendix A.

There are two (small) dimensionless parameters in the problem, that are the resistivity $\eta$, which is the inverse of the magnetic Reynolds number $S$, and the kinematic viscosity $\nu$, which is the inverse of the viscous Lundquist number $M$. Here the modelling of the resistive and viscous effects has been made in the (usual) tractable scalar way. There are also two driving parameters in the model equations associated with the magnitude of the external curl-free toroidal magnetic field, $B_0$, and to the toroidal loop voltage, $E_0$. The equations were solved under their weak form using the finite element method with FreeFem++ (Hecht Reference Hecht2012; Oueslati et al. Reference Oueslati, Bonnet, Minesi, Firpo and Salhi2019).

3. Simulation results with localized heating

Compared with previous frameworks, temperature is considered here as a variable of the MHD visco-resistive system of equations. In order to study the impact of a localized heating in breaking the symmetry of the problem, we shall mainly use a vertically linear profile for the temperature on the boundary. This external temperature gradient can be characterized by the relative difference $\Delta T$ between the top temperature, $T_{\mathrm {up}}$, and bottom temperature, $T_{0}$, of the plasma through $\Delta T = (T_{\mathrm {up}}-T_{0})/T_{0}$. Physically speaking, we consider positive $\Delta T$ corresponding to an inhomogeneous extra heating directed on the top of the tokamak plasma. On the numerical side, the introduction of a non-vanishing external temperature gradient interestingly happens to stabilize the code compared with its fixed-temperature original version (Oueslati et al. Reference Oueslati, Bonnet, Minesi, Firpo and Salhi2019) and allows larger Hartmann numbers to be reached.

Let us define by $\eta _{0}=\eta (T_{0})$ the value of the resistivity at the bottom of the tokamak plasma. In the numerical simulations, we used the dimensionless parameters $B_{0}=0.87$ and $E_{0}=3 \times 10^{-9}$ relevant to JET and a value of the maximal, bottom, resistivity $\eta (T_{0})=6.9 \times 10^{-7}$. It is well known that the order of magnitude of the effective (scalar) plasma viscosity in a tokamak is uncertain (Kamp & Montgomery Reference Kamp and Montgomery2004), therefore, we run simulations for a large span of values of $\nu$, or, identically, a large span of Hartmann numbers, where the Hartmann number $H$ is defined as a combination of $\eta _{0}$ and $\nu$ through $H=1/\sqrt {\eta _{0}\nu }$. Alternatively, we shall consider the mean Hartmann number $\langle H \rangle$ defined as the average value of $H$ in the tokamak domain.

Let us first investigate the fate of the temperature field. One may actually think that the introduction of the stationary heat equation in the model would not have any significant effect on the value of $T$ in the bulk, compared with simply maintaining an affine temperature everywhere. Indeed, only convective transport of temperature can make $T$ differ from the Laplace equation solutions. Because poloidal speed was typically very low in previous simulations (Oueslati et al. Reference Oueslati, Bonnet, Minesi, Firpo and Salhi2019), this would seem to suggest a rather limited effect. Figure 1 partly confirms this intuition. For $H = 10$, convection is clearly not visible globally. Yet, for $H = 10^5$, the temperature profile happens to be deviated by convective effects, meaning that, in this case, poloidal flows are sufficient to provide some heat transport. One observes there some tendency for a temperature homogenization in the tokamak core at large $H$.

Figure 1. Temperature profiles (in eV units) for $\Delta T = 10$ and Hartmann numbers (a) $H = 10$, (b) $H = 10^5$.

Our main interest is to examine the impact of the inhomogeneity of the heating on the steady-state plasma toroidal speed. In the present study, this generation of the toroidal flow is essentially due to viscous effects, as emphasized in the seminal study (Montgomery et al. Reference Montgomery, Bates and Li1997), meaning that, in the Navier–Stokes equation, the toroidal contribution of the $\boldsymbol {\omega } \times \boldsymbol {v}$ term increases with $H$ while remaining negligible in front of the toroidal contribution of the $\boldsymbol {J} \times \boldsymbol {B}$ term, that is almost independent on $H$.

Figure 2 represents the root mean square (r.m.s.) value of the toroidal velocity field as a function of the mean Hartmann number $\langle H \rangle$ for $\Delta T = 1$ in Alfvèn units. The numbers $M$ label different meshes. Up to $\langle H \rangle \times 10^{6}$, only one curve can be seen, which means that the simulations yield identical results no matter what mesh is used, illustrating the numerical robustness of the code (see also Appendix B for a discussion on numerical issues at large Hartmann numbers). For reference, the curve corresponding to a uniform temperature on the border ($\Delta T =0$) is also plotted. At any given mean Hartmann number, one observes that the r.m.s. toroidal speed is larger in the inhomogeneous heating case than for constant temperature with a difference increasing with $\langle H \rangle$. Figures 3 and 4 show velocity profiles that depart from previous constant temperature results (Oueslati et al. Reference Oueslati, Bonnet, Minesi, Firpo and Salhi2019) as the Hartmann number increases. One observes that the toroidal velocity field goes from being almost odd with respect to the midplane at $H=10$ to being almost even at $H=10^{5}$ with a net toroidal flow. The toroidal velocity in itself, although higher than in previous cases, is not particularly impressive at these Hartmann numbers. It is rather the rate at which the velocity increases with $\langle H \rangle$, as shown on figure 2, that catches our interest.

Figure 2. The r.m.s. value of $v_{\varphi }$ (in Alfvèn velocity $v_{{A}}$ units) depending on the mean Hartmann number $\langle H \rangle$ for $\Delta T = 0$ and $\Delta T = 1$ using different mesh sizes (labelled by the number of vertices on the edge of the domain).

Figure 3. Toroidal velocity field $v_{\varphi }$ in Alfvèn velocity $v_{{A}}$ units with $\Delta T =1$ for (a) $H = 10$ and (b) $H = 100$.

Figure 4. Toroidal velocity field $v_{\varphi }$ in Alfvèn velocity $v_{{A}}$ units with $\Delta T =1$ for (a) $H = 10^3$ and (b) $H = 10^5$.

We now attempt to map the results of the simulations for different $\Delta T$, even if the largest $\Delta T$ values may be physically unrealistic. One of our objectives is to obtain faster flows than were previously possible for a uniform temperature. The mean quadratic (r.m.s.) toroidal velocities for different edge temperature gradients are presented in figure 5 for a toroidal free-slip boundary condition and figure 6 for a toroidal no-slip boundary condition. The toroidal free-slip boundary condition allows somehow higher speeds to be reached than the no-slip one for a given $\langle H \rangle$ which is visible from a comparison between figures 5 and 6. Otherwise, the general behaviour of the curves in both plots is quite similar. The toroidal free-slip boundary condition has been taken in our simulations unless otherwise stated. As expected, for a given $\Delta T$, the toroidal velocity always increases with $\langle H \rangle$, as the viscous dissipation diminishes. Moreover, the mean value of $\eta$ has an easily identifiable effect on velocity: for low Hartmann numbers, the more intense the heating, the higher the toroidal velocity. Yet, for higher $\langle H \rangle$, the flow behaves differently. The plots corresponding to the most intense heatings are rapidly saturating: the slope of $\langle v_{\varphi } \rangle _{\mathrm {rms}}$ becomes small compared with simulations with less heating and the greatest temperature gradients do not yield the fastest flows. Conversely, curves with softer heating follow a linear regime on figure 6 over a large band of Hartmann numbers, with a slope equal to 2 (in the log-log scale). This corresponds to $\langle v_{\varphi } \rangle _{\mathrm {rms}}$ being proportional to $\langle H \rangle ^2$. This slope is steeper than what could be achieved without heating, and it provides a physically interesting regime. However, the higher the temperature gradient, the sooner this regime ends (in terms of $\langle H \rangle$). In order to quantify these effects, the instantaneous slope $\beta (\langle H \rangle , \Delta T) \equiv \partial (\log \langle v_{\phi } \rangle _{\mathrm {rms}})/ \partial (\log \langle H \rangle )$ has been plotted on figure 7. In the infinite viscosity limit, corresponding to $\langle H \rangle \rightarrow 0$, one recovers that the quadratic mean (or r.m.s.) of the toroidal velocities scales as $\langle H \rangle ^{4}$ (i.e. $\beta \rightarrow 4$) which was analytically predicted for $\Delta T=0$ in Kamp, Montgomery & Bates (Reference Kamp, Montgomery and Bates1998). One observes the aforementioned intermediate scaling $\langle H \rangle ^{2}$ with $\beta$ plateauing about 2 for moderate heating gradients compared with an intermediate scaling with $\beta$ approximately 3/2 taking place for $\Delta T \rightarrow 0$.

Figure 5. The r.m.s. of the toroidal velocity, $\langle v_{\varphi } \rangle _{\mathrm {rms}}$, (in Alfvèn velocity $v_{{A}}$ units) as a function of $\langle H \rangle$ for different values of $\Delta T$.

Figure 6. The r.m.s. of the toroidal velocity, $\langle v_{\varphi } \rangle _{\mathrm {rms}}$, (in Alfvèn velocity $v_{{A}}$ units) as a function of $\langle H \rangle$ for different values of $\Delta T$. Here the toroidal no-slip condition ($u_{4}=0$ on the border of the plasma domain) is used.

Figure 7. Instantaneous slopes of the logarithm of the r.m.s. of the toroidal velocity, $\langle v_{\varphi } \rangle _{\mathrm {rms}}$, as a function of the logarithm of $\langle H \rangle$ for different values of $\Delta T$ for the toroidal no-slip results plotted in figure 6.

All this makes it not easy to point out an optimal set of parameters in order to reach fast-flow states, especially since the real, physical value of the viscosity coefficient – and therefore the Hartmann number – remains unknown. Besides, simulations with high $\Delta T$ present little physical interest besides drawing a general picture of the situation, yet moderate $\Delta T$-values provide some interesting results: as seen on figures 2, 5 and 6, several orders of magnitude in velocity can be gained over the uniform-temperature simulations. On the basis of figure 1, one may propose an explanation for the dual role played by the temperature gradient $\Delta T$. Indeed, when both $\Delta T$ and $\langle H \rangle$ are large, the effective temperature field will tend to homogenize in the core of the tokamak. The effect of this homogenization is to make less effective or to counterbalance the symmetry breaking induced primarily by the inhomogeneous heating imposed by boundary conditions. One may then understand that there exists some optimum in the choice of $\Delta T$ in terms of the maximization of mean toroidal plasma speed and that a ‘moderate’ up-down temperature gradient $\Delta T$ may eventually be more efficient than a large $\Delta T$ to make plasma rotate in the large $\langle H \rangle$ values relevant to fusion tokamak plasmas.

Finally, it may be interesting to capture the effect of the additional, vertically ($y$)-inhomogeneous, heating in terms of the toroidal plasma rotation. Figure 8 represents the ratios $\langle v_{\varphi } \rangle _{\mathrm {rms}}(\Delta T \neq 0)/\langle v_{\varphi } \rangle _{\mathrm {rms}}(\Delta T=0)$ as a function of the Hartmann number at the bottom of the tokamak, $H$, for several values of $\Delta T$. This indicates that, for large fusion-relevant values of $H$ (H between $10^6$ and $10^8$), the toroidal plasma rotation is maximized by an optimal, moderate, differential heating corresponding to a $\Delta T$ of order one. At those large fusion-relevant values of $H$, numerical simulations become challenging due to the emergence of a boundary layer having a thickness decreasing with $H$. This is visible on the plots of the toroidal velocity profile of figure 9. A point on the numerical results at large $H$ is presented in Appendix B.

Figure 8. For the same data as in figure 6, plots of the ratio $\langle v_{\varphi } \rangle _{\mathrm {rms}}(\Delta T \neq 0)/\langle v_{\varphi } \rangle _{\mathrm {rms}}(\Delta T=0)$ as a function of $H$ (Hartmann number at the bottom of the tokamak) for different values of $\Delta T$.

Figure 9. Profiles of the toroidal velocity field (in Alfvèn velocity $v_{{A}}$ units) in the case of the toroidal no-slip boundary condition for $\Delta T=0.1$ and large-$H$ values equal to (from a–c) $2\times 10^{6}$, $5\times 10^{6}$ and $8\times 10^{6}$. Colours are varying with the relative magnitude of $v_{\varphi }$ for better visualization.

4. Spatial symmetry of solutions

4.1. Preliminary remarks

In order to characterize more precisely the flow modes described earlier, we shall now study the symmetry of solutions with respect to the tokamak horizontal $y=0$ midplane. By introducing an inhomogeneous heating along the vertical $y$ axis, the up-down symmetry of the problem has been broken. For the uniform-temperature problem, symmetry properties were studied in Oueslati & Firpo (Reference Oueslati and Firpo2020) and the toroidal velocity field was shown to be naturally antisymmetric, i.e. an odd function of $y$.

In the present simulations, for high Hartmann numbers we obtain some rather symmetric profiles (see e.g. Figure 4) in stark contrast to the non-heated case. This happens to be a somewhat paradoxical situation since breaking the symmetry of the problem seems to make the solutions go from antisymmetry to symmetry, rather than merely eliminating any symmetry properties. Nevertheless, the picture is more intricate as for instance, in the case of intense heating, there are no longer any obvious symmetry properties for high enough Hartmann numbers (see figure 10).

Figure 10. Toroidal velocity field $v_{\varphi }$ for $\Delta T = 50$ and $H = 100$.

4.2. Symmetry analysis

To go beyond the aforementioned qualitative observations and obtain more precise information, one must quantitatively characterize the symmetry or antisymmetry of our solutions. To that end, we consider the $L^2$ norm of the symmetric and antisymmetric parts of the toroidal velocity. Defining by $\varOmega$ the tokamak cross-section and by $S_{\varOmega }\equiv \int _{\varOmega }\,\mathrm {d}x\, \mathrm {d}y$ its surface, we put

(4.1)\begin{equation} V_{\mathrm{sym}} = \left( \frac{1}{S_{\varOmega}} \int \left( \frac{v_{\varphi}(x,y) + v_{\varphi}(x,-y)}{2} \right)^2\,\mathrm{d}x\,\mathrm{d}y \right) ^{1/2} \end{equation}

and

(4.2)\begin{equation} V_{\mathrm{anti}} = \left( \frac{1}{S_{\varOmega} } \int \left( \frac{v_{\varphi}(x,y) - v_{\varphi}(x,-y)}{2} \right)^2\,\mathrm{d}x\,\mathrm{d}y\right) ^{1/2}. \end{equation}

One has the identity

(4.3)\begin{equation} \langle v_{\varphi}^2 \rangle = V_{\mathrm{sym}}^2 + V_{\mathrm{anti}}^2 \end{equation}

which leads us to introduce a symmetry argument, $\chi$, given by $\chi \equiv \arctan ( V_{\mathrm {sym}}/ V_{\mathrm {anti}} )$. This scalar quantity allows us to easily visualize the balance between symmetry and antisymmetry in a given velocity field. This argument does not depend on the norm of the velocity itself (which can vary by many orders of magnitude) and it is valued between $0$ and ${\rm \pi} /2$. For simulations with uniform temperature, we get $\langle v_{\varphi }\rangle _{\mathrm {rms}} = V_{\mathrm {anti}}$ and $V_{\mathrm {sym}}=0$, and therefore the symmetry argument is $\chi =0$.

4.3. Link between velocity and symmetry

Figures 11 and 12 show a transition from antisymmetry to symmetry in the case $\Delta T = 0.1$ as the Hartmann number increases. For low Hartmann numbers, the toroidal velocity profile is indistinguishable from that of the simulations done with uniform temperature: it appears to be antisymmetric and the two plots of $\langle v_{\varphi } \rangle _{\mathrm {rms}}$ are identical. We then have a rather abrupt transition around $H = 1000$, and the symmetry argument tends towards ${\rm \pi} /2$ (total up-down symmetry) for large $H$. What catches our interest is that the transition from antisymmetry to symmetry occurs simultaneously with the separation of the r.m.s. velocities of the homogeneous and non-homogeneous heating cases; it appears that the flow goes from a uniform-temperature-like state to a faster, more symmetric state.

Figure 11. Toroidal velocity field $v_{\varphi }$ for $\Delta T = 0.1$ and (a) $H=10$, (b) $H = 10^5$.

Figure 12. (a) Symmetry argument $\chi =\arctan ( V_{\mathrm {sym}}/V_{\mathrm {anti}} )$ for $\Delta T=0.1$; (b) $\langle v_{\varphi }\rangle _{\mathrm {rms}}$ for $\Delta T$ equal to 0 and 0.1 as a function of $\langle H \rangle$.

This phenomenon is also encountered for other values of $\Delta T$, although the transition occurs at different $H$. The plots do not reach sufficiently high Hartmann numbers to feature the asymptotic symmetry, due to limitations in computing capabilities (somewhat paradoxically, for small values of $\Delta T$, the less intense the inhomogeneity of the heating, the more numerically demanding the code becomes) as visible on figure 13.

Figure 13. (a) Toroidal velocity $v_{\varphi }$ field for $H = 10^4$ and $\Delta T = 0.01$; (b) symmetry argument $\chi =\arctan ( V_{\mathrm {sym}}/V_{\mathrm {anti}} )$ for $\Delta T=0.01$; (c) $\langle v_{\varphi }\rangle _{\mathrm {rms}}$ for $\Delta T$ equal to 0 and 0.01 as a function of $\langle H \rangle$.

For higher values of $\Delta T$, however, other phenomena occur at high $H$ as depicted in figure 14. After transitioning towards symmetry, this new behaviour (that cannot be easily interpreted in terms of symmetry) seems to be correlated with the velocity saturation mentioned earlier.

Figure 14. (a) Symmetry argument $\chi =\arctan ( V_{\mathrm {sym}}/V_{\mathrm {anti}} )$ for $\Delta T=1$; (b) $\langle v_{\varphi }\rangle _{\mathrm {rms}}$ for $\Delta T$ equal to 0 and 1 as a function of $\langle H \rangle$.