Keywords: magnétohydrodynamique, tokamak, reversed field pinch, méthode de pénalisation en volume, méthode Fourier pseudo-spectrale, moment angulaire, influence de la toroïdicité magnetohydrodynamics, tokamak, reversed field pinch, volume penalization method, spectral method, toroidal momentum, quasi single helicity, toroidicity effects . . . . . . . . . . . . . . . dimensional instationary diffusion problem . . . . .

R ÉSUM É

La description magnétohydrodynamique est utilisée pour étudier les plasmas de fusion par confinement magnétique dans deux configurations: tokamak et reversed field pinch. Une méthode de Fourier pseudo-spectrale et une technique de pénalisation en volume sont employées pour résoudre les équations. La méthode de pénalisation permet d'introduire des conditions aux limites de Dirichlet et donc de faire varier facilement la géométrie considérée. Les simulations dans des géométries toroïdales de type tokamak montrent l'apparition spontanée de vitesses. Une importante composante toroïdale se développe si le système est peu dissipatif. Il est aussi montré que la brisure de symétrie dans la forme de la section du tore fait apparaître un moment angulaire toroïdal. Pour le Reversed Field Pinch on montre l'émergence de structures hélicoïdales. La forme de ces structures varie en fonction des coefficients de transport ainsi que du paramètre de pincement du champ magnétique imposé. Pour compléter l'étude on compare les résultats du tore aux calculs dans un cylindre périodique. Les différences dans la dynamique des deux cas sont mises en avant. Finalement les simulations sont confrontées à des expériences et un meilleur accord est observé entre simulation et expérience pour la géométrie toroïdale que pour la géométrie cylindrique.

Nuclear fusion is a possible source of energy. Indeed, it is the reaction that powers the sun. The process involves lighter elements than those used in fission reactions. Where in the latter uranium atoms are split, in fusion mainly hydrogen and its isotopes deuterium and tritium are combined to produce other particles like helium. In this process a large amount of energy becomes available. For the last 60 years actively research has been conducted to achieve controlled thermonuclear fusion on earth. Succeeding, this could bring to the planet a source of energy with almost unlimited resources, while generating only a very small amount of radioactive waste. During the past decades various devices and methods have been developed. Two different ways to obtain a controlled reaction are inertial fusion and magnetic fusion. In the former, high energy laser beams hit a target containing the fuel. The resulting shock waves compress and heat the deuterium-tritium mixture making the atoms fuse. In a magnetic fusion configuration to start the reactions very high temperatures have to be reached, of the order of 150 • 10 6 K [START_REF] Freidberg | Plasma Phys. and Fusion Energy[END_REF], hotter than the center of the sun. For practical reasons the pressure in a magnetic fusion reactor is far less than during a fusion reaction in the stars, and this lack of high pressure must be compensated by an important thermal agitation that increases the probabilities of having collisions fusing the particles.

The main idea behind magnetic fusion is that the very high temperature ionizes the media, forming a plasma that can be controlled, in principle, by means of a magnetic field. This is essential, because no material on earth can support the temperatures reached in a fusion reactor. The conducting fluid, or plasma, can therefore be isolated from the walls by a strong magnetic field. Different magnetic configurations have been studied for the past decades to achieve a controlled magnetic fusion reaction. Two of them, on which we will focus here, are the tokamak and the RFP (Reversed Field Pinch). These two devices have many points in common. The magnetic field is produced by external coils and is combined with the magnetic field induced by an electric current flowing through the plasma itself. Both configurations have a toroidal shape and the imposed magnetic field that confines the plasma is helical and surrounds the toroidally shaped chamber (see Fig. 1.1). The tokamak is technologically the most advanced configuration and has a high probability to achieve a controlled nuclear reaction. In this device a strong magnetic field confines the plasma. The toroidal component of this field is usually larger than the poloidal component. The reason for this is that, in this geometry it is tried to avoid magnetohydrodynamic instabilities that develop if the ratio of toroidal over poloidal field is lower than unity. In contrast, for the RFP the poloidal and toroidal components of the magnetic field are usually of the same order of magnitude. The plasma is then unstable, and several instabilities can develop making the conducting fluid more difficult to confine.

One of the key issues to make a fusion reactor work is the ability to increase the confinement time of the plasma. In practice this is equivalent to reducing the thermal losses, which are mainly caused by turbulent fluctuations. These fluctuations degrade the quality of the confinement and thereby reduce the performance of the fusion reactor. In the tokamak configuration it was discovered three decades ago (in the ASDEX experiment) [START_REF] Wagner | Regime of improved confinement and high beta in neutral-beam-heated divertor discharges of the ASDEX tokamak[END_REF] that, under certain circumstances, the turbulent activity is reduced, leading to a better confinement. This improved regime is called high confinement mode (H mode) in contrast with the commonly found low confinement state (L mode). In H mode the confinement time is roughly two times bigger than in L mode. This increment might seem small but it can be of great importance for the performance of future reactors like ITER (the international experimental reactor being constructed in Cadarache, France).

In existing fusion reactors, large bulk velocities are observed [START_REF] Rice | Central impurity toroidal rotation in ICRF heated Alcator C-Mod plasmas[END_REF][START_REF] Rice | Observations of impurity toroidal rotation suppression with ITB formation in ICRF and ohmic H mode Alcator C-Mod plasmas[END_REF]. In particular the H mode is characterized by toroidal velocities that are several times superior to the ones measured in L-mode [START_REF] Rice | Inter-machine comparison of intrinsic toroidal rotation in tokamaks[END_REF]. It is not yet well understood if these velocities are either at the origin, or a consequence of the low-to-high confinement transition (L-H transition). Different theories have been proposed to explain the spontaneous generation of toroidal velocity, mostly based on the turbulent transport of toroidal momentum generated at the tokamak edge (e.g., in Refs. [START_REF] Mattor | Momentum and thermal transport in neutralbeam-heated tokamaks[END_REF][START_REF] Peeters | Overview of toroidal momentum transport[END_REF][START_REF] Diamond | Physics of non-diffusive turbulent transport of momentum and the origins of spontaneous rotation in tokamaks[END_REF]), but no complete understanding of the phenomenon is available at present. Toroidal velocities, and the understanding of plasma velocities in general will be one of the main subjects of the present thesis.

The RFP possesses a great advantage compared to the tokamak which is the necessity of a lower magnetic field to confine the plasma, approximatively ten times smaller. The toroidal field has a weak reversal in the plasma edge, which is the reason for the name "Reversed Field Pinch". The weak RFP magnetic field yields several advantages: the use of normal magnets instead of superconductors, high mass density and the possibility of reaching thermonuclear temperatures without additional heating besides the Ohmic one [START_REF] Martin | Overview of RFX-mod results[END_REF]. The drawback of the device is the broad spectrum of magnetohydrodynamic resistive instabilities that increases the turbulence and thereby lowers the energy confinement performance. However, recently quasi-single helicity (QSH) states were observed in experiments. In these states the full turbulent regime is avoided [START_REF] Escande | Quasi-single-helicity reversed-field-pinch plasmas[END_REF]. There is a decrease of magnetic chaos and the formation of a coherent helical structure within the plasma. Further, as for the tokamak, in the RFP important velocities have been measured [START_REF] Prager | Dynamo and anomalous transport in the reversed field pinch[END_REF] and their origin also remains unclear.

In this thesis we will focus on the numerical study of these two magnetic fusion configurations (tokamak and RFP) using one of the simplest models describing a plasma, the magnetohydrodynamic (MHD) approximation. Far more sophisticated descriptions of fusion plasmas are at hand nowadays, such as kinetic or gyrokinetic approximations and at a detailed level these descriptions must be superior. However, it seems that at present, if one is interested in global dynamics of the full fusion plasma over timescales relevant to the bulk dynamics, MHD is the only possible description which allows to be treated numerically. With respect to previous fusion related studies our work differs in that we consider both RFP and tokamak dynamics in a fully toroidal domain taking into account nonzero viscosity and resistivity without assuming any macroscopic equilibrium. Several assumptions, such as uniform transport coefficients, isothermal, incompressible plasma dynamics, order unity magnetic Prandtl number, are obviously oversimplifications compared to reality, but they should be considered as a logical first step in understanding the basic viscoresistive MHD behavior, in toroidal geometry, while retaining the freedom to easily change the precise shape of the domain.

In order to allow an efficient computation of the full viscoresistive MHD equations, we use a Fourier pseudo-spectral method. To take into account boundary conditions we have implemented and validated a volume penalisation technique originally developed for hydrodynamics [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF]. The large advantage of this method is the flexibility on the considered geometry and the ease of implementation. This flexibility facilitates the comparison of different geometries and the assessment of the boundary effects on the dynamics of the considered physical systems. The strength of the combination of a Fourier pseudospectral solver with the penalization method is that the numerical code uses fast Fourier transforms, that can be easy parallelized, and at the same time a large variety of relatively complex geometries can be taken into account.

In Chap. II we present the MHD equations and detail the numerical algorithm. The penalization method is described and several validation test-cases are exposed. Different implementations of the penalization technique are compared and the order of convergence of the numerical method is presented. We consider two-and three-dimensional test-cases for hydrodynamic and magnetohydrodynamic systems † .

In Chap. III the study of the magnetohydrodynamic flows generated in toroidal geometries is presented. The magnetic configuration (as measured by the safety factor) is here close to that of the tokamak device. Two different toroidal geometries, one with a symmetric and the other with an asymmetric cross section are compared. The diffusivity of the system is varied and a transition in the flow topology is observed. The differences resulting of the change in the geometry are highlighted as well as the influence of the change of the ratio of the toroidal over the poloidal imposed magnetic field. In particular we present a mechanism which explains the spontaneous generation of toroidal angular momentum, which is also observed in existing tokamaks. The mechanism we observe seems to be generic since it is present even in the coarsest description of a fusion plasma: viscoresistive magnetohydrodynamics § .

In Chap. IV the toroidal dynamics are presented for a magnetic configuration close to a RFP reactor. The transition from a calm axisymmetric state to a fluctuating helical state are presented for a toroidal and for a cylindrical geometry. The differences in the bulk velocities and other features of the plasma dynamics of these two configurations are stressed. Finally the numerical results are compared to experimental RFP data.

Most of the more technical details are given in the appendices. The dimensionless quantities, that should be introduced in the MHD simulations to reproduce the physical values of a tokamak like JET, are shortly discussed in Appendix E. Also, in Appendix F can we find a short report on preliminary simulations of a plasma dynamo experiment.

II NUMERICAL METHOD AND VALIDATION

INTRODUCTION

Magnetohydrodynamics is the discipline that studies the interaction between conducting fluids and magnetic fields. Depending on the topology and the intensity of the magnetic field, as well as on the values of the kinematic viscosity and the magnetic diffusivity, numerous different flow behaviors can be observed. As a function of the viscosity the flows can vary from a laminar to a highly turbulent state and the magnetic diffusivity allows to change the dynamics from a highly diffusive transport to an almost frozen-in advection of the magnetic field. Therefore, in the turbulent state, even in the statistically homogeneous case, a large range of dynamically active scales can be observed. Most of the interesting applications of MHD are however not statistically homogeneous due to the presence of solid walls. Examples are the planetary dynamo mechanism, magnetically confined fusion plasmas and industrial applications involving liquid metals [START_REF] Davidson | An introduction to magnetohydrodynamics[END_REF][START_REF] Suplee | The plasma universe[END_REF]. In order to study these phenomena, either experiments need to be carried out, or a set of nonlinear differential equations must be solved, combined with adequately chosen boundary conditions. In most cases, these equations cannot be solved analytically, so that numerical integration is needed in order to describe the dynamics. A wide range of MHD solvers have been developed over the last decades and an exhaustive listing is beyond the scope of this chapter. Nevertheless we will mention several approaches which are used to compute MHD in wall-bounded geometry. MHD codes applied to tokamak geometry are proposed in [START_REF] Luciani | XTOR-2F: a fully implicit Newton-Krylov solver applied to nonlinear 3D extended MHD in tokamaks[END_REF][START_REF] Huysmans | Non-linear MHD simulations of edge localized modes (ELMs)[END_REF] and to Reversed Field Pinch geometry in [START_REF] Ón | An optimal, parallel, fully implicit Newton-Krylov solver for three-dimensional viscoresistive magnetohydrodynamics[END_REF]. A more complete review of MHD solvers developed to compute fusion-plasma-related flows is given in [START_REF] Jardin | Review of implicit methods for the magnetohydrodynamic description of magnetically confined plasmas[END_REF]. Solvers aiming at a description of the dynamo effect are, for example, given in [START_REF] Guermond | An interior penalty Galerkin method for the MHD equations in heterogeneous domains[END_REF][START_REF] Iskakov | An integro-differential formulation for magnetic induction in bounded domains: boundary element-finite volume method[END_REF][START_REF] Reuter | A parallel implementation of an MHD code for the simulation of mechanically driven, turbulent dynamos in spherical geometry[END_REF] and computations investigating the magnetorotational instability in bounded domain were reported by Rüdiger and Shalybkov [START_REF] Shalybkov | Stability of axisymmetric Taylor-Couette flow in hydromagnetics[END_REF], Gissinger et al. [START_REF] Gissinger | Instabilities in magnetized spherical Couette flow[END_REF] and Willis and Barenghi [START_REF] Willis | Hydromagnetic Taylor-Couette flow: numerical formulation and comparison with experiment[END_REF]. An early numerical approach to study MHD in cylindrical geometry was proposed and validated by Shan et al. [START_REF] Shan | Nonlinear magnetohydrodynamics by Galerkin-method computation[END_REF] and more recently applied to spherical geometry by Mininni et al. [START_REF] Mininni | Magnetohydrodynamic activity inside a sphere[END_REF]. Most of these solvers are adapted to either a single geometry or a particular application. To change the geometry substantial effort must be invested to adapt the numerical mesh or to change the basis functions used in the numerical algorithm.

We present in this chapter the implementation and validation of the volume penalization method applied to magnetohydrodynamics. The strength of this approach is the high flexibility in the geometry and the ease of implementation. The volume penalization method is an immersed boundary method [START_REF] Peskin | The immersed boundary method[END_REF], in which both the fluid region and the confining boundaries are part of the same computational domain. The influence of the boundaries is then modeled by adding a force or drag term to the dynamical equations in the part of the domain in which the boundaries are to be present. In the volume-penalization method, for the hydrodynamic case, the solid bodies are modeled as porous media whose permeability tends to zero. This so called Navier-Stokes/Brinkman model, where the penalization source term in the momentum equation corresponds to the Darcy drag, was first proposed by Arquis and Caltagirone [START_REF] Caltagirone | Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide-milieux poreux: application à la convection naturelle[END_REF] in the context of the natural convection flow inside a cavity with porous walls. It was then generalized to study fluid -porous wallsolid boundary systems [START_REF] Caltagirone | New graphical and computational architecture concept for numerical simulation on supercomputers[END_REF]. In addition to being physically motivated, this model is mathematically justified, since Angot et al. [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] rigorously proved that the method converges to the Navier-Stokes equations combined with no-slip boundaries, when the porosity in the part of the domain corresponding to the boundaries is taken infinitesimally small. A first use of the method in combination with a pseudo-spectral Navier-Stokes solver was reported in [START_REF] Schneider | Numerical simulation of the transient flow behaviour in chemical reactors using a penalisation method[END_REF]. An extensive validation of the method for three dimensional fixed and moving boundaries is reported in [START_REF] Schneider | A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles[END_REF].

The strength of the combination of a pseudo-spectral Navier-Stokes solver with the penalization method is the compromise between accuracy and ease of implementation while retaining a great flexibility in the choice of the geometry of the boundaries. The method has been used to study two-dimensional MHD [START_REF] Bos | Rapid generation of angular momentum in bounded magnetized plasma[END_REF][START_REF] Neffaa | The decay of magnetohydrodynamic turbulence in a confined domain[END_REF][START_REF] Bos | Self-organization and symmetrybreaking in two-dimensional plasma turbulence[END_REF] which allowed to compare square, circular and periodic boundaries using the same Cartesian grid and numerical method. In the present chapter we will present a detailed validation of the method for two-and three-dimensional confined hydrodynamic and MHD flows.

The chapter is organized as follows. In section 2 -3, we expose the physical model and its numerical discretization. Section 4 assesses the parallel performance of our implementation. As a first validation, section 5 presents two-dimensional kinematic and magnetic test cases together with a comparison to analytical results. Validation of the threedimensional periodic MHD calculations is exposed in section 6.1. Subsequently Ohmic decay is considered in section 6.2. Sections 6.3 and 6.4 compare our three-dimensional results for 3D Taylor-Couette hydrodynamic and MHD flows with those available in the literature, as further validation. Section 6.5 reports on the flow induced in a conducting fluid by the presence of an imposed helical magnetic field.

MHD EQUATIONS

The media we study are isothermal, incompressible and we consider constant permeability µ, permittivity ε and conductivity σ of the material. The MHD equations for this case are the Navier-Stokes equation (including the Lorentz force) and the induction equation (that combines Ohm's law, Faraday's equation and Ampère's law). Introducing conventional normalization of the velocity by the Alfvén velocity C a = B 0 / √ ρ µ, a reference magnetic field B 0 and a conveniently chosen lengthscale L, the normalized equations read,

∂ u ∂t -ν∇ 2 u = -∇Π + u × ω + j × B, (2.1)
and

∂ B ∂t -λ ∇ 2 B = ∇ × [u × B] , (2.2)
where ν is the dimensionless kinematic viscosity (inverse of the Reynolds number), λ the dimensionless magnetic diffusivity (inverse of the magnetic Reynolds number) and ρ = 1 is the density. The vorticity ω and current density j are given by

ω = ∇ × u (2.3) j = ∇ × B, (2.4)
and Π = P + 1 2 u 2 is the modified pressure. The velocity field u is considered incompressible and the magnetic field B divergence free, ∇ • u = 0, (2.5) ∇ • B = 0.

(2.6)

To complete the problem one needs to specify the initial and the boundary conditions corresponding to the physical system that we are interested in. In particular the boundary conditions will be discussed in more detail in the rest of this chapter.

THE NUMERICAL CODE

The penalization method was introduced into two independently developed pseudo spectral MHD solvers, the L-code from Lyon and the M-code from Marseille. The cross-check of the results obtained by the two codes allowed a careful debugging and implementation of the method.

PSEUDO-SPECTRAL DISCRETIZATION

A classical Fourier pseudo-spectral method is used for the spatial discretization of a cubic periodic domain Ω of size 2π [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF]. The physical size of the domain can be modified rescaling the box by multiplying by L x , L y and L z . Spatial derivatives are evaluated in Fourier space and multiplications are computed in physical space. In the following we denote the Fourier Transform by the symbol ˆor F { }. All fields are represented as truncated Fourier series and here we show this representation explicitly for the velocity:

u(x,t) = N x /2-1 ∑ k x =-N x /2 N y /2-1 ∑ k y =-N y /2 N z /2-1 ∑ k z =-N z /2 u(k,t)e ik•x , (2.7) u(k,t) = 8π 3 N x N y N z N x -1 ∑ n x =0 N y -1 ∑ n y =0 N z -1 ∑ n z =0 u(x n ,t)e -ik•x n , (2.8)
with the wave vector k

= (k x , k y , k z ) where -N x /2 ≤ k x ≤ N x /2 -1, -N y /2 ≤ k y ≤ N y /2 - 1, -N z /2 ≤ k z ≤ N z /2 -1, x n = (n x 2π/N x , n y 2π/N y , n z 2π/N z ) ∈ [0, 2π] 3 with n x = 0, ..., N x -1, n y = 0, ..., N y -1, n z = 0, ..., N z -1.
The number of grid points in x, y and z-direction, N x , N y and N z , respectively, can be adapted to obtain the accuracy needed in the different directions. To avoid aliasing errors, i.e., the production of small scales due to nonlinear terms which are not resolved on the grid, the velocity and magnetic fields are dealiased at each time step by truncating its Fourier coefficients using the 2/3 rule [START_REF] Canuto | Spectral Methods in Fluid Dynamics[END_REF]. For the transformation between physical and Fourier space two different Fourier transforms were used in the two codes, firstly the P3DFFT routine, based on the FFTW library, secondly the JMFFT library. Both Fourier Transforms have an order of complexity of N log 2 N with N = N x N y N z .

The pressure term can be eliminated from the equations in spectral space by using the incompressibility condition of the medium. This introduces the projector P i j = δ i j - k i k j /k 2 in front of the nonlinear term. The Eqs. (2.1) and (2.2) in spectral space are then written:

∂ ûi ∂t + νk 2 ûi = P i j u × ω + j × B j , (2.9)
∂ Bi ∂t + λ k 2 Bi = ik × (u × B) i , (2.10)
where

k 2 = |k 2 |.

PENALIZATION METHOD

The volume penalization method is based on the idea of modeling solid bodies as porous media whose permeability tends to zero. The flow is considered in a domain in which both fluid and solid domains are embedded. The difference between the fluid and solid subdomain is the permeability. In the fluid domain the permeability is infinite and in the solid domain the permeability tends to zero. The method allows to consider an arbitrary shape and number of obstacles. The equations are modified by adding the penalization term:

∂ u ∂t = u × ω -∇Π + ν∇ 2 u + j × B - χ η (u -u wall ) (2.11) ∂ B ∂t = ∇ × (u × B) + λ ∇ 2 B - χ η (B -B wall ), (2.12)
with u wall and B wall the imposed values of the velocity and magnetic field in the solid domain and η is the permeability of the solid domain, i.e., the penalization parameter, which could be different for each equation, and χ(x,t) the mask function (see Fig. 2.1):

χ(x,t) = 0 for x ∈ Ω f , the fluid domain 1 for x ∈ Ω s , the solid domain.

(2.13)

Note that B wall can be freely chosen, we are not obliged to penalize all components.

For instance choosing B wall = B with B the component of B parallel to the wall, only penalizes the normal component and leaves the parallel component free. According to Eqs. (2.11) and (2.12), the flow is governed by the Navier-Stokes and induction equations in Ω f , and by Darcy's law in Ω s for small η. As mentioned in the introduction, the convergence of the velocity of the penalized equation in the limit of vanishing η to the solution given by the Navier-Stokes equations with no-slip boundary conditions was rigorously proven by Angot et al. [START_REF] Angot | A penalization method to take into account obstacles in incompressible viscous flows[END_REF] for fixed obstacles. The estimates were then refined by Carbou and Fabrie [START_REF] Carbou | Boundary layer for a penalization method for viscous incompressible flow[END_REF], who demonstrated that the solution of the penalized equations converges in the L 2 -norm with √ η towards the solution of the non-penalized equations with Dirichlet boundary conditions. Similar results are anticipated for the induction equation. To use a pseudo-spectral solver we need to Fourier-transform Eqs. (2.11) and (2.12) and we obtain

Ωs

∂ ûi ∂t + νk 2 ûi = P i j u × ω + j × B -F χ η (u -u wall ) j . (2.14)
We stress here that the Riesz projector P i j is also applied to the penalization term. Indeed, this form straightforwardly appears in the Fourier-transformed equations when the pressure is eliminated by solving a Laplace equation. The fact that the projector also acts on the penalization term is important to ensure incompressibility, since the penalization term is not necessarily divergence free at the fluid-solid boundary. This is also the case for the magnetic field. Due to the penalization term, the magnetic field is no longer divergence free. One way to cure this is to add an auxiliary pressure to the magnetic field

∂ B ∂t -λ ∇ 2 B = ∇ × [u × B] -∇Ξ Auxiliary pressure - χ η (B -B wall ) Penalization term . (2.15)
In the absence of boundaries in the domain, this pressure gradient would be equal to zero, as can be directly seen by solving a Laplace equation for Ξ and using the solenoidality constraint, Eq. (2.6). Indeed the (curl)-term ik × (u × B) is necessarily divergence free.

Eliminating the pressure from Eq. (2.15), we find for the Fourier-transformed equation for the magnetic field,

∂ Bi ∂t + λ k 2 Bi = P i j ik × (u × B) -F χ η (B -B wall ) j , (2.16)
which guarantees the incompressibility. The penalization method allows for a simple implementation of complex geometries, since to change the shape of the walls, one only needs to redefine the mask function. This is a huge advantage, because almost no effort is required to modify the shape of the flow geometry during an investigation and arbitrarily complex shapes can be considered. Several limitations should however be mentioned. First, the dynamical equations are solved in both the fluid domain and the penalized domain, so that, if the penalized domain is large, an important part of the numerical resources is used to compute the dynamically unimportant flow inside the walls (see for example section 6.4). Second, here no mesh refinement near the wall is used, so that, to capture small scale dynamics near boundaries, one needs to globally increase the resolution. Third, the boundary conditions imposed by the numerical method need to be satisfied in the computational domain. Let us explain this latter point in more detail.

Intrinsically, the boundary conditions of the Fourier pseudo-spectral solver are periodic in the three directions. Thus in the computational domain this periodicity must be satisfied. This imposes certain constraints on the geometries and especially on the boundary conditions that can be considered. If in the geometry sketched in Fig. 2.1 the solid domain corresponds to no-slip walls, i.e. u wall = 0, the periodicity condition is met automatically. However, if the outer-walls are chosen to move in solid-body rotation anticlockwise, the left border of the domain will move downwards whereas the right border will move upwards. In that case the periodicity condition is not satisfied. A solution to this problem is to add a third, unpenalized, subdomain outside the walls which will allow to respect the periodic boundary conditions of the computational domain. This solution is sketched in Fig. 2.2 (left). However it has an inconvenience which we will describe below.

Discontinuities in the velocity field, or in its gradients, are a source of Gibbs oscillations. These oscillations are an unavoidable feature in the present approach, and as long as their amplitude is small compared to physical effects, they do not constitute a serious problem in most cases. When discontinuities become strong, the Gibbs-oscillations also increase in size. Considering Fig. 2.2 (left), it is clear that if the solid domain turns and the outer fluid domain is very small, the velocity gradient becomes strong in the outer fluid domain and the discontinuity of the velocity gradient will become large on the interface between Ω s and the outer Ω f . Gibbs oscillations might get strong in this case. One solution is the following: instead of imposing in Ω s a solid body rotation, we impose a velocity profile which gives the correct boundary condition at the solid-fluid interface, and which smoothly tends to zero towards the edges of the computational domain Fig. 2.2 (right). The latter solution is slightly more complicated since an interpolation needs to be computed, using a Hermite's interpolating polynomial for instance. Its advantage is that the continuity of the solution and its derivatives between the boundary value and zero value can be imposed in a smooth way, which yields an improved order of convergence of the algorithm, as we will see in section 5 (an analytical analysis can be found in Appendix A).

Another drawback of the penalization method is that it is not yet possible to impose inhomogeneous Neumann conditions at the boundaries using a Fourier spectral code. In Appendix B a one-dimensional penalization method for taking into account non-homogeneous Neumann boundary conditions is presented. The lack of a threedimensional implementation makes it not yet possible to impose arbitrary values of the velocity gradient or magnetic gradient, for example, to impose the vorticity and the current density at the walls. This would in particular be important to model the influence of solid boundaries with arbitrary magnetic properties on the magnetic field generated in the fluid. A recent investigation by Kadoch et al. [START_REF] Kadoch | A volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles[END_REF] presents a technique for implementing homogeneous Neumann conditions using the penalization method with a spectral method. The extension to three-dimensional inhomogeneous Neumann conditions is an important perspective for further research.

TIME-DISCRETIZATION

Two different implementations of the penalization method in the time-advancing scheme will be compared. The first is explicit and constrains the penalization parameter η to be bigger than the time step ∆t to avoid numerical instabilities. The second is a semi-implicit implementation that allows the penalization parameter to be independent of the time step. Second and third order time schemes are used. In all approaches an exact integration of the viscous and magnetic diffusion term is used. In the following two sections these different methods are detailed.

EXPLICIT IMPLEMENTATION OF THE PENALIZATION TERM

In this section we detail the time integration of the equations using an explicit treatment of the nonlinear and penalization terms. It must be noted that along with its simplicity and robustness this approach has a drawback: the explicit treatment of the penalization term imposes a stability condition, in addition to the usual CFL condition. An analytical analysis of the magnetic part of the method has yet to be done to check if it adds another stability condition. Up to now no problem occurred if the same stability criterions were used for the velocity and the magnetic field.

The basic time-stepping schemes that are implemented are an adaptive second and third order Adams-Bashforth method (denoted by AB2 and AB3 respectively). The use of these schemes fits well into our general concept of compromise between the ease of implementation and computational efficiency. Exact integration of the diffusion term is feasible because the Laplace operator is diagonal in Fourier space and hence no linear system has to be solved. It improves stability of the scheme, avoiding the stability condition ∆t < ∆x 2 /ν. The remaining terms are discretized explicitly to avoid the solution of nonlinear equations, however it implies a CFL condition on the time step size ∆t and also a condition due to the explicit discretization of the penalization term , i.e., ∆t < η for AB2 and ∆t < 6 11 η for AB3, as linear stability analysis shows [START_REF] Schneider | A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles[END_REF]. For illustration, the equations will be given for the case of the velocity only. The discretization of the magnetic field is handled analogously, the only difference is the exact form of the nonlinear and penalization terms. First the Navier-Stokes equation is rewritten in the form of a nonlinear evolution equation and transformed into Fourier space,

∂ t u -ν∇ 2 u = N(u)
(2.17)

∂ t u + νk 2 u = N( u). (2.18)
For the initial condition u (k,t n ), the above equation has the following solution

u (k,t n+1 ) = e -ν∆t n+1 k 2 u (k,t n ) + t n+1 t n e -ν(t n+1 -s)k 2 N ( u (k, s)) ds, ( 2.19)
which can be discretized

AB2 -→ u (k, t n+1 ) =e -ν∆t n+1 k 2 u (k,t n ) + β 10 N n + β 11 e -ν∆t n k 2 N n-1 AB3 -→ u (k, t n+1 ) =e -ν∆t n+1 k 2 u (k,t n ) + β 20 N n +e -ν∆t n k 2 β 21 N n-1 + β 22 e -ν∆t n-1 k 2 N n-2 , (2.20)
with N n denoting the value of the nonlinear term at the time instant t n , the second order Adams-Bashforth coefficients

β 10 = 1 2 ∆t n+1 ∆t n (∆t n+1 + 2∆t n ) β 11 = - 1 2 ∆t 2 n+1 ∆t n , (2.21)
and the third order Adams-Bashforth coefficients

β 20 = ∆t n+1 (2∆t 2 n+1 + 6∆t n ∆t n+1 + 3∆t n-1 ∆t n+1 + 6∆t 2 n + 6∆t n-1 ∆t n ) 6∆t n (∆t n + ∆t n-1 ) β 21 = -∆t 2 n+1 (2∆t n+1 + 3∆t n + 3∆t n-1 ) 6∆t n-1 ∆t n β 22 = ∆t 2 n+1 (2∆t n+1 + 3∆t n ) 6∆t n-1 (∆t n + ∆t n-1 ) , (2.22)
where ∆t n = t nt n-1 [START_REF] Schneider | Numerical simulation of the transient flow behaviour in chemical reactors using a penalisation method[END_REF]. For start-up a first order scheme is used, as two time steps are required to start a second-order scheme. Similarly a first order and a second order scheme are used to start the third order scheme. The time step size control is based on the CFL stability limit of the explicit discretization of the nonlinear term, with addition of the stability criterion due to the penalization. Therefore, at each time step t n , the maximal point-wise velocity is computed and the new time step is given by ∆t n+1 = C∆x/U max where C < 1 is the CFL constant and ∆x is the minimal spatial grid size. Moreover, the time step has to verify the condition ∆t n+1 < η (AB2) or ∆t n+1 < 6 11 η (AB3) due to the presence of the penalization term. The same method is applied to the magnetic field and the time step is chosen to be small enough to verify the stability criteria of both the magnetic field and velocity field discretization.

SEMI-IMPLICIT IMPLEMENTATION

As noted in the previous section, the stability condition for a third order time scheme constrains ∆t < 6 11 η. To avoid this limitation, an implicit implementation was introduced by Kolomenskiy and Schneider [START_REF] Schneider | A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles[END_REF] for Burgers' equation and extended to Navier-Stokes' equation by Jause-Labert et al. [START_REF] Jause-Labert | Numerical validation of the volume penalization method in three-dimensional pseudo-spectral simulations[END_REF]. In this case the penalization term is evaluated at the time step t n+1 . The penalization is thus no longer treated together with the nonlinear term. The diffusion term, as in the explicit method, is exactly integrated.

This method is more time-expensive because two additional Fourier transforms are required. In addition to the projection of the nonlinear term, a second projection is realized (that includes the penalization term at t n+1 ) to ensure the solenoidal nature of the two considered fields. The fact that the time step can be adaptive (taking into account the CFL condition) makes this technique more suitable for unsteady simulations. The penalization term is introduced using a first order time scheme, which does not influence negatively the precision as long as boundaries are fixed.

The magnetic equations being handled analogously we present the new time scheme for the velocity field

ûi (k,t n+1 ) = P i j    F F -1 {Q n i } + ∆t η χ u wall i (x,t n+1 ) 1 + ∆t η χ j    . (2.23)
The third order Adam-Bashforth scheme is retained for the nonlinear terms in this formulation

Q n i = e -ν∆t n+1 k 2 u (k,t n ) + β 20 N n + e -ν∆t n k 2 β 21 N n-1 + β 22 e -ν∆t n-1 k 2 N n-2 . (2.24)
This numerical scheme for a penalization parameter η sufficiently small (η << ∆t) converges towards an explicit modified scheme where the time step ∆t replaces the penalization parameter η and where the nonlinear term vanishes in the solid region. We therefore call this method semi-implicit. This is further explained in Appendix C. For this case with very small penalization parameter, the permeability of the solid media is given by the value of the time step. The asymptotic convergence of the porous boundaries towards a solid wall, if η is sufficiently small, is in that case limited by the value η effective ≈ ∆t.

PARALLEL PERFORMANCES

The numerical code is parallelized using MPI libraries. The parallelization performances are evaluated on the calculator Vargas of the French high performance computing center IDRIS. The performances are estimated for a hydrodynamic and MHD calculation for three different grid resolutions. The test cases used are a three-dimensional Taylor-Couette flow for the hydrodynamic calculation (see section 6.3) and a MHD Taylor-Couette for the MHD calculation (see section 6.4). The comparison is made to ensure that the extra MPI exchanges and cache memory effects of the MHD calculations do not produce important slowdowns. We find that the MHD calculation is roughly twice as long as the hydrodynamic one. This is expected as we use for a MHD run the double number of fast Fourier transforms at each iteration. Note that per iteration, using a semi-implicit implementation of the penalization term, we need 12 Fourier transforms for each vector field. The results are shown in table 2.1 and in Fig. 2 The time per iteration per processor scales with the number of processors following an approximate power-law dependence close to the ideal scaling law, 1/(number of processors).

TWO-DIMENSIONAL VALIDATION

In this section we present a purely hydrodynamic test-case, the two-dimensional Taylor-Couette flow and a purely magnetic test-case, the z-pinch configuration. For both cases analytical solutions are known, which allows a careful convergence study and which allows to check different ways to introduce the boundary conditions.

TWO-DIMENSIONAL TAYLOR-COUETTE FLOW

We consider the classical two-dimensional hydrodynamic problem of a flow between two coaxial rotating cylinders (e.g., Taylor [START_REF] Taylor | Stability of a viscous liquid contained between two rotating cylinders[END_REF]). The inner cylinder rotates at constant speed, while the outer cylinder is kept at rest. The steady flow solution of the problem is

U θ (r) = Ω 2 R 2 2 -Ω 1 R 2 1 R 2 2 -R 2 1 r + (Ω 1 -Ω 2 )R 2 1 R 2 2 R 2 2 -R 2 1 1 r , (2.25)
where Ω (1,2) are the angular velocities of the cylinders, R (1,2) the radii of the cylinders and r the cylindrical coordinate (see Fig. 2.4). The relative L 2 error in the fluid domain || f numericalf analytical || L 2 /|| f analytical || L 2 with f being the considered field, is calculated for different penalization parameters η and number of grid points N [START_REF] Mininni | Magnetohydrodynamic activity inside a sphere[END_REF]128,256,512), in one direction with N = N x = N y and N z = 4.

R 1 R 2  1  2
As mentioned above, the present calculation is entirely hydrodynamic. The simulations are carried out until a steady state is obtained, so that the error is independent of the time discretization. A cubic domain with size-length 2π is considered, the time step is fixed to a value ∆t = 5.10 -5 and the kinematic viscosity ν = 1. The radii R 1 , R 2 are (0.32π,0.82π) respectively. At t = 0 the fluid domain is at rest and the inner-cylinder is set into movement with a fixed velocity U θ (R 1 ) = 1 while the velocity U θ (R 2 ) is kept equal to zero. The runs are stopped when the time t max = 5 is reached. At this time instant, the difference in the kinetic energy between two iterations is less than 10 -9 (for a kinetic energy of order unity), which indicates that a steady state is satisfactorily achieved.

The velocity profile imposed in the mask is chosen in two different ways, corresponding to the discussion in section 3.2. In the first case, the velocity in the inner cylinder is straightforwardly set to a solid-body rotation, U θ = Ω 1 r, in the inner cylinder and U θ = 0 in the outer cylinder. This is the most obvious choice. The component U r is set to zero at the boundaries. The velocity field is hereby continuous, but there exists a discontinuity of the velocity field derivative at the fluid boundaries (which is also the case in the real, physical situation). The error evolution with the penalization parameter and the convergence of the error with the resolution are shown in Figs. 2.5 and 2.6. For these calculations the expected √ η convergence order [START_REF] Carbou | Boundary layer for a penalization method for viscous incompressible flow[END_REF] is found and the convergence is second order in space as a function of the resolution N, confirming the results in [START_REF] Schneider | A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles[END_REF]. We also observe a saturation of the error for large N, corresponding to the penalization error. Error norm L A second way to impose the velocity in the mask will now be described. In this particular test-case the analytical solution is known (Eq. (2.25)) and we can use this information to increase the precision of the method. As mentioned in section 3.2, Gibbs oscillations are created due to discontinuities in the fluid variables or their derivatives. The discontinuity in the velocity gradient field can here be removed by using a 4 th order Hermite interpolation at the boundaries at r = R 1 , R 2 . Any purely azimuthal, axisymmetric flow is solenoidal so we can freely choose the velocity U θ (r) in the mask, as long as it respects the correct boundary conditions at r = R 1 , r = R 2 and r = π, the latter condition being imposed by the periodicity of the pseudo-spectral method. The velocity U θ (r) in the mask is chosen such that velocity and derivative at the fluid-solid boundaries are continuous. Subsequently the velocity field is interpolated to decay smoothly from the analytical solution at r = R 1 , R 2 to zero at r = 0 and r = π respectively, using an interpolating Hermite polynomial. Any discontinuity on the derivative of the velocity field at the fluid boundaries is thus avoided and the Gibbs oscillations are hereby significantly reduced. In principle even higher order velocity derivatives could be smoothed in this way using higher order Hermite interpolation. Note that a similar regularization is used in [START_REF] Laizet | High-order compact schemes for incompressible flows: A simple and efficient method with quasi-spectral accuracy[END_REF] for the velocity field, where the imposed velocity in the mask is called an internal flow.

2 of u θ N η = 1 • 10 -4 η = 1 • 10 -5 N -2
The error as a function of the penalization parameter and the convergence of the error with the resolution are shown respectively on Figs. 2.7 and 2.8. The numerical error is only calculated in the fluid domain. It is observed that the convergence of the error with the penalization parameter is close to third order. An optimum for the penalization parameter depending on the resolution appears, when the gridscale becomes of order √ νη. At this scale the viscous term becomes of the order of the penalization term. In Figs. 2.7 and 2.8 results for both explicit and semi-implicit methods are presented.

The regularization of the boundary conditions using an interpolation clearly improves the numerical convergence of the solution with the penalization parameter. Also the convergence with the grid resolution is improved. If the Hermite interpolation is used, a fourth order convergence with N is found for both the explicit (dashed line) and semiimplicit (solid line) implementations (see Fig. 2.8). We recall that if no regularization of the velocity field is introduced, second order convergence is recovered (see Fig. Error norm L

2 of u θ N η = 1 • 10 -3 η = 5 • 10 -4 η = 1 • 10 -4 N -2 N -4

THE Z-PINCH

The second validation test is the reproduction of the z-pinch phenomenon, a well-known textbook example of a confined plasma situation [START_REF] Guthmann | Physique des Plasmas[END_REF]. This configuration is illustrated in Fig. 2.9. Two ideal electrodes drive an axial current in the z-direction producing a purely azimuthal magnetic field (in the θ -direction). The current density in the z-direction, which together with the induced azimuthal magnetic field yield a radially pinching Lorentz force, is the motivation for the name z-pinch. In this configuration we set the velocity to zero so the code is entirely magnetic. We impose the boundary conditions B θ = B C and B r = 0 at r = R 1 the radius of the fluid domain. The component B z is not penalized and can freely evolve.

Periodic conditions are set in the axial direction. In this configuration the governing equations reduce to

∂ t B = λ ∇ 2 B (2.26)
In cylindrical coordinates, the steady state solution is a linear evolution of the azimuthal magnetic field

B θ (r) = B C r R 1 .
(2.27)

The quantity B C /R 1 is linked to the constant axial current density, using Ampère's law,

j z = 2B C R 1 .
(2.28)

The computational domain is chosen similar to Fig. 2.2 (left). The mask is chosen to be annular, leaving the outer domain free to adapt to the periodic boundary conditions of the computational domain. If a uniform azimuthal magnetic field is imposed inside the mask, the discontinuity in the profile of the radial derivative of the azimuthal magnetic field at the boundary causes Gibbs' oscillations in the current density, analogously to what was observed in the Taylor-Couette case. To avoid this, a linear profile of B θ (r) = B C r R 1 is imposed inside the mask to ensure a continuity with the analytical solution. This feature eliminates the discontinuity at r = R 1 in the derivative of the magnetic field and thus greatly reduces the oscillations for j z . The convergence of the method can be further enhanced by using a Hermite polynomial to smoothly interpolate the magnetic field in the mask to zero at the outer boundaries of the computational domain. With this method, the continuity of the derivative of the magnetic field is assured through the entire domain (see Fig. 2.2 (right)). For these simulations the number of grid points are the same as for the Taylor-Couette case (N ∈ {64, 128, 256, 512}). The parameters are a cubic domain with size-length 2π with magnetic diffusivity λ = 1, the time step is fixed to ∆t = 5.10 -5 and t max = 5. For t = t max the difference in the magnetic energy between two iterations is less than 10 -9 so we have reached the steady state. The inner radius of the annulus is R 1 = 0.65π and the outer radius is R 2 = 0.78π. If the Hermite polynomial is used, the radius where it reaches the value 0 is R 3 = 0.94π. The boundary condition B C = 1.

In Fig. 2.10 the convergence of the method is shown as a function of the penalization parameter η and the resolution N. An improved convergence, proportional to η 4 or N -4 is observed. The relative importance of the smoothing of the magnetic field in the mask is illustrated in Fig. 2.10 where the results of these computations with and without regularization are compared. Without the regularization technique, the convergence reduces to second order in N and order √ η for the penalization parameter. Error norm L

2 of B θ N η = 5 • 10 -4 -Hermite η = 1 • 10 -4 -Hermite η = 5 • 10 -4 -No Hermite η = 1 • 10 -4 -No Hermite N -2 N -4

ASSESSMENT OF THE REGULARIZATION METHOD TO ENHANCE THE PERFORMANCE OF THE PENALIZATION METHOD

For these two-dimensional test cases, either Taylor-Couette or z-pinch, the error of convergence as a function of the penalization parameter and the resolution are determined. The regularization of the different fields in the solid domain (or mask) allows an enhancement of the accuracy of the numerical solution in the fluid region. The fact that the continuity of the velocity derivative field inside the solid domain affects the error inside the fluid domain is an intrinsic feature of the pseudo-spectral method. Such methods use periodic trial functions and the Gibbs oscillations introduced in one point of the domain decay only inversely proportional to the distance from the discontinuity. The Hermite interpolation method regularizes and yields fields which are C 1 in the whole domain (see Appendix A for an analytical analysis of the Hermite regularization). Gibbs oscillations are thus reduced and consequently the numerical errors are also decreased considerably. The limitation of this method is that an analytical solution must be known, or a baseflow, which is not far from the expected developed flow. Without such regularization the convergence is reduced to second order in resolution, which can be sufficient for many applications, as illustrated in the following sections.

THREE-DIMENSIONAL VALIDATION

In this section the code will be validated by considering three-dimensional test-cases. First a periodic MHD case is considered, without using the penalization method, subsequently the magnetic part of the code is validated by studying the Ohmic decay in a cylindrical cavity. Then the three-dimensional Taylor-Couette flow is studied with and without the presence of a magnetic field and to conclude we investigate the instabilities in a cylinder with helical magnetic boundary conditions.

PERIODIC MHD VALIDATION

To validate the capacity of the numerical code to simulate the three-dimensional nonlinear MHD equations, we reproduce first a classical test-case with periodic boundary conditions. This case is the generalization of the Orszag-Tang vortex to three dimensions.

The results are compared with those of Mininni et al. [START_REF] Mininni | Small-Scale Structures in Three-Dimensional Magnetohydrodynamic Turbulence[END_REF].

The initial condition used for the simulation is given analytically and yields:

u(x, y, z,t = 0) = [-2 sin(y), 2 sin(x), 0] for x, y, z ∈ [0, 2π] 3 (2.
max |j| = max j 2 x + j 2 y + j 2 z (2.31)
and the total dissipation rate is

ε(t) = ν ω 2 + λ j 2 , (2.32)
where ν and λ are respectively the kinematic viscosity and the magnetic diffusivity. Three runs are performed: the first with N 3 = 64 3 and ν = λ = 0.01, the second with N 3 = 128 3 and ν = λ = 0.005 and the third N 3 = 256 3 with ν = λ = 0.001. The evolution of the maximum of the current density and total dissipation rate in the domain are shown in Fig. 2.11. Results are compared using a second-and third-order Adams-Bashforth time advancing scheme described in section 3.3. Both schemes give the same results. The results agree well with the computations presented in reference [START_REF] Mininni | Small-Scale Structures in Three-Dimensional Magnetohydrodynamic Turbulence[END_REF]. The same exponential growth followed by a self similar growth ∼ t 3 is found in our calculations for the evolution of the maximum current density (see Fig. 2.11 (left)). With an increasing Reynolds number Re the maximum of |j| is also found to be reached at later Dissipations

Time ν = λ = 1 • 10 -3 ν = λ = 5 • 10 -3 ν = λ = 1 • 10 -2 Figure 2
.11: Comparison of the solution of the numerical code for the time evolution of the maximum of the current (left) and total dissipation rate (right) using second (dashed line) and third order time-schemes (solid line). The inset (left) shows the evolution at early times in lin-log units.

times. The evolution of the total dissipation rate (see Fig. 2.11 (right)) shows the same delay in the onset of the formation of small scales with increasing Re as exposed in the cited article. This test allows us to evaluate the full MHD code and validate the numerical results for relatively high Reynolds numbers (up to Re = 3000).

OHMIC DECAY IN A PERIODIC CYLINDER

In this test case we compute the evolution of the magnetic field in a periodic three dimensional cylinder coated with an insulator [START_REF] Guermond | A new Finite Element Method for magneto-dynamical problems: two-dimensional results[END_REF][START_REF] Laguerre | Induction effects in isolated axisymmetric conductors using a new finite element method[END_REF]. In the induction equation we set the velocity to zero so the equation for the magnetic field reduces to the diffusion equation:

∂ B ∂t = ∇ 2 B. (2.33)
We consider an axisymmetric case, z-independent and the magnetic field has no r component. This case is not completely three-dimensional since we use the three components of the magnetic field (B x , B y , B z ), but the derivatives are zero in the z direction. In cylindrical coordinates the set of equations is:

   ∂ B θ ∂t = ∂ 2 B θ ∂ r 2 + 1 r ∂ B θ ∂ r -B θ r 2 ∂ B z ∂t = ∂ 2 B z ∂ r 2 + 1 r ∂ B z ∂ r .
(2.34)

Using separable elementary solutions, the magnetic field can be written in the following form:

B θ (r,t) = A f (r)e -ω 2 θ t B z (r,t) = B g(r)e -ω 2 z t . (2.35)
Introducing the following change of parameter: s = ωr the system of equations writes:

s 2 θ f ′′ + s θ f ′ + (s θ 2 -1) f = 0 s 2 z g ′′ + s z g ′ + s 2 z g = 0.
(2.36)

The solutions of these equations are Bessel functions. Imposing at the radius R 0 of the cylinder the azimuthal and axial field vanishing the general solution is:

     B θ (r,t) = J 1 j 1 R 0 r e - j 1 R 0 2 t B z (r,t) = J 0 j 0 R 0 r e - j 0 R 0 2 t .
(2.37)

Here j 0 = 2.4048... and j 1 = 3.8314... are the first zeros of the Bessel functions J 0 and J 1 respectively.

In our simulation we set R 0 = 1. The decay rate is determined doing a least square fitting of the azimuthal and axial magnetic energy time evolutions (Figs. 2.12 and 2.13). In table 2.2 we present our results for the decay rate of the azimuthal and axial fields and we compare them to the analytical values (see Eq. (2.37)). All the simulations are done with 96 3 grid points, the penalization parameter η = 5 • 10 -4 and the computational domain size is L x = L y = L z = 2π. We calculate the error for different fixed time steps. In none of these calculations the regularization of the magnetic field in the solid region (mask) is used. For these calculations of the diffusion of a magnetic field in a periodic cylinder the decay rate of the azimuthal and axial component of the magnetic field agree quite well with the analytical values. We find the same relative error in the azimuthal and in the axial direction at each considered time discretization. This test-case yields a validation of the magnetic part of the numerical code and shows that the magnetic boundary conditions are well taken into account via the penalization method. Also we note that the time scheme is well implemented as it allows to recover the analytical decay rates for the considered components of the magnetic field with good accuracy. The numerical solution converges towards the analytical solution if the time step is decreased. Azimuthal

E k Time ∆t = 1 • 10 -2 ∆t = 1 • 10 -3 ∆t = 5 • 10 -4 ∆t = 1 • 10 -4
Figure 2.12: Azimuthal magnetic energy decay for different time steps in an periodic cylinder.

1.10 -6

1.10 -5

1.10 -4

1.10 -3

1.10 -2 1.10 -1 1 0 0.2 0.4 0.6 0.8 1 Axial E k Time ∆t = 1 • 10 -2 ∆t = 1 • 10 -3 ∆t = 5 • 10 -4 ∆t = 1 • 10 -4
Figure 2.13: Axial magnetic energy decay for different time steps in an periodic cylinder.

THREE-DIMENSIONAL TAYLOR-COUETTE FLOW

In this test case we aim to determine the critical Reynolds number for the first instability of the Taylor-Couette flow with periodic boundaries in the axial direction. Different values of the Reynolds number are explored with several calculations with a resolution of 128 3 grid-points, the penalization parameter η = 5 • 10 -4 and the computational domain size is L x = L y = 5π/2 and L z = 2π. To assess the influence of the Hermite interpolation method, in one of these calculations the regularization of the velocity in the solid region (mask) is used (for the case with Re = 120). The reference length scale is the gap between the inner and outer cylinder L = R ext -R int and the reference velocity is the inner rotation speed U = Ω int R int . The outer cylinder is fixed. We define the Reynolds number and also a radius ratio ζ and a aspect ratio Γ as follows:

Re = UL ν , ζ = R int R ext , Γ = L z L , (2.38)
where L z is the axial length. To be able to compare with the literature we take the same values as in [START_REF] Guermond | Nonlinear magnetohydrodynamics in axisymmetric heterogeneous domains using a Fourier/finite element technique and an interior penalty method[END_REF] for the dimensionless values, radius ratio and aspect ratio, ζ = 0.5 and Γ = 4. The base flow consists of an azimuthal velocity only, as in the two-dimensional case (section 5.1). The first Taylor-Couette instability is centrifugal and is characterized by vortices that appear and break the axial invariance. Velocities in the radial and axial directions thereby appear. To determine the critical Reynolds number we analyze the evolution of the axial kinetic energy. The analysis of the evolution of the axial kinetic energy allows us to assess the critical Reynolds number (when the instability is triggered). The value of the critical Reynolds is compared with a theoretical value of Re = 68.23 that has been determined by Chandrasekhar [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF]. We present in Fig. 2.14 the axial kinetic energy evolution for Reynolds numbers varying from Re = 65.7 to Re = 69.7. We start the simulations with a small random perturbation so the initial axial kinetic energy is non zero. The axial energy either grows or decays exponentially. The critical Reynolds number can be determined from Fig. 2.14. Increasing the Reynolds number from Re = 67.3 to Re = 68.1 the axial kinetic energy changes from decaying to increasing. A simple linear interpolation of the growth and decay rates (that are determined with a least square method fitting) gives the value of the critical Reynolds Re c ≈ 67.9. The estimated error compared with the theoretical result is ∼ 0.44%. The theoretical estimate is thus well approached by our numerical simulations. The Taylor vortices appear early in the simulation and they grow or decay in strength depending on the Reynolds number. Fig. 2.15 presents the axial kinetic energy evolution of a Taylor-Couette flow for Re = 120 comparing two simulations, one with Hermite regularization and one without. The difference between the two considered methods is the onset of the instability, it is slightly earlier if the Hermite regularization is used. The growth rates are similar, when the nonlinear saturation is reached, a steady state is obtained that almost coincides for both methods. In Figs. 2.16 to 2.19 we visualize the Taylor vortices in the steady state. Two pairs of counter-rotating vortices appear. The aspect ratio is Γ = 4 so that four vortices form. In Fig. 2.17 we distinguish the boundary layer, in which the azimuthal vorticity is contrary to the vorticity of the Taylor-vortices.

The flow structure of Figs. 2.16 to 2.19 can be compared with the one presented by Guermond et al. (Fig. 5 in [START_REF] Guermond | Nonlinear magnetohydrodynamics in axisymmetric heterogeneous domains using a Fourier/finite element technique and an interior penalty method[END_REF]). The same flow topology with four vortices is found. In Figs. 2.18 and 2.19 we note that there is a transport of azimuthal momentum by the radial flow. The azimuthal velocity isosurfaces are not axial invariant but they are dragged inand outwards by the radial flow. A positive radial velocity increases the azimuthal velocity near the inner cylinder. To compare quantitatively both simulations, in table 2.3 the maxima of the three velocity components at the steady state are compared with those of Guermond et al. [START_REF] Guermond | Nonlinear magnetohydrodynamics in axisymmetric heterogeneous domains using a Fourier/finite element technique and an interior penalty method[END_REF].

The azimuthal velocity is not exactly unity because with the penalization method the precise value at the boundary is not necessarily coinciding with the numerical grid as is the case in [START_REF] Guermond | Nonlinear magnetohydrodynamics in axisymmetric heterogeneous domains using a Fourier/finite element technique and an interior penalty method[END_REF]. Nevertheless the numerical method yields good agreement for all three components of the velocity with the results of the code described in [START_REF] Guermond | Nonlinear magnetohydrodynamics in axisymmetric heterogeneous domains using a Fourier/finite element technique and an interior penalty method[END_REF]. The three numerical methods yield similar relative variations, there is an improvement in the error with the regularization for the radial component of the velocity field. The improvement is limited because the tangent imposed with the Hermite polynomial is calculated using the Guermond et al. [START_REF] Guermond | Nonlinear magnetohydrodynamics in axisymmetric heterogeneous domains using a Fourier/finite element technique and an interior penalty method[END_REF] 3rd order semiimplicit 2nd order explicit 3rd order semiimplicit no regularization no regularization with regularization max u r 0.1935 0.19434 (∼ 0.43%) 0.1969 (∼ 1.75%) 0.19355 (∼ 0.03%) max u θ 1 0.99693 (∼ 0.31%) 0.9980 (∼ 0.20%) 0.99669 (∼ 0.33%) max u z 0.1454 0.14639 (∼ 0.68%) 0.1506 (∼ 3.57%) 0.14632 (∼ 0.63%) analytical base flow, which is known in this case. Since the development of the instability makes the flow change, discontinuities appear at the fluid-solid interface and the Gibbs oscillations can grow. The improvement can be substantial if the saturated state is not far from the calculated analytical base flow.

With this test-case the well known linear first instability threshold of the Taylor-Couette flow is found. Also the nonlinear saturation is comparable to what is reported in the literature. The same topology of the flow is observed and quantitatively similar velocity magnitudes at the steady state are established. This makes us confident about the accuracy of the method in taking into account centrifugal, pressure and nonlinear effects in a bounded domain.

MAGNETOHYDRODYNAMIC TAYLOR-COUETTE FLOW

We now extend our validation to the test-case of the instability of an axisymmetric MHD Taylor-Couette flow with periodic boundary conditions in the axial direction. In this case the instability studied in the previous section is modified due to the presence of an axial magnetic field. An imposed constant axial field B 0 is added to the magnetic field. It is well known that such an axial magnetic field has a significant stabilizing effect. This phenomenon of delay in the appearance of the first Taylor-Couette instability was found by Chandrasekhar [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF], confirmed by linear numerical calculations in [START_REF] Shalybkov | Stability of axisymmetric Taylor-Couette flow in hydromagnetics[END_REF] and by spectral numerical simulations [START_REF] Willis | Hydromagnetic Taylor-Couette flow: numerical formulation and comparison with experiment[END_REF][START_REF] Willis | Hydromagnetic Taylor-Couette flow: wavy modes[END_REF]. The fluid flow will try to bend the axial magnetic field lines but the restoring Lorentz force will prevent the fluid motion, which stabilizes the flow profile.

To assess our numerical code in this context, we evaluate the evolution of the axial kinetic energy as a function of the Hartmann number Ha, which measures the ratio between electrodynamic forces and viscous forces. The presence of a uniform magnetic field in the axial direction does not affect the profile of the stable azimuthal velocity that exists without the magnetic field, Eq. (2.25), which we choose as initial condition. The dimensionless geometric parameters, radius ratio ζ and aspect ratio Γ are the same as in the previous section. We introduce here the magnetic Prandtl number which is the ratio of viscosity and magnetic diffusivity. The dimensionless numbers describing the problem are then

Pr = ν λ , Re = UL ν , Ha = B 0 L µ 0 ρνλ , ζ = R int R ext , Γ = L z L , (2.39)
where µ 0 is the magnetic constant and ρ is the fluid density. The simulations are performed for Pr = 1, Re = 100, ζ = 0.5 and Γ = 4. The resolution used is N 3 = 128 3 , the penalization parameter η = 5 • 10 -4 and the computational domain size is L x = L y = 5π/2 and L z = 2π. For none of these calculations regularization of the velocity or magnetic field in the solid region (or mask) is used. The boundary conditions described in [START_REF] Shalybkov | Stability of axisymmetric Taylor-Couette flow in hydromagnetics[END_REF] are a fixed azimuthal velocity on the inner cylinder (U = Ω int R int = 1), no-slip on the outer cylinder and perfectly conducting walls, so that the normal magnetic field at the wall and the axial current density vanish (b r = 0 and j z = 0) [START_REF] Shalybkov | Stability of axisymmetric Taylor-Couette flow in hydromagnetics[END_REF]. With these parameters the Taylor-Couette flow is hydrodynamically unstable, as was observed in the previous section. With the penalization method we can impose the vanishing radial magnetic field but the current density is not constrained. The evolution of the axial kinetic energy varies as a function of the imposed magnetic field (or Hartmann number), for a fixed Reynolds number, as is shown in Fig. 2.20. These calculations allow us to determine the threshold for the instability. For Re = 100 the critical Hartmann number found is Ha c ≈ 7.9. Like in the previous section the threshold is determined by linear interpolation of the growth and decay rates. For Re = 100 the flow is hydrodynamically unstable and Taylor vortices should appear, but for Ha > 7.9, the instability is suppressed by the magnetic field. With this test-case we therefore found the well known stabilizing effect of an axial magnetic field on the Taylor-Couette flow. Qualitatively the flow behavior is very similar to what is found in other investigations. We were not able to quantitatively compare with the literature since our boundary conditions on the magnetic field are not the same as those considered in previous studies on magnetohydrodynamic Taylor-Couette flow and the parameter ζ was different. To compare with more accuracy our code with the literature we treat in the following section a test case with boundary conditions and geometry which are adapted to our numerical method.

FLOW INDUCED BY A HELICAL MAGNETIC FIELD

Shan, Montgomery and Chen [START_REF] Shan | Nonlinear magnetohydrodynamics by Galerkin-method computation[END_REF] studied numerically a conductive fluid confined in a periodic cylinder where an axial electric and magnetic field are imposed, which results in a helical magnetic field (see Fig. 2.22). They used a spectral code which decomposed the fields into Chandrasekhar-Kendall orthonormal eigenfunctions of the curl. In that study they discovered a transition between an axisymmetric state with a zero velocity to a laminar helical state where a dynamic equilibrium appears, i.e., a steady state with non zero velocity. The parameters chosen for the numerical study are selected to closely reproduce the simulations in [START_REF] Shan | Nonlinear magnetohydrodynamics by Galerkin-method computation[END_REF]. A fixed axial magnetic field B 0 = 4.5 is imposed and the fluid has a constant magnetic diffusivity and kinematic viscosity λ = ν = 0.045. The computational domain size is L x = L y = 0.8π and L z = 8. The typical length scale is the cylinder radius, L = R 0 = 1, the axial length is L z = 8R 0 , the resolution used for the simulations is N 3 = 128 3 grid-points and the penalization parameter η = 5 • 10 -4 . Three dimensionless numbers characterize the system: the Lundquist number (S), the Hartmann number (Ha) and the pinch ratio for the axisymmetric zero flow state (Θ 0 ):

S = C A L λ , Ha = B 0 L ρ µ 0 λ ν , Θ 0 = B θ B z , (2.40)
with C A the axial Alfvén velocity C A = B 0 / √ ρ µ 0 , B θ is the wall-averaged poloidal magnetic field and B z is the volume-averaged axial magnetic field.

The transition between states is determined as a function of one of these dimensionless numbers, the pinch ratio Θ 0 , which is varied by adjusting the imposed average axial electric field E 0 and which is directly linked to the poloidal magnetic field. The other dimensionless numbers are maintained constant, S = Ha = 100. For the parameters given above, the linear theory predicts a transition for an imposed electric field E 0 = 0.33 [START_REF] Shan | Nonlinear magnetohydrodynamics by Galerkin-method computation[END_REF][START_REF] Montgomery | Helical, dissipative, magnetohydrodynamic states with flow[END_REF].

To compare with Shan et al. we impose the same boundary conditions. The walls are treated as perfect conductors and are coated with a thin layer of insulator. Hereby both the radial magnetic and current density field vanish (B r = j r = 0). For the velocity field only the radial component vanishes at the wall. The penalization method is used to impose a vanishing normal component of the velocity (u • n = 0) without any regularization technique. Also an azimuthal magnetic field B 0 θ is imposed via the penalization term. In this case the vanishing radial current density ( j r = 0) is automatically satisfied because the azimuthal magnetic field does not generate a radial current density (the r-component of the curl of the imposed magnetic field is zero). The boundary conditions are thus satisfied.

The way the electric field is imposed in our simulations differs from the simulations by Shan et al.. In their simulations the electric field explicitly appears in the discretized equations, whereas in our case the electric field is indirectly imposed through the magnetic field at the wall (which can be related to the electric field using Stokes' theorem). This can lead to small differences in the transients, but is not expected to greatly affect the steady state solutions.

The azimuthal magnetic field B 0 θ (r) is imposed with the volume penalization method in the solid region using the regularization technique, like for the z-pinch case (see section 5.2). The azimuthal magnetic field increases linearly with r from r = R 0 (the fluid-solid frontier) to r = 0.34π and then smoothly tends to zero using a Hermite's interpolating polynomial from r = 0.34π to r = 0.385π. The periodicity of the computational region is hereby satisfied and the magnetic derivative of the base-field is continuous at the wall.

To validate the code we perform the same calculations done by Shan et al.. We vary the axial current density ( j z ) and we calculate E 0 when the simulation reaches a steady state using Ohm's law

E 0 = (-u × B) z + j z σ . (2.41)
We find (see Fig. 2.24) that the instability threshold between the axisymmetric and helical state is situated between E 0 = 0.302 and E 0 = 0.355 as found in [START_REF] Shan | Nonlinear magnetohydrodynamics by Galerkin-method computation[END_REF]. The kinetic energy starts to grow when the imposed electric current is E 0 = 0.355. The growth rate of the energy is calculated using Alfvén time units (t A = t num C A /L z ). A least-square fitting gives the growth rate of the kinetic energy as 0.54, corresponding to a growth rate of 0.27. This can be compared to the analytical value 0.279 calculated for an applied electric field E 0 = 0.35. Taking into account that our imposed magnetic field is slightly different (since the electric field is imposed indirectly in our case), the two different growth rates match in good order. In Fig. 2.23 the excited mode is visualized, which is a helical mode with azimuthal and axial mode numbers m = k z = 1, respectively, as in [START_REF] Shan | Nonlinear magnetohydrodynamics by Galerkin-method computation[END_REF]. Increasing the pinch ratio to E 0 = 0.402 the flow returns to its axisymmetric copperwire solution, which is also observed for certain values in [START_REF] Shan | Nonlinear magnetohydrodynamics by Galerkin-method computation[END_REF]. An explanation for this behavior is the shape of the instability curves in the Θ 0 -Ha plane. By increasing the pinch ratio, different (m, k z ) helical modes appear at a fixed Ha, but they can disappear by increasing Θ 0 to even higher values. This was investigated in [START_REF] Shan | Global searches of Hartmann-numberdependent stability boundaries[END_REF].

In the next figures, 2.25 to 2.28, we compare our different simulations with the ones performed by Shan et al.. We find quantitatively the same evolution of the average current density and the total dissipation rate (ε T = λ j 2 + ν ω 2 ). These quantities are time averaged during the dynamical steady state, since the flow becomes chaotic, if the pinch ratio (or E 0 ) is large. Some quantitative differences are observed in Figs. 2.25 and 2.26 for values around E 0 = 0.6. At that point both methods might not trigger exactly the same helical modes, since the electric field is imposed in a slightly different way. Both methods might therefore give results corresponding to different multi-mode states. When the pinch ratio is increased further, the deviations become smaller, as can be seen in Figs. 2.27 and 2.28.

This test-case allows us to validate the nonlinear MHD code with boundary conditions imposed on both the velocity and magnetic field. A linear analytical result is confirmed. The "multi-mode" and turbulent states that are observed in literature also appear in our simulations.

CONCLUDING REMARKS CONCERNING THE NUMER-ICAL METHOD

An extension and implementation of the penalization method into a pseudo spectral Fourier code solving the MHD equations is presented. This penalization method, which allows the introduction of obstacles and walls in the computational domain, is implemented in different ways with respect to the numerical scheme and definition of the fields within the solid domains.

The numerical code is validated by comparison with several test-cases and theory. First in two dimensions the convergence of the results towards an analytical solution by decreasing the penalization parameter and increasing the resolution are shown. The method converges faster than second order if a regularization technique in the solid domain is applied which removes the discontinuities in the derivatives of the velocity and magnetic field at the solid-fluid boundary.

Then in three dimensions the first instability threshold of the hydrodynamic Taylor-Couette flow is found with good accuracy. Also the nonlinear saturation of this flow is compared and validated with the literature. For the MHD Taylor-Couette flow the current inability to impose non-homogeneous Neumann boundary conditions using the penalization method, makes the comparison of our numerical results with the literature difficult. A more appropriate test-case to validate the three-dimensional implementation of the penalization method to compute MHD flows is the flow induced by a helical magnetic field. This case is correctly reproduced. The linear threshold of the transition between an axisymmetric and a helical state is found. Also the evolution of the average axial current and the total dissipation rate as a function of the average electric field are compared with the literature and are in good agreement.

All these test-cases allowed us to validate the numerical method to solve correctly the MHD equations in a confined domain. The limitations are the restricted magnetic boundary conditions that can be applied. Presently, the current density can only be imposed indirectly via the magnetic field. A modified volume penalization method that allows the introduction of non-homogeneous three-dimensional Neumann boundary conditions is currently being developed, preliminary one-dimensional results are presented in Appendix B.

III MAGNETOHYDRODYNAMICALLY GENERATED VELOCITIES IN CONFINED PLASMA 1 INTRODUCTION

Toroidal magnetic plasma confinement has been under investigation since the 1940s when it was recognized as a promising geometry for controlled thermonuclear fusion. Despite all the attention devoted to the idea, there are aspects of it that must be regarded as incomplete, even in theory. The difficulties in many cases reduce to the fact that there is no mathematical description of a magnetically active, dissipative plasma that is tractable, by use of even the fastest supercomputers. Time dependent electromagnetic fields combined with the particle kinetics of plasmas having the range of mass ratios represented among the various charges is simply a too large system to be susceptible to a complete treatment. Enormous simplifying assumptions have to be made to achieve any analytical/numerical progress. A common assumption has been that of an unstable ideal equilibrium whose numerous linear instabilities may reveal insight into the nonlinear dynamical behavior that is observed. It must be conceded that any description that is manageable at a detailed level will omit certain important features of a real plasma and at this stage it is to some extent a matter of taste as to which incomplete theoretical description is adopted for study.

In the following pages, we report the investigation of one such description: a voltage driven, dissipative magnetohydrodynamic (MHD) fluid with non-ideal toroidal boundaries. We omit some features that would be desirable and which seem reasonable to inject, at a later date, one at a time, into the numerical recipe we use. The principal unrealistic assumptions we make are those of uniform mass density and incompressibility, a scalar valued Newtonian viscosity, a scalar valued electrical conductivity, and the omission of a finite thermal conductivity (it will be seen that in effect an infinite thermal conductivity has been assumed, since no thermal effects are allowed to develop except those associated with the incompressible velocity field). Despite what appear to be these gross oversimplifications, what remains is at the very perimeter of what is computable if we intend to stay with arbitrary initial configurations which are not in equilibrium, and to follow through with enforcing viscous and resistive boundary conditions.

What is of particular interest is the spontaneous development of both toroidal and poloidal rotation of the bulk magnetofluid as a whole. It is not physically obvious that this should happen, even though it has been known for some time to occur in toroidal laboratory devices [START_REF] Rice | Inter-machine comparison of intrinsic toroidal rotation in tokamaks[END_REF]. The importance of non-zero velocities in the MHD description of toroidally confined plasma was realized by Pfirsch and Schlüter [START_REF] Schl Üter | Der Einfluss der elektrischen Leitfähigkeit auf das Gleichgewichtsverhalten von Plasmen niedrigen Drucks in Stelleratoren[END_REF], however without taking into account all the different terms in the force balance. We will take into account all these terms. The resulting flow pattern is presented here as a computational fact. The degree of the two types of rotation are seen to depend upon several things, such as the Reynolds-like dimensionless numbers assumed for the magnetofluid; the geometry of the toroid, which is allowed to have variable cross sections; the safety factor of the magnetofluid; and perhaps others.

The pseudo-spectral algorithm combined with the volume penalization method are used. This numerical method was presented in detail in the Chap. II. In this chapter a toroidal geometry is carved out within the periodic volume.

In Sec. 2, we fix the geometry of the confined magnetofluid and write down the system of equations and boundary conditions that will govern the dynamics. An external forcing of the magnetic field provides the toroidal electric field which initiates and drives the current. A vacuum toroidal dc magnetic field, regarded as externally supported from outside the system, is also assumed to be present. In addition another toroidal component of the magnetic field is allowed to develop in time if the dynamics so dictate.

In Sec. 3, the results are presented. They are divided in four different parts. The first discusses the generation of toroidal velocities for a dissipative system. The second exposes the results where the nonlinear term is dominant and a comparison is made between different toroidal geometries. In the third section we study the effect of the variation of the safety factor and in the last part, how the system evolves if the imposed toroidal magnetic field is inverted. We illustrate in detail the development of the driven magnetofluid configurations and the development of spontaneous toroidal rotation.

GEOMETRY AND GOVERNING EQUATIONS

In the MHD approximation the plasma is described as a charge-neutral conducting fluid. Despite its low complexity compared to kinetic descriptions it can give rise to a wealth of intricate phenomena and its analytical treatment is only possible in some simplified cases, either in the absence of velocity fields [START_REF] Berkowitz | Proceedings of the Second United Nations International Conference on the Peaceful Uses of Atomic Energy[END_REF][START_REF] Shafranov | Plasma equilibrium in a magnetic field[END_REF] or in the absence of non-linear interactions [START_REF] Bates | Toroidal visco-resistive magnetohydrodynamic steady states contain vortices[END_REF]. If one considers the complete problem one necessarily needs to consider a discretized numerical approximation of the full nonlinear system. Here we recall the equations we consider for the reader's convenience: the dimensionless incompressible viscoresistive MHD equations for the velocity field u and for the magnetic field B, in 'Alfvénic' units [START_REF] Kamp | Toroidal steady states in visco-resistive magnetohydrodynamics[END_REF].

∂ u ∂t -M -1 ∇ 2 u = -∇ P + 1 2 u 2 + u × ω + j × B, (3.1)
∂ B ∂t = -∇ × E, (3.2)
E = S -1 j -u × B, (3.3) ∇ • u = 0, ∇ • B = 0, (3.4)
with the current density j = ∇ × B, the vorticity ω = ∇ × u, the pressure P and the electric field E. These equations are non-dimensionalized using the toroidal Alfvén speed C A = B re f / √ ρ µ 0 as typical velocity, with B re f = 1.2 the reference toroidal magnetic field at the center of the torus (R = R 0 = 0.55π ≈ 1.73 and Z = 0 for both considered geometries), ρ the density and µ 0 the magnetic constant. We will exclusively consider two toroidal geometries with differently shaped cross-sections (see Fig. 3.1). The reference length L is the diameter of the cross section for the circular case and is the minor diameter for the asymmetric 'D' shape (L = 0.6π ≈ 1.88 for both geometries). The 'D' shape parametric equation is a modified version of the formula given by Manickam [START_REF] Manickam | Stability of n=1 internal modes in tokamaks[END_REF],

R(t) = R 0 + L 2 [cos(t -α + δ sin(t)) cos(ζ ) -κ sin(t) sin(ζ )] , (3.5
)

Z(t) = L 2 [cos(t -α + δ sin(t)) sin(ζ ) + κ sin(t) cos(ζ )] , (3.6)
with t ∈ [0, 2π], δ the triangularity, κ the ellipticity, α the asymmetry and ζ the rotation angle. For the considered asymmetric cross section the following values of these parameters are chosen: δ = 0.5, κ = 2.1, α = 0.4 and ζ = 0.15.

The MHD equations are completed by the initial and boundary conditions of the problem, and two dimensionless quantities: the viscous Lundquist number (M) and the Lundquist number (S) defined as

M = C A L ν , S = C A L λ , (3.7)
with λ the magnetic diffusivity and ν the kinematic viscosity. The ratio of these two quantities is the magnetic Prandtl number Pr = ν/λ , which we have chosen unity in the present study, thereby reducing the number of free parameters, which characterize the magnetofluid, to one, the viscous Lundquist number, M. Previous investigations indicate that it is the geometric mean of the viscosity and the magnetic diffusivity which is important to the dynamics [START_REF] Cappello | Bifurcation in viscoresistive MHD: The Hartmann number and the reversed field pinch[END_REF][START_REF] Shan | On the role of the Hartmann number in magnetohydrodynamic activity[END_REF]. In setting the Prandtl number to one, a change in the viscous Lundquist numbers, M or S, is equivalent to a change in the Hartmann number. In the ideal MHD framework a scalar-pressure equilibrium state is assumed in which

u = 0, j × B = ∇P. (3.8)
This equilibrium is possible in a cylindrical geometry, for instance in z-and θ -pinches. It is shown in [START_REF] Montgomery | Toroidal resistive MHD equilibria[END_REF][START_REF] Bates | A toroidal boundary-value problem in resistive magnetohydrodynamics[END_REF] that in the case of finite conductivity such an equilibrium is not possible in a toroidal geometry if irrotational toroidal magnetic and electric fields are applied. A steady state in Faraday's law imposes the toroidal electric field to be irrotational in the region of interest. The chosen spatial dependence for E 0 T is ∝ 1/R. In the simple case of a space-uniform conductivity, which we consider in the present study, the current density has the same dependence. The form for the imposed toroidal magnetic field, which is also proportional to 1/R, comes from the integration of Ampère's law on a toroidal loop. So the externally imposed magnetic field and toroidal, laminar, voltage-driven current density are given by,

B 0 T (R) ∝ R 0 R e T , j 0 T (R) ∝ R 0 R e T . (3.9)
The toroidal magnetic and current density profiles give the imposed three-dimensional helical magnetic field

B 0 = B 0 T + B 0 pol , with B 0 pol = B 0 R e R + B 0 Z e Z .
The poloidal magnetic field is calculated from the current density distribution j 0 T (R). The reference toroidal current density at the center of the geometry (at R = R 0 and Z = 0) is j T re f = 0.5. For the details of generating the poloidal magnetic field in general geometries numerically we refer to Appendix D. Here e T , e R and e Z are unit vectors in the toroidal/azimuthal, radial and vertical directions respectively (Fig. 3.1).

The toroidal magnetic field magnitude is tuned to have an edge safety factor q = rB 0 T | wall /R 0 B 0 P | wall = 5.7 for the asymmetric geometry and q = 3.3 for the symmetric cross section. A bar over a symbol indicates an average over the entire boundary. These safety factor values will be used for the majority of studied cases. The pinch-ratio associated to these values of q, defined as the ratio between the wall-averaged poloidal and the volume-averaged toroidal imposed magnetic field, Θ 0 = B 0 P / B 0 T = 0.16, is the same for both geometries. The resulting three-dimensional magnetic field lines are visualized for both geometries in Fig. 3.2. The Lorentz force resulting from the calculated poloidal field B 0 pol and the imposed toroidal current density j 0 T is not curl-free [START_REF] Montgomery | Toroidal resistive MHD equilibria[END_REF][START_REF] Bates | A toroidal boundary-value problem in resistive magnetohydrodynamics[END_REF]. Since the curl of a pressure gradient is necessarily zero, the equilibrium described by (3.8) becomes impossible and additional terms of Eq. (3.1) need to be taken into account to balance the equation. Since all other terms in (3.1) are proportional to (or quadratic in) the velocity, the resulting state must be dynamic. In other words if we take the curl of Eq. (3.1) we end with the vorticity equation,

∂ ω ∂t -M -1 ∇ 2 ω -∇ × (u × ω) = ∇ × (j × B) = 0, (3.10)
we observe that if the Lorentz force term is not curl-free, it acts as a source of vorticity: a toroidal plasma, described by viscoresistive MHD, confined by curl-free toroidal electric and magnetic fields, necessarily moves! It is true that the rationale described above depends on the choice of the electric conductivity, which was assumed to be uniform. It was however shown [START_REF] Montgomery | Resistive magnetohydrodynamic equilibria in a torus[END_REF] that to satisfy Eq. (3.8) in a torus, very unusual profiles of the electrical conductivity must be assumed. The simple case of constant magnetic resistivity is then treated in this study. The case of non-uniform resistivity profiles is one important perspective.

It follows from the foregoing that it is necessary to take into account all other terms in the MHD equations, and analytical treatment becomes impossible unless symmetries are assumed. To study the full dynamics we are obliged to solve numerically the system and this is what is done in the present investigation. Equations (3.1)-(3.4) are discretized with a Fourier pseudo-spectral method on a Cartesian grid. To impose the boundary conditions we use the volume-penalization technique. The method is presented in detail for threedimensional viscoresistive MHD equations in Chap. II.

The total magnetic field is decomposed into a base component and a perturbation,

B = B 0 + B ′ . (3.11)
Numerically only the perturbation of the magnetic field is computed, the base magnetic field, B 0 , computed from (3.9) is fixed and it is introduced in the Navier-Stokes equation and in the induction equation as follows,

∂ u ∂t -M -1 ∇ 2 u = -∇ P + 1 2 u 2 + u × ω + j ′ + j 0 × B ′ + B 0 (3.12) ∂ B ′ ∂t -S -1 ∇ 2 B ′ = ∇ × u × B ′ + B 0 (3.13)
To close the equations we have the incompressibility of the velocity field and the solenoidal constraint on the perturbed part of the magnetic field,

∇ • u = 0, ∇ • B = 0. (3.14)
The boundary conditions are to be no-slip, u| wall = 0, for the velocity. For the magnetic perturbation, the poloidal component and the component normal to the wall vanish, B ′ P wall = B ′ ⊥ wall = 0, while the toroidal component is free. The normal component B ⊥ vanishing at the wall physically corresponds to perfectly conducting boundary conditions. The zero poloidal fluctuations B ′ P wall are imposed for numerical convenience. Since the perturbed magnetic field remains small compared to the field B 0 in the present investigation, we do not think that this simplification significantly influences the results.

The initial condition for the simulations is zero magnetic perturbations and zero velocity. The simulations are carried out on a cubic domain of size (2π) 3 for the asymmetric and (2π × 2π × π) for the symmetric cross section consisting of 256 3 grid points. We fix the penalization parameter to η = 5 • 10 -4 . The time step is adaptive and the chosen CFL coefficient is 0.1.

RESULTS AND DISCUSSION

The results are divided into four different parts. The first shows the solution of the simulations at a low viscous Lundquist number, where an illustration of the generation of toroidal velocities is presented. The second exposes the calculations at higher viscous Lundquist, where the flow behavior of the plasma changes towards a dominantly toroidal flow. In the third section we compare, at fixed transport coefficients, simulations carried out for different safety factors and in the fourth section we show the results when the toroidal magnetic field is reversed.

GENERATION OF TOROIDAL VELOCITIES AT LOW VISCOUS LUNDQUIST NUMBER

In this section the calculations are performed for a low viscous Lundquist number, M = 23, in the geometry with symmetric cross section and q = 3.3. All the results are presented when the system has reached a statistically stationary state. Fig. 3.3 shows the presence of a poloidal flow, a pair of counterrotating vortices in the poloidal plane. In this case the flow topology is almost axisymmetric with respect to the Z-axis. To visualize more clearly the toroidal velocities and the double poloidal recirculation, the azimuthally averaged velocity field is presented in Fig. 3.4. We distinguish four different zones, where the toroidal velocity changes sign, and the already mentioned "double smoke ring". Indeed, in the limit of vanishing nonlinearity, Bates and Montgomery [START_REF] Bates | Toroidal visco-resistive magnetohydrodynamic steady states contain vortices[END_REF] showed analytically that the steady state solution is a pair of poloidally rotating vortices, aligned with the toroidal direction. The origin of toroidal velocities was demonstrated for vanishing viscous Lundquist in a rectangular cross section [START_REF] Kamp | Toroidal flows in resistive magnetohydrodynamic steady states[END_REF]. For a circular cross section and at low M number we will illustrate the generation of this velocity component. First, we illustrate that the forcing appearing in the vorticity equation (3.10) creates a toroidal vorticity with opposite sign in relation to the mid-plane of the torus (see Fig. 3.5 (a)). This creates automatically a radial velocity that will interact with the imposed toroidal magnetic field (Fig. 3.5 (b)). The interaction will produce a perturbation to the toroidal magnetic field (B ′ T ). Notice that this magnetic field will have positive and negative areas located in a similar position as the radial velocity (Figs. [START_REF] Caltagirone | Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide-milieux poreux: application à la convection naturelle[END_REF].5 (b) and (c)). It was shown [START_REF] Kamp | Toroidal flows in resistive magnetohydrodynamic steady states[END_REF] that the equation giving the first order perturbed toroidal magnetic component

B ′(1)
T is,

∇ 2 (B ′(1) T e T ) ∼ -u R B T re f R 2 e T .
(3.15)

The sign of the right hand side, will only depend on the sign of u R and of the imposed toroidal field B T re f . It follows that the curl of the perturbed toroidal magnetic field (B ′ T ) will produce a poloidal current density, j' pol = ∇ × B ′ T (Fig. 3.6 (a)). The imposed poloidal magnetic field B 0 pol will then interact with the perturbed current density j' pol to create a toroidal Lorentz force (see Figs. [START_REF] Caltagirone | Sur les conditions hydrodynamiques au voisinage d'une interface milieu fluide-milieux poreux: application à la convection naturelle[END_REF].6 (b) and (c)). The Lorentz force will finally induce toroidal velocities. Note that there is a similarity in the negative and positive zones between the toroidal velocity and the toroidal Lorentz force fields (see Figs. 3.4 and 3.6 (c)). We note that the sign in the toroidal Lorentz force depends exclusively on the angle between j' pol and B 0 pol . As a consequence this angle influences directly the toroidal velocity direction. Another way, to apprehend the fact that the poloidal flow interacts first with the magnetic field creating subsequently toroidal velocities, is to see the time evolution of the different velocity components. The velocities in the poloidal plane (in the poloidal direction P and in the minor radius direction r) grow first. After that the toroidal velocity is generated (see Fig. 3

ω T -5 • 10 -2 -3 • 10 -2 0 3 • 10 -2 5 • 10 -2 (a) u R (b) B ′ T (c)

.7).

At low viscous Lundquist number the dominant velocities are in the poloidal plane and form two counterrotating vortices. Small toroidal velocities appear and they form a quadrupole with alternating positive and negative directions. The analytical results published by Bates and Montgomery [START_REF] Bates | Toroidal visco-resistive magnetohydrodynamic steady states contain vortices[END_REF] are in good agreement. Also the numerical generation of toroidal velocities agrees with the calculations made by Kamp et al. [START_REF] Kamp | Toroidal flows in resistive magnetohydrodynamic steady states[END_REF].

SIMULATIONS FOR HIGHER VISCOUS LUNDQUIST NUMBERS

In this section the calculations are made for the asymmetric cross section with fixed q = 5.7 and for the circular cross section, q = 3.3. The viscous Lundquist numbers are modified changing the transport coefficients ν and λ (with Pr = 1), keeping the geometry and the reference toroidal magnetic field unchanged, B re f = 1.2.

With higher viscous Lundquist numbers it takes longer for the system to reach the saturated state. In the first instants an oscillatory behavior is present (see for example the different energy evolutions in Figs. 3.8 and 3.9). The kinetic and the fluctuating magnetic energy oscillate in opposition of phase, but these oscillations are damped out in a finite time. In the following section we will analyse and compare the different simulations when the system has reached this non-oscillatory steady state. The flow evolution is quantified in Fig. 3.12, where we observe that the principal direction of the flow is toroidal if M is raised beyond ∼ 40. The square toroidal velocity saturates for increasing M at a value of ∼ 80% of the total square velocity for the asymmetric cross section and at ∼ 60% for the circular profile. This toroidal organization of the flow is consistent with the tendency of the velocity field to align with the magnetic field, as is illustrated in Fig. 3.13, where we compute the average (over the toroidal domain) of the absolute value of the cosine of the angle between the velocity and magnetic field. This quantity is equal to one if the velocity and the magnetic field are perfectly aligned or antialigned. The evolution of the ratio u 2 T / u 2 with M shows the same trend as the alignment between the magnetic and the velocity field.

An important difference is observed between the flows that are generated in the two geometries. The volume averaged toroidal angular momentum is defined by

L T = 1 V V Ru T dV. (3.16)
For the torus with circular cross section, this quantity is zero to a good computational approximation (< 10 -15 ). The up-down anti-symmetry of the velocity field is responsible for this absence of toroidal angular momentum. However, for the torus with asymmetric cross section this is not the case. There is a symmetry breaking of the flow and the volume integral of the toroidal velocity is nonzero. In our calculations this can be visualized in the azimuthally averaged velocity fields in Fig. 3.14. It is more clear for the last case, at M = 4524, that the positive toroidal velocity occupies a larger part of the poloidal plane than the negative toroidal velocity. To quantify the amount of dissymmetry in the flow we present the evolution of the normalized toroidal angular momentum with M (see Fig. 3.15). This quantity increases with the viscous Lundquist number. This up-down symmetry effect is in agreement with time-independent computations [START_REF] Kamp | Toroidal steady states in visco-resistive magnetohydrodynamics[END_REF] and also with gyrokinetic simulations and experiments [START_REF] Camenen | Transport of parallel momentum induced by current-symmetry breaking in toroidal plasmas[END_REF][START_REF] Camenen | Experimental demonstration of an up-down asymmetry effect on intrinsic rotation in the TCV tokamak[END_REF].

Furthermore in Fig. 3.14 we can observe the two counterrotating vortices. They are still present at higher viscous Lundquist but undergo a deformation and their center is shifted outwards. The larger toroidal velocities concentrate near the boundaries as well as the poloidal speeds (this can be seen from the stream function isocontours that tend to converge near the boundaries). Nevertheless the velocity magnitude is globally less important for high M. In fact the kinetic energy has a maximum and then decreases if the viscous Lundquist number is raised (see Fig. 3.16). This behavior is explained by the decrease of the magnitude of the Lorentz force with the viscous Lundquist number in the center of the domain. Indeed, the plasma seems to self-organize to a state with a forcefree region in the center, a behavior also observed for straight-cylinder computations at high pinch ratio [START_REF] Dahlburg | Driven, steady-state RFP computations[END_REF].

The evolution with M of the root mean square (RMS) value of this force is presented in Fig. 3.17. The spatial distribution of the norm of the Lorentz force vector in the poloidal plane is visualized for the asymmetric geometry in Fig. 3.18. The vanishing of the Lorentz force in the core comes from the alignment between the magnetic and current density fields. A measure giving the alignment between these three-dimensional quantities is the volume-averaged current helicity defined as

H j = j • B j B . (3.17)
We observe (Fig. 3.19) that for increasing viscous Lundquist number the global current density and magnetic field tend to be oriented in the same direction, the quantity in the figure approaches the unit value. This causes the Lorentz force term to decrease for higher M in the center of the domain, the magnitude of the imposed toroidal current density and magnetic fields remaining constant. Mainly, the variation of the alignment between j and B occurs in the poloidal plane. To quantify the alignment among the poloidal current density and the poloidal magnetic field we compute the volume-averaged absolute value of the cosine of the angle between these two fields,

|cos Φ| = j pol • B pol j pol B pol , (3.18)
where J pol and B pol are the projections of J and B on the poloidal plane. This quantity at low viscous Lundquist is smaller compared to the value of the current helicity at the same M number (see Figs. 3.19 and 3.20). With increasing viscous Lundquist the cosine of this angle grows and approaches unity. There is a stronger change in the alignment between the current density and magnetic field in the poloidal plane. This poloidal alignment makes the toroidal Lorentz force vanish in the core of the domain. Whether or not the Lorentz force term reaches an asymptote at higher M or if a transition to another state exists remains an open question.

The system is almost axisymmetric around Z but small fluctuations around the toroidally averaged fields exist, defined as

ũ = u -u T , B′ = B ′ -B ′ T . (3.19)
The most important normalized fluctuations around the axisymmetric state are in the velocity field (Fig. 3.21), they are localized at the boundaries (see Fig. Lundquist number are presented respectively in Figs. 3.21 and 3.22. For the highest viscous Lundquist, M = 4524 and asymmetric cross section, we have the maximum ratio ũrms /u rms ∼ 0.14. The greatest normalized departure from axisymmetry for the perturbed magnetic field is also at M = 4524 for the 'D' cross section, B′ rms /B ′ rms ∼ 0.015. In fact for the magnetic field, the fluctuations are of the same order of magnitude as the velocity field, but the magnitude of the perturbed magnetic field is larger, hence the normalized quantities are smaller. The distribution of the perturbations in the two-dimensional plane (Figs.

E k τ q = 5.7 q = 4.8 q = 3.8 q = 2.9 1 • 10 -6 1 • 10 -5 1 • 10 -4 1 • 10 -3 0 4 8 12 16 20
E k τ q = 5.7 q = 4.8 q = 3.8 q = 2.9

INFLUENCE OF THE SAFETY FACTOR ON THE DYNAMICS

The study of the influence of the safety factor q is presented in this section where we consider only the asymmetric cross section geometry and the transport coefficients are kept constant (ν = λ = 2 • 10 -3 ). We recall that for all the simulations presented in this manuscript the magnetic Prandtl number is equal to one, Pr = 1. In this case the viscous Lundquist number varies because the reference magnetic field used for its calculation is the imposed toroidal component and to modify the safety factor the magnitude of this field is changed, as also done in experiments [START_REF] Rice | Observations of core toroidal rotation reversals in Alcator C-mod ohmic L-mode plasmas[END_REF][START_REF] Rice | Ohmic energy confinement saturation and core toroidal rotation reversal in Alcator C-mod plasmas[END_REF]. The parameter q takes four different values. We recall that the safety factor is defined as the ratio between the wall-averaged toroidal and poloidal imposed magnetic fields,

q = rB 0 T | wall R 0 B 0 P | wall . (3.20)
The values of the viscous Lundquist number associated to each safety factor are presented in Tab. The evolution of the total kinetic energy and the magnetic energy of the perturbation is similar for all the studied cases (see Figs. 3.25 and 3.26). The main difference is the

1 • 10 -6 1 • 10 -5 1 • 10 -4 1 • 10 -3 1 • 10 -2 0 20 40 60 80
100 120

E mag τ q = 5.7 q = 4.8 q = 3.8 q = 2.9 1 • 10 -6 1 • 10 -5 1 • 10 -4 1 • 10 -3 1 • 10 -2 0 4 8 12 16 20
E mag τ q = 5.7 q = 4.8 q = 3.8 q = 2.9 magnitude of the energies that is higher if the safety factor is small. At the steady state the dependence of the kinetic energy on the safety factor is visualized in Fig. 3.27.

The growth of the kinetic energy with decreasing q is in agreement with the reduction of the current helicity value (Fig. 3.28). Hence the Lorentz force term is stronger for a low safety factor. It is also observed that in the toroidal direction the Lorentz force increases, since the alignment between the poloidal current density and the poloidal magnetic field is less important for small q (Fig. 3.29). This variation is smaller compared to the variation caused by the transport coefficients modification, as shown in Sec. 3.2.

As in the previous section the toroidal velocity dominates, but the ratio u 2 T / u 2 decreases with decreasing q (Fig. 3.30). Also the alignment between the magnetic and velocity field is less important (inset Fig. 3

.30).

An important feature is the change of sign in the volume averaged toroidal angular momentum, found also in experimental observations [START_REF] Rice | Observations of core toroidal rotation reversals in Alcator C-mod ohmic L-mode plasmas[END_REF][START_REF] Rice | Ohmic energy confinement saturation and core toroidal rotation reversal in Alcator C-mod plasmas[END_REF][START_REF] Scarabosio | Toroidal plasma rotation in the TCV tokamak[END_REF], when the toroidal magnetic field, hence the safety factor, is varied (Fig. 3.31). In our case the averaged angular momentum changes completely in sign, it passes from negative to positive for increasing q. The two-dimensional azimuthally averaged toroidal velocities (Fig. 3.32) show the increase of the area in which the toroidal velocity is negative when the safety factor is decreased. For the lowest value of q that we consider, the vertical cut (Fig. 3.32 (d)) shows larger velocities and a small downward shift of the position where the toroidal velocity changes sign. This displacement enlarges the negative velocity area. The growth of the negative toroidal velocity is better visualized in the cuts along the direction of the big radius (Fig. 3.33). For decreasing q the velocities tend to be more peaked and near the center of the torus a region appears where the toroidal velocity is negative. We notice that the change of sign of the toroidal velocity mainly occurs in the center of the geometry. Close to the boundaries the toroidal component grows but does not reverse sign.

As presented in the previous section small fluctuations around the azimuthal average exist. We see in Fig. 3.34 that the magnitude of these fluctuations is relatively insensitive to the value of the safety factor. The change is just of a few percent for the normalized T / u 2 as a function of q. In the inset we show the average over the domain of the absolute value of the cosine of the angle between the velocity field and magnetic field. velocity fluctuations. It is larger for the normalized magnetic fluctuations, but it remains below ∼ 20% (Fig. 3.35). Hence the safety factor variation, in the considered range, does not increase substantially the non-axisymmetric perturbations.

INFLUENCE OF THE REVERSAL OF THE IMPOSED TOROIDAL MAGNETIC FIELD

The simulation with inverted toroidal magnetic field is performed for q = 5.7 and M = 1131. The results show that the velocity reverses sign in the whole domain (Fig. 3.36).

The counterrotating poloidal vortices are unchanged, only the toroidal velocities are affected. In Fig. 3.36 (c) the profiles are exactly symmetric with respect to the vertical axis.

Basically, what happens is that the perturbed toroidal magnetic field reverses its sign and this generates an inverse poloidal current density. The existing poloidal magnetic field associated with the inverted poloidal current density field gives an opposite toroidal Lorentz (c) -2.5 • 10 -3 0 2.5 • 10 -3 uT q = 5.7 q = 4.8 q = 3.8 q = 2.9

(d) u T R q = 5.7 q = 4.8 q = 3.8 q = 2.9 force. Finally, this Lorentz force will make the toroidal velocities reverse in all the domain. We can write the three components of the Lorentz force in cylindrical coordinates:

     F R = j T B Z -1 R ∂ (RB T ) ∂ R B T , F T = ∂ B T ∂ Z B Z + 1 R ∂ (RB T ) ∂ R B R , F Z = -∂ B T ∂ Z B T + j T B R . (3.21)
The inversion of the sign of B T transforms the original Lorentz force vector

(F R , F T , F Z ) into (F R , -F T , F Z ).
Only the toroidal component is affected. Hence the poloidal velocities are unchanged but the toroidal velocities are inverted.

CONCLUDING REMARKS ABOUT THE GENERATION OF VELOCITIES IN TOROIDAL GEOMETRIES

In this chapter it was demonstrated numerically, solving the fully nonlinear time-dependent resistive MHD equations, that in a toroidal geometry, assuming constant transport coefficients, if the imposed toroidal magnetic and toroidal electric fields are irrotational, the conducting flow inside a torus necessarily moves. The reason for this is that the curl of the Lorentz force resulting of the imposed fields is nonzero. It follows that the gradient of a scalar (in this case the pressure) can not balance the equation. Consequently vorticity is created. This vorticity in the toroidal direction creates poloidal velocities. The poloidal velocities interact with the imposed toroidal magnetic field creating a perturbation that gives rise to a poloidal current density. This current density associated with the existing poloidal magnetic field produces a toroidal Lorentz force. As a consequence toroidal velocities appear. The angle between the poloidal current density and poloidal magnetic field plays an important role in the determination of the toroidal velocity direction. For a low viscous Lundquist number the system tends to produce small toroidal velocities, the dominant flow being a pair of counterrotating vortices in the poloidal plane. A dramatic change occurs when the viscous Lundquist number is increased. There is a transition from a dominantly poloidal to dominantly toroidal flow. This transition is in agreement with the tendency of the velocity field to align with the magnetic field.

-2 • 10 -3 0 2 • 10 -3 u T B 0 T ⊕ B 0 T ⊖ (c)
Two different toroidal geometries are considered in the present study, one with an updown symmetric and the other with an asymmetric cross section. A fundamental difference exists between both studied cases: the volume-averaged angular momentum is zero for the symmetric case, while for the asymmetric cross section a finite volume-averaged angular momentum appears. There is a breaking in the up-down symmetry of the flow and a toroidal preferred direction emerges. This volume-averaged normalized angular momentum tends to increase with the viscous Lundquist number.

Nevertheless the kinetic energy decreases with increasing nonlinearity, since the total magnetic and current density fields tend to align in the center of the domain. The limitation in the numerical resources prevents the study of this system for larger viscous Lundquist numbers. It remains an open question if there will be a continuous increase of the alignment between the magnetic and current density field or if a transition exists.

When the safety factor is decreased while maintaining the transport coefficients constant, the kinetic and fluctuating magnetic energy become higher. The main qualitative effect is the influence on the toroidal velocity direction. There is a change in the volumeaveraged angular momentum that reverses sign. For low q it is negative and at large safety factor it becomes positive. Mainly in the center of the domain, for decreasing q, the region in which the velocity is negative becomes larger, at expense of the region with positive toroidal velocity. Near the boundaries the toroidal velocity direction remains unchanged.

The last part of the chapter was dedicated to the influence of the reversal of the toroidal magnetic field. It is shown that it plays a role only in the toroidal velocities. The reversal changes the sign of the poloidal current density, that gives rise to the toroidal Lorentz force. In consequence the toroidal force reverses in the whole volume making the toroidal velocities reverse their direction compared to the original case.

IV MAGNETOHYDRODYNAMICS FOR HIGH PINCH RATIOS IN TOROIDAL AND CYLINDRICAL GEOMETRIES

INTRODUCTION

The Reversed Field Pinch (RFP) is characterized by pinch ratios bigger than unity. This device, together with the tokamak and the stellarator, is one of the magnetic fusion configurations that could achieve thermonuclear fusion. A practical advantage of the RFP as compared to the tokamak geometry is that to confine the plasma the RFP needs a smaller externally imposed magnetic field. The toroidal field is enhanced through the self-organization of the plasma in the center, via a dynamo mechanism induced by the plasma itself, and decreases at the plasma edge. The reduction of the toroidal field can be very large and the value at the edge can reverse sign.

The dynamo mechanism can be explained if we consider the following simple example which we adopted from Escande et al. [START_REF] Escande | Single helicity: a new paradigm for the reversed field pinch[END_REF]. Let us consider a wire placed in the interior of a solenoid (Fig. 4.1). The current in the wire and the magnetic field created by the solenoid are in the same direction. If this wire is exactly at the geometrical center, it is in unstable equilibrium, because the currents in the solenoid attracting the wire radially compensate one another. If this wire is deformed by a small perturbation, it will bend, forming a helical geometry that will follow the original solenoid. As a consequence there is a second 'solenoid' which is created in the interior of the first one. This second winding will produce a magnetic field in the same direction as the original in the center, but it will induce an opposite magnetic field on its exterior. In the center it will therefore enhance the magnetic field but near the edges the magnetic field is decreased and can even reverse sign.

As a consequence of the high pinch ratio, magnetohydrodynamic instabilities such as the above discussed kink-instability appear in the RFP, which increase the turbulence and lower the energy confinement performance. Generally this chaotic regime is called multiple-helicity (MH) state because a multitude of helically extended modes interact. However in the last decades quasi-single helicity (QSH) states were observed in experiments where the full turbulent regime is avoided and one helical mode predominates above the others [START_REF] Escande | Quasi-single-helicity reversed-field-pinch plasmas[END_REF][START_REF] Martin | Quasi-single helicity states in the reversed field pinch: Beyond the standard paradigm[END_REF][START_REF] Nordlund | The structure and stability of the reversed field pinch magnetic equilibrium in Extrap T1[END_REF][START_REF] Brunsell | Coherent magnetic field fluctuations and locked modes in a reversed-field pinch[END_REF][START_REF] Sarff | Increased confinement and β by inductive poloidal current drive in the reversed field pinch[END_REF]. In the QSH state there is a decrease of magnetic chaos and the formation of a coherent helical structure within the plasma. The reduction of particle transport has been observed in QSH states with respect to MH states [START_REF] Frassinetti | Experiments and modelling of active quasi-single helicity regime generation in a reversed field pinch[END_REF][START_REF] Predebon | Particle-transport analysis in reversed field pinch helical states[END_REF][START_REF] Gobbin | Numerical studies of transport mechanisms in RFX-mod low magnetic chaos regimes[END_REF]. Also it has been found that increasing the toroidal current makes the QSH regime more persistent [START_REF] Martin | A new paradigm for RFP magnetic self-organization: results and challenges[END_REF][START_REF] Lorenzini | Single-helical-axis states in reversed-fieldpinch plasmas[END_REF][START_REF] Lorenzini | Selforganized helical equilibria as a new paradigm for ohmically heated fusion plasmas[END_REF]. The QSH state is responsible for the creation of internal transport barriers (ITB) that improve the confinement time by a factor of two [START_REF] Carraro | Improved confinement with internal electron transport barriers in RFX-mod[END_REF][START_REF] Piovesan | Magnetic order and confinement improvement in high-current regimes of RFX-mod with MHD feedback control[END_REF]. It has also been discovered relatively recently that external transport barriers (ETB) are created for low plasma densities that increase in 30% the confinement time. The origin of this barrier is not well understood but it could be related to a shear flow [START_REF] Puiatti | Internal and external electron transport barriers in the RFX-mod reversed field pinch[END_REF].

The majority of numerical studies of the RFP, using the MHD approximation, consider periodic cylindrical geometries, ignoring the influence of the curvature of the magnetic field. In the eighties the first MHD calculations were performed using the nonviscous MHD model (see e.g. [START_REF] Dahlburg | Driven, steady-state RFP computations[END_REF][START_REF] Aydemir | Compressibility as a feature of field reversal maintenance in the reversed-field pinch[END_REF]) and in the nineties viscoresistive calculations were performed [START_REF] Shan | Nonlinear magnetohydrodynamics by Galerkin-method computation[END_REF][START_REF] Shan | On the role of the Hartmann number in magnetohydrodynamic activity[END_REF][START_REF] Shan | Global searches of Hartmann-numberdependent stability boundaries[END_REF]. In the last decades the viscoresistive MHD model has been taken into account using a constant-pressure constant-density approximation (SpeCyl code) [START_REF] Cappello | Bifurcation in viscoresistive MHD: The Hartmann number and the reversed field pinch[END_REF][START_REF] Cappello | Bifurcation in the MHD behaviour of a self-organizing system: the reversed field pinch (RFP)[END_REF][START_REF] Bonfiglio | Dominant electrostatic nature of the reversed field pinch dynamo[END_REF][START_REF] Veranda | Impact of helical boundary conditions on nonlinear 3D magnetohydrodynamic simulations of reversed-field pinch[END_REF]. Also in a cylinder, but using the full compressible viscoresistive MHD equations, the effect of the compressibility has been studied for a constant uniform [START_REF] Onofri | Compressibility effects in the dynamics of the reversed-field pinch[END_REF][START_REF] Onofri | Effects of compressibility and heating in magnetohydrodynamics simulations of a reversed field pinch[END_REF] and nonuniform resistivity [START_REF] Onofri | Single-helicity states in compressible magnetohydrodynamics simulations of the reversed-field pinch with nonuniform resistivity[END_REF].

Three dimensional codes solving the MHD equations exist, but there is not yet, at our knowledge, a complete study of the dynamics of a toroidal magnetofluid in the RFP context. A recent three dimensional code with a toroidal geometry implementation (PIXIE3D) has been compared with the cylindrical SpeCyl code [START_REF] Bonfiglio | Nonlinear three-dimensional verification of the SPECYL and PIXIE3D magnetohydrodynamics codes for fusion plasmas[END_REF][START_REF] Cappello | Equilibrium and transport for quasi-helical reversed field pinches[END_REF]. Also, a MHD toroidal code (MIPS) has been used to study the dynamics of the RFP RELAX [START_REF] Mizuguchi | Modeling of formation of helical structures in reversed-field pinch[END_REF].

In this chapter we study the helical plasma configurations using the viscoresistive MHD model in a torus and in a periodic cylinder. In general in real experiments the β P (average ratio of the plasma pressure over the poloidal magnetic pressure, 2 P /B 2 P (a)) is considerably below unity so it is believed that the toroidal outward shift of the magnetic surfaces is small such that equilibrium and stability are well described in a cylindrical approximation. Here we perform cylindrical and toroidal simulations to find out if there are important changes between these two geometries and if the cylindrical approximation can be justified.

First we discuss, in Sec. 2, the results of simulations in toroidal geometry. We investigate in some detail the transition from an axisymmetric state towards different helical regimes in a toroidal geometry. The pinch ratio of the imposed magnetic field will be varied and we will consider two different viscous Lundquist numbers. For the cylindrical case (Sec. 3) the self-organisation of the magnetofluid is investigated. As for the toroidal geometry, the pinch ratio is varied and two viscous Lundquist numbers are considered. In Sec. 4 a comparison is made between the toroidal and cylindrical results, where it is shown that even if the two configurations share some characteristics, a fundamental difference between the two flows exist, in particular when the torus and cylinder results are compared with experimental data (Sec. 4.2).

TOROIDAL SIMULATIONS

In this section we will first briefly illustrate the RFP dynamo effect. Then the transition from an axisymmetric to a helical state with increasing pinch ratio is discussed. The toroidal mode energies at the steady state are presented as well as the velocity field topology and different volume averaged quantities as a function of the pinch ratio. Finally the time evolution of the different toroidal modes for the velocity field are assessed and the diversity of dynamics in this system are underlined.

The toroidal geometry has a minor radius a = 0.3π and a big radius R 0 = 0.55π, the aspect ratio is Γ = R 0 /a ≈ 1.83 (Fig. 4.2). The size of the periodic box where the calculations are performed is (2π × 2π × π). We fix the penalization parameter to η = 5•10 -4 . The resolution for all the toroidal calculations is 128 3 are four, two of them are the pinch ratio and the reversal parameter,

Θ = B P (a) B T , F = B T (a) B T , (4.1)
with a the small radius. The remaining two are the viscous Lundquist number and the magnetic Prandtl number,

M = C A L ν , Pr = ν λ , (4.2)
with C A the poloidal Alfvén velocity, as reference length we take the diameter L = 0.6π, ν the viscosity and λ the magnetic diffusivity. We vary the viscous Lundquist number (M) with constant magnetic Prandtl number (Pr) and evaluate the system as a function of the pinch ratio at the final state (Θ). The imposed magnetic field for the toroidal geometry has the same topology as the one in the previous chapter. The imposed current in the center of the domain at R = R 0 (Fig. 4.2) is fixed at, J 0 = 1. The poloidal magnetic field is calculated with a precalculation as presented in Appendix D.

To vary the pinch parameter (Θ) numerically, we modify the toroidal magnetic field. Here we recall the profile of this quantity,

B 0 T ∝ R 0 R e T . (4.3)
The imposed poloidal component of the magnetic field is kept constant. For the toroidal calculations presented in this chapter the reference magnetic field (B re f ) is the surface averaged poloidal magnetic field at the boundaries, B re f = B P (a) = 0.35. This value of B re f is invariant for all the toroidal simulations.

DYNAMO EFFECT AND TRANSITION TO A HELICAL STATE

The toroidal simulations are started with zero velocity field and zero magnetic perturbations, and we fix the magnetic Prandtl number to Pr = 3. This value of the Prandtl number is chosen because Pr > 1 is expected to be found in fusion plasmas and in previous numerical studies Pr ∈ [1, 20] [18]. However, since the value of the viscosity in fusion plasmas is not so well defined, this choice should be considered of arbitrary nature and parametric studies on the influence of the Prandtl number constitute a very interesting perspective of the present investigation.

As shown in the previous chapter an axisymmetric flow develops, characterized by two counter-rotating vortices with their vorticity aligned with the toroidal direction. If the pinch ratio is increased an instability appears. This helical instability deforms the pre-existing axisymmetric state. This instability was absent in the previous chapter.

The instability consist in the growth of a helical-kink perturbation [START_REF] Escande | Single helicity: a new paradigm for the reversed field pinch[END_REF] with a toroidal mode number n. The value of this mode depends on the considered case: the variation of the diffusion coefficients and the pinch ratio modifies the value of the dominant toroidal mode. After a linear phase the instability saturates through nonlinear effects and different dynamics in the time evolution of the system exist. The mode, that grows at the early stages of the nonlinear phase, can be damped leaving its place to other modes. Oscillatory behavior appears for high values of the viscous Lundquist number and pinch ratios, as will be shown later.

F τ Θ = 1.6, M = 658 Θ = 1.5, M = 658 Θ = 2.0, M = 329
We now present the evolution of the volume averaged quantities for three different simulations. First we show the evolution of the pinch ratio and the reversal parameter, respectively, in Figs 4.3 and 4.4. From the Θ and F time evolution we can see the system reaching a steady state. Fluctuations in this state are present for the highest value of the viscous Lundquist number and pinch ratio. We note in Fig. 4.4 that at the steady state we do not have a complete reversal of the toroidal magnetic field at the edge, i.e., F remains positive. The reversal parameter has a transient period where it is negative but at the final times it grows and becomes positive for the considered cases. Notice that the Θ and F time evolution are similar. We can explain the evolution of these two quantities if we consider the dynamo effect that enhances the toroidal magnetic field. Fig. 4.5 shows the evolution of the perturbed toroidal magnetic field that is enhanced by the dynamo mechanism (phenomenon illustrated in Fig. 4.1). The increase of this selfgenerated magnetic field decreases the pinch ratio and at the same time increases the F parameter.

The dynamo mechanism will limit the maximum value that can be reached by Θ. Numerically, to rise the pinch parameter we decrease the imposed toroidal magnetic field (B 0 T ). For B 0 T tending to zero, the generated toroidal magnetic field tends to a constant value. The perturbed toroidal component generated by the dynamo mechanism dominates the system if B 0 T is sufficiently decreased. As a consequence the steady state pinch ratio reaches an asymptote when the initial pinch ratio tends to infinity † . The transition to a dominant helical system appears for a pinch parameter superior to one. This phenomenon is visualized in Fig. 4.6. Here the color of the streamlines, representing the toroidal velocity, is azimuthally invariant (Fig 4 .6 (left)). On the other hand, a helical deformation in the streamlines appears in Fig. 4.6 (right). This is a transition from a dominant axisymmetric state to a helical configuration. The streamlines form two helices, one with positive and the other with negative toroidal velocity. The energy contained in the different toroidal modes allows us to highlight this transition. In Fig. 4.7 (left) is presented the ratio of the energy of the 0, 3 and 4 toroidal modes over the total energy, for a viscous Lundquist number M = 329. It is shown that for a pinch ratio inferior to unity the energy is contained in the zero toroidal mode (the axisymmetric field). For a value of Θ ≈ 1.4 the zero mode is not longer dominant, the energy of the mode n = 3 contains ∼ 70% of the total energy. For increasing pinch ratio there are two other states. The first for Θ ≈ 1.8 where the dominant mode is not longer n = 3 but n = 4, with also ∼ 70% of the total energy. And for the highest pinch parameter (Θ ≈ 2) there is a regrowth of the n = 3 mode (∼ 70% of the total energy), the energy of the fourth mode being reduced.

For the highest viscous Lundquist number (M = 658), we have only considered Θ values larger than one (Fig. 4.7 (right)). As for the previous case we have the toroidal modes three and four that dominate alternatively. In this case for Θ ≈ 1.45, the dominant mode is n = 3 with ∼ 75% of the total energy. And for the highest values of the pinch ratio, Θ ≈ 1.57, we have the n = 4 mode dominating with ∼ 50% of the total energy. For this last case the energy is spread over a larger number of modes. Whether a completely turbulent state can be reached at higher values of the viscous Lundquist number remains an open question.

The helical modes, n = 3 and n = 4, dominating our system are visualized in Fig. 4.8 for different pinch parameter values but constant viscous Lundquist M = 658. In this figure a positive and a negative isocontour of the toroidal velocity are presented. Three and four lobes appear respectively, for Θ = 1.45 and Θ = 1.57.

Another visualization is presented in Fig. 4.9, where a perturbed magnetic isocontour in the core is colored by the toroidal velocity. For the first figure (left) two opposite toroidal velocities are clearly visualized as well as the toroidal mode deformation n = 3 in the toroidal magnetic isocontour. For increasing pinch ratio (Fig. 4.9 (right)) the mode n = 4 dominates the velocity field but with oscillations at the frontier between the positive and negative toroidal velocities. Now if we look at the shape of the magnetic isocontour, a correlation between the toroidal velocity pattern and the deformation of the toroidal magnetic component seems to exist. This isocontour takes a square shape (n = 4) that corresponds to the dominant mode in the velocity field. Experimentally (in the RFP MST) similar correlations have been measured between the fluctuation amplitude of the toroidal magnetic field and the velocity field projected in the direction of the experimental measure (chord projected) [START_REF] Piovesan | Measurements of the MHD dynamo in the quasi-single-helicity reversed-field pinch[END_REF].

Taking into account the spectra at the final state, Fig. 4.10 for Θ = 1.45 and Fig. 4.11 for Θ = 1.57, we see that there is clear correlation between the dominant toroidal modes of the velocity field and the toroidal magnetic field component. The spectrum maximum for the different two quantities is located at the same toroidal mode number, n = 3 for the low, and n = 4 for the high value of Θ. Notice that in Fig. 4.10 the mode n = 1 is important for the toroidal magnetic field, but is not large for the velocity field. It is not clear how the interaction between the two fields generates this n = 1 mode. It is probably a complex nonlinear effect composed of a dynamo action and a Lorentz force feedback.

In Fig. 4.12 the azimuthally averaged toroidal and poloidal components of the velocity field are visualized for three different values of the pinch parameter and a constant viscous Lundquist number, M = 329. Before the transition, Fig. 4.12 (a), the poloidal flow is composed of the already presented (see Chap. III) pair of counter-rotating vortices, and in the poloidal plane four different zones can be distinguished with positive and negative toroidal velocities. After the transition, Fig. 4.12 (b), the poloidal vortices are considerably enhanced and the toroidal velocity is distributed in just two regions with opposite toroidal directions. For the highest value of the pinch ratio, Fig. 4.12 (c), the 'double smoke ring' is still present and in the center of the poloidal plane a circular zone emerges with small toroidal velocity values. This is a consequence of the helical deformation of the magnetofluid. In the poloidal velocities we can notice large values close to the boundaries. There is a steep poloidal velocity gradient in this region. The double vortex flow pattern has been experimentally found in the RFX-mod device [START_REF] Bonomo | Flow measurements and modelling in helical RFX-mod equilibria[END_REF]. In the literature a flow is shown with poloidal mode number m = 1 with the same shear at the plasma edge (for the experiment at r/a ≈ 0.75). This steep gradient zone could be related to internal transport barriers (ITB) that are observed on a statistical basis in RFX-mod at r/a ≤ 0.7 [START_REF] Puiatti | Internal and external electron transport barriers in the RFX-mod reversed field pinch[END_REF]. For the toroidal velocity component, recent measurements have been performed in the MST experiment that show the toroidal velocity inversion for r/a ∼ 0.5 [START_REF] Ding | Kinetic stress and intrinsic flow in a toroidal plasma[END_REF].

-4 • 10 -3 -2 • 10 -3 0 2 • 10 -3 4 • 10 -3 (a) -4 • 10 -3 -2 • 10 -3 0 2 • 10 -3 4 • 10 -3 (b) -4 • 10 -3 -2 • 10 -3 0 2 • 10 -3 4 • 10 -3 (c)
The azimuthally averaged perturbed toroidal magnetic field is presented in Fig. 4.13 for three different pinch parameters as in the previous paragraph. In this figure we see, with increasing Θ, the enhancement of the toroidal field by the dynamo effect. For the axisymmetric simulation a very weak toroidal field develops (Fig. 4.13 (a)). For increasing pinch ratio the helical transition occurs and we observe a more important toroidal field (Fig. 4.13 (b)). In the last image (Fig. 4.13 (c)) the toroidal magnetic field continues to grow and becomes more peaked in the center.

EVOLUTION OF GLOBAL QUANTITIES

The kinetic energy and the magnetic energy of the perturbation are presented, respectively, in Figs. 4.14 and 4.15. The saturation, as well as the oscillations for high viscous Lundquist and pinch ratio values are visualized. At the steady state the kinetic energy of the system varies considerably with the pinch ratio (see Fig. 4.16). For increasing values of Θ the kinetic energy grows at the steady state. We observe that for low values of the pinch parameter the kinetic energy is more important if the viscous Lundquist number is low (as also observed in the previous chapter). But for a certain value of Θ (≈ 1.4) the high viscous Lundquist curve crosses the other curve and the system with lower viscosity has a larger kinetic energy. Also we notice that after the transition (between Θ = 0.65 and Θ = 1.44) the kinetic energy increases by a factor of ten (Fig. 4.16, curve M = 329). The magnetic energy of the fluctuation is calculated at the steady state, and its evolution with the pinch parameter is presented in Fig. 4.17. This energy increases with the pinch ratio. Its magnitude and its growth rate with Θ are large for the higher viscous Lundquist simulations.

0 1 • 10 -4 2 • 10 -4 3 • 10 -4 4 • 10 -4 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2.2 E k Θ M = 658 M = 329
For the case M = 329 we note a large increase of the perturbed magnetic energy between the first two Θ values.

The root mean square value of the Lorentz force is shown in Fig. 4.18. We observe that there is not a significant change in its value for increasing pinch ratio and for the considered viscous Lundquist numbers. There is a slight increase of the magnitude of the Lorentz force with the pinch parameter. This is compatible with the decrease of the alignment between the current density and magnetic field with Θ. This misalignment is stressed in Fig. 4.19 where the current helicity is presented. It is shown that for increasing pinch ratio the alignment between the magnetic and current density field is less important.

The spatial distribution of the norm of the Lorentz force (Fig. 4.20) shows that for increasing pinch ratio the central region of the torus becomes force free, as also observed in the previous chapter. But as stated above the root mean square of this quantity slightly grows, and this is because at the boundaries the magnitude of the Lorentz force strongly increases.

TIME EVOLUTION OF THE DIFFERENT TOROIDAL MODES

In the present study we have seen that the energy is distributed over a relatively small number of modes. In this section we will focus on the analysis of the system when it has arrived at the saturated state. Generally the simulations show a large variety of behaviors depending on the choice of the parameters. Small changes in the imposed parameters can generate significant modifications in the magnetofluid dynamics. With the evolution of the kinetic energy contained in the different toroidal modes, we want to illustrate the diversity in the dynamics and the sensibility of the magnetofluid to the variation of the pinch parameter and the viscous Lundquist number. The evolution of different mean square toroidal spectra is presented in Fig. 4.21 for different pinch ratios and fixed M = 329. Here the transition is clearly visible, it occurs between the cases Θ = 0.649 and Θ = 1.437. The first two figures (at the top) show no fluctuations and the zero toroidal mode is at least two orders of magnitude bigger than the second largest. For Θ = 1.437 the zero mode is not longer dominant. If we look at this figure the axisymmetric mode (n = 0) shows a non oscillatory behavior unlike the other modes that oscillate with a large time period. The dominant mode n = 4 alternates with others and at the final time a steady state can be suspected.

The dynamics of the magnetofluid for a higher pinch parameter (Θ = 1.820) are quite different. For the early times the system oscillates importantly and then for a time τ ≈ 700 (in poloidal Alfvén times) the system is completely steady, all the temporal oscillations are damped. A similar phenomenon appears for the two highest values of Θ presented in Fig. 4.21. The system oscillates strongly and then settles into a quiet steady state.

For increasing viscous Lundquist number (M = 658) a different behavior is present (see Fig. 4.22). After the transition to a helical state between Θ = 1.227 and Θ = 1.454 the system starts to strongly oscillate. Only for the case Θ = 1.454 we have a similar To quantify and compare the oscillations in the system we compute the standard deviation of the kinetic and magnetic energies in the saturated state (see Fig. 4.23). For both energies the standard deviation for the biggest viscous Lundquist number is one or two orders of magnitude bigger than for the low viscous Lundquist number. The figure also shows a transition with increasing pinch ratio. For a particular value of Θ there is a strong increase in the value of the standard deviation (at least one order of magnitude). This threshold is at the same position for the kinetic and magnetic energies. Its value is Θ ≈ 1.9 for M = 329 and Θ ≈ 1.55 for M = 658.

SUMMARY OF THE TOROIDAL SIMULATIONS

In this section we have presented the transition of the toroidal magnetohydrodynamic flow between an axisymmetric and a helical state. This transition appears for increasing pinch ratio and sufficiently large viscous Lundquist and Prandtl numbers. The helical state is characterized by different toroidal modes that dominate the system depending on the considered parameters. The kinetic and perturbed magnetic energies increase considerably when the instability is triggered. The root mean square average of the Lorentz force is almost invariant but its spatial distribution changes considerably. For increasing pinch ratio at the center of the domain a force free region appears but at the edge of the torus the magnitude of this force is considerably enhanced.

The analysis of the time evolution of the mean square velocity for each toroidal mode displays a complex behavior. Depending on the Θ parameter and the viscous Lundquist number the dynamics change dramatically. For a low M a non oscillatory steady state is obtained at the final times. On the other hand for the biggest viscous Lundquist numbers, strong oscillatory dynamics appear. This state is turbulent in the sense that we observe chaotic spatio-temporal fluctuations involving different modes, but the dynamics is still dominated by a single toroidal mode. It is therefore a quasi-single helicity state [START_REF] Escande | Quasi-single-helicity reversed-field-pinch plasmas[END_REF]. It is still an open question if multiple helicity states will be obtained with increasing viscous Lundquist number.

CYLINDRICAL SIMULATIONS

In this section the helical states of a cylindrical magnetofluid are presented. We will roughly consider the same quantities as in the preceding section on the toroidal simulations. This will allow us to highlight the differences and similarities between the dynamics in the two geometries. First the analysis of the dominant axial modes and the study of the flow topology as a function of the pinch ratio are carried out. Then the global quantities varying with Θ are shown. Finally the time evolution of the axial modes is analysed.

The imposed magnetic field for the cylinder is helical. It is composed of a poloidal linear magnetic field (B 0 P (r)) and a constant axial magnetic component (B 0 z ) (Fig. 4.24). The linear profile for the poloidal magnetic field comes from an imposed axial constant electric field in a system with constant resistivity (see Chap. II Sec. 5.2 and [START_REF] Shan | Nonlinear magnetohydrodynamics by Galerkin-method computation[END_REF]). The profile of the poloidal magnetic field is the following,

B 0 P (r) = σ E 0 z r 2 . (4.4)
The value of this field is fixed at the boundary of the cylinder and it is kept constant, B 0 P (a) = 0.32 (we recall that a is the cylinder radius). The reference magnetic field for the cylindrical calculations is fixed to B re f = B 0 P (a) = 0.32.

To vary the pinch ratio we modify the value of the axial magnetic field (B 0 z ) and the poloidal component is constant for the considered simulations. The initial velocity field is a small perturbation. It is a white noise with kinetic energy ≈ 4.3 • 10 -7 . This perturbation is necessary to trigger the helical instability. The perturbed initial magnetic field is zero. Note that no perturbation was needed in the toroidal geometry where the spontaneous appearance of the n = 0 mode allowed to trigger the transition.

The cylinder has a radius a = 1 and height L z = 4π, the curvature effects are neglected which is the major difference with the toroidal case. This corresponds to a "straight torus" where the aspect ratio, major radius over minor radius is Γ = 2. The size of the periodic box where the calculations are performed is (0.8π × 0.8π × 4π). We fix the penalization parameter to η = 5•10 -4 . The resolution for all the cylindrical calculations is 128 3 and we consider for both cases constant resistivity and viscosity. We vary the viscous Lundquist number (M) and we consider two Prandtl numbers (Pr). We evaluate the system as a function of the pinch ratio at the final state (Θ). The dimensionless numbers characterizing the dynamical system are defined as follows,

M = C A L ν , Pr = ν λ , Θ = B P (a) B z , F = B z (a) B z , (4.5)
with C A the poloidal Alfvén velocity. With respect to the toroidal simulations the main difference is the toroidal magnetic field which has been replaced by the axial magnetic field B z . As a reference length we take the diameter L = 2, ν is the viscosity, λ the magnetic diffusivity and a is the radius of the cylinder.

B0z B0P

Figure 4.24: Scheme of the helical magnetic field imposed for the cylindrical simulations.

The periodic cylindrical simulations are performed for various pinch ratios and two viscous Lundquist numbers, as in the previous section. We take into account the following dimensionless numbers M = 634 with a magnetic Prandtl number Pr = 3 and M = 1902 with Pr = 2. As in the previous section this choice is physically quite arbitrary but allows for a numerical investigation of a wide range of different pinch ratios.

Θ τ Θ = 2.7, M = 1902 Θ = 2.1, M = 634 Θ = 2.8, M = 634
F τ Θ = 2.7, M = 1902 Θ = 2.1, M = 634 Θ = 2.8, M = 634

HELICAL FLOWS

The time evolution of the reversal parameter and the pinch ratio are presented in Figs. [START_REF] Aydemir | Compressibility as a feature of field reversal maintenance in the reversed-field pinch[END_REF].25 and 4.26 respectively. The two quantities reach an asymptote, and there is a final steady state. For the highest viscous Lundquist number (M = 1902) there exist oscillations in the saturated state. This effect is quite similar as the one found in the toroidal simulations.

The F parameter evolution shows a decrease for the early times followed by a regrowth. This is linked to the dynamo action that enhances the axial magnetic field (see Fig. 4.27).

For the cylindrical case a bigger viscous Lundquist number value (M = 1902) is considered than for the study of the toroidal geometry. For this large M a reversal of the toroidal magnetic field is found if the the pinch ratio is increased sufficiently (see Fig. 4.26). The increase of the viscous Lundquist number seems to be important to make the toroidal magnetic field reverse at the edge. As in the previous toroidal study, the flow in the cylinder has a transition towards a helical configuration. The difference is that the state preceding the transition is not longer a toroidally invariant state (in the case of a cylinder we would say z-invariant or axially invariant state) but it is a zero velocity field, also called 'copper wire' solution [START_REF] Shan | Nonlinear magnetohydrodynamics by Galerkin-method computation[END_REF]. After the transition a helical deformation of the magnetofluid appears. It is characterized by different axial mode numbers (k z ). For the majority of studied cases the dominant axial mode is k z = 4. For some cases higher dominant modes are present such as, k z = 5 and k z = 7.

B ′ z τ Θ = 2.7, M = 1902 Θ = 2.1, M = 634 Θ = 2.8, M = 634
In Fig. 4.28 the ratio of the energy contained in the axial modes 0, 4, 5 and 7 over the total energy is presented as a function of the pinch ratio and for two different viscous Lundquist numbers. In both cases the zero axial mode has a negligible fraction of the total energy for all the pinch ratio values considered.

For the low viscous Lundquist number, M = 634 (Fig. 4.28 (left)) there are always one or two modes that possess more than 30% of the total energy. Therefore this system is dominated by a few helical modes. If M = 1902 (Fig. 4.28 (right)) at low pinch ratio (Θ = 0.9) the mode k z = 4 dominates clearly with more than 50% of the total energy. Then, for a pinch ratio greater than 2.2, the energy of the considered modes converges, in each mode the energy is below 20% of the total kinetic energy. This means that the energy spectrum is more flat and we have a state that could be described as a multiple helicity state.

To illustrate these helical states we show in figure 4.29 a positive and a negative isocontour of axial velocity for two different values of the pinch parameter. For the case with Θ = 2.1 the fourth axial mode dominates the final state. On the other hand for Θ = 2.8 the dominant mode is k z = 5. The two axial velocity isocontours share the same magnitude but are in opposite directions. These oppositely signed isocontours have a similar topology, and they are separated by a short axial distance.

As for the toroidal simulations there is a correlation between the helical deformation of the perturbed axial magnetic field and the dominant toroidal velocity modes. This is visualized in Fig. 4.30 where we present a perturbed toroidal magnetic isocontour colored by the axial velocity. In this figure the deformation of the magnetic field (the shape of the If we now look at the kinetic and axial magnetic spectra (Figs. 4.31 and 4.32) we notice that the maximum values are located at the same axial mode numbers. In the case with Θ = 2.1 (Fig. 4.31 (left)) the dominant mode for the kinetic field is k z = 4 followed by the third mode k z = 3. For the axial magnetic field (Fig. 4.31 (right)) the dominant mode is k z = 3 but the second in magnitude is k z = 4. The same two dominant modes in the system appear both in the kinetic and in the axial magnetic fields spectra.

For a larger value of the pinch ratio (Θ = 2.8) the similarity between the kinetic and magnetic spectra is more pronounced (see Fig. 4.32). The dominant and the secondary modes are similar. The biggest axial modes correspond to small values of k z , varying from k z = 3 to k z = 5, with a strong predominance of the fifth mode for the two considered fields.

The velocity field has helical perturbations and varies as a function of the pinch ratio. To illustrate this variation we take at one instant, in the steady state, the average in the z direction of the axial and poloidal plane velocities. The Fig. 4.33 shows this velocity average, two distinct flow topologies as a function of Θ can be distinguished. For the lowest value of the pinch parameter (Fig. 4.33 (a)) the axial average shows small velocities compared to the others two figures ((b) and (c)). This is a consequence of the strong dominance of the axial mode k z = 4 for the case (a). It was shown in Fig. 4.28 that the energy contained in the fourth mode for the pinch ratio Θ = 0.9 is almost 60% of the total energy, so that we have a laminar quasi-single helicity state. Previously we saw that a positive axial velocity with a helical structure comes along with a negative one (Fig. 4.29). In this laminar helical state, the axial averaging adds both contributions up to zero.

The flow topology of the cases (b) and (c) in Fig. 4.33 is similar but its axially average is qualitatively different. There exist positive and negative axial velocity regions, but the flow is more chaotic for the largest value of the pinch ratio (Θ = 2.7). A central circular zone can be distinguished where the flow is oriented in one axial direction and where the poloidal velocities are relatively small.

EVOLUTION OF GLOBAL QUANTITIES

The time evolution of the kinetic and perturbed magnetic energy are presented in Fig. 4.34 for a viscous Lundquist number M = 1902. There is a similar growth in both quantities when the instability is triggered. The increase in the pinch parameter make the two energies oscillate. At the saturated state there are fluctuations for the calculations with Θ = 1.88 and Θ = 2.70. The evolution of the kinetic and perturbed magnetic energy with the pinch parameter for the two considered M numbers is almost linear (Fig. 4.35). The kinetic energy for the largest value of Θ is four times bigger than for the lowest. The influence of the viscous Lundquist number is not very important in the evolution of the kinetic energy. On the other hand in the perturbed magnetic energy the computations show different evolutions depending on the considered M number. There is a linear evolution but at high pinch ratio the energy for M = 1902 is twice as big as for M = 634. The dynamo action seems to be stronger for large viscous Lundquist values. Given that the kinetic energy is comparable, this seems logical since the magnetic perturbation will be more smoothened out by higher diffusion at lower viscous Lundquist number. The evolution in the kinetic and perturbed magnetic energies shows the same tendency as the evolution of the root mean square of the total Lorentz force in the volume (see Fig. 4.36). The Lorentz force has a linear evolution with Θ, and there is a growth of ∼ 50% between the lowest and the highest pinch parameter value.

The current helicity is presented in Fig. 4.37. This quantity is far from unity, indicating the absence of alignment between the total magnetic and current density fields. Both fields have a tendency to align with increasing pinch ratio, as indicated by the current helicity value that grows, but this change is relatively small. The spatial distribution of the norm of the Lorentz force is presented in Fig. 4.38 for three different pinch ratios and constant viscous Lundquist number (M = 1902). As for the toroidal simulations the field tends to be force free in the center of the geometry. There is also an important increase of the value of its magnitude near the boundaries. Notice also that for the case with Θ = 2.7 helical perturbations in the center of the domain are visualized. These will probably average out if in addition to the axial average a time average would be applied.

TIME EVOLUTION OF THE DIFFERENT AXIAL MODES

The evolution of the square velocity for the different axial mode numbers is presented in Fig. 4.39 for M = 634 and in Fig. 4.40 for M = 1902. The axial mode k z = 0 is not presented because for all the considered simulations the square velocity magnitudes for this mode are several orders of magnitude smaller than for the other modes exposed here.

In the first figure (Θ = 0.93) there is clear temporal frequency that dominates the system with a constant period of ∼ 100 poloidal Alfvén unit times. The dominant modes in this case are: k z = 4 and k z = 5. For increasing pinch ratio the temporal frequency vanishes and there is a steady predominance of the fourth axial mode (Θ = 2.09).

With a further increase in the pinch ratio (Θ = 2.76) the dominant mode changes, it is not longer k z = 4 but k z = 5. For the largest value, Θ = 3.16, a higher mode (k z = 7) starts dominating the system together with the lower mode, k z = 4.

The system at this viscous Lundquist number shows dynamics different from the toroidal geometry. The mode k z = 4 is only one of the dominant modes. A difference with the torus case is the appearance of higher mode numbers such as k z = 5 and k z = 7. For the present viscous Lundquist number and pinch ratios there are no chaotic temporal oscillations.

Now if we look at the time evolution of the axial modes for the highest viscous Lundquist number (M = 1902) completely different dynamics appear (Fig. 4.40). For

u 2 τ k z = 1 k z = 2 k z = 3 k z = 4 k z = 5 k z = 6 k z = 7 k z = 8 Θ = 0.93 10 -9
10 -8 10 -7

10 -6 10 -5

10 -4 10 -3 0 400 800 1200 1600 the low pinch parameter considered (Θ = 0.90) the dynamics are dominated by the fourth mode. However if the pinch ratio is increased oscillations around the saturated state appear (Θ = 1.88). For this case the dominant mode is still k z = 4, but it has an oscillatory behavior along with the others toroidal modes. For the highest pinch ratios (Θ = 2.38 and Θ = 2.70), oscillations are presented for a large spectrum of modes. The difference is that the dominant mode is not unique anymore, the amplitudes of the oscillations are larger and, in time, different modes possess the biggest kinetic energy. The dominant axial modes are 4, 5, 6 and 8 for the case with Θ = 2.38. For the greatest pinch parameter (Θ = 2.70) the dominant ones are 4, 5 and 7. The kinetic spectrum becomes flatter and we can say that a fully turbulent state starts to emerge.

u 2 τ k z = 1 k z = 2 k z = 3 k z = 4 k z = 5 k z = 6 k z = 7 k z = 8 Θ = 2.09
u 2 τ k z = 1 k z = 2 k z = 3 k z = 4 k z = 5 k z = 6 k z = 7 k z = 8 Θ = 2.
u 2 τ k z = 1 k z = 2 k z = 3 k z = 4 k z = 5 k z = 6 k z = 7 k z = 8 Θ = 3.16
u 2 τ k z = 1 k z = 2 k z = 3 k z = 4 k z = 5 k z = 6 k z = 7 k z = 8 Θ = 0.
u 2 τ k z = 1 k z = 2 k z = 3 k z = 4 k z = 5 k z = 6 k z = 7 k z = 8 Θ = 1.
u 2 τ k z = 1 k z = 2 k z = 3 k z = 4 k z = 5 k z = 6 k z = 7 k z = 8 Θ = 2.
u 2 τ k z = 1 k z = 2 k z = 3 k z = 4 k z = 5 k z = 6 k z = 7 k z = 8 Θ = 2.70

COMPARISON BETWEEN THE RESULTS IN TOROIDAL

AND CYLINDRICAL GEOMETRY

INTERCOMPARISON OF THE SIMULATIONS

The two studied geometries show similarities in the final dynamics of the magnetofluid. Both cases at the final times develop helical configurations if the pinch ratio is sufficiently increased, and the fourth toroidal (or axial) mode prevails in the majority of computed helical cases.

The fundamental difference appearing between the two geometries is the absence of the zero mode for the cylindrical cases. For the toroidal geometry, at low pinch parameter, the flow is an axisymmetric velocity field that with increasing Θ has a transition to a state, where a dominant helical flow is superposed on the original axisymmetric field. To illustrate the difference between the flow in the cylinder and the torus we take one simulation in each geometry. For the velocity field we perform a time average at the steady state and then a second spatial average in the azimuthal direction for the toroidal geometry and in the axial direction for the cylindrical case. A time average is performed because for both considered cases at the steady state we have an oscillatory evolution of the magnetofluid. The image resulting of this average is presented in Fig. 4.41 for the toroidal geometry (M = 658, Pr = 3 and Θ = 1.57) and in Fig. 4.42 for the cylinder (M = 1902, Pr = 2 and Θ = 1.88). In the toroidal result we can visualize clearly the two axisymmetric regions with opposite toroidal velocity. The maxima are located near the boundaries and surround a circular region with small toroidal velocity. On the other hand, in the cylindrical case, the time and axial average give an axial velocity field which is negligible compared to the toroidal case. Another way to apprehend this difference between the toroidal and cylindrical flow is 88 to see the ratio of the kinetic energy of the zero mode over the total kinetic energy (see Fig. 4.43). For M = 329 in the toroidal configuration the energy is mainly concentrated in the zero toroidal mode at low pinch ratio, after the transition between Θ = 0.65 and Θ = 1.44 the zero mode energy falls but does not reach a negligible value. The same behavior is observed for the toroidal case with M = 658, where for the greatest Θ the energy in the zero toroidal mode is small, but still constitutes several percent of the total kinetic energy. For the cylindrical case for all the considered pinch ratios the energy of the zero axial mode is always negligible compared the other modes of the spectrum. It is always inferior to 1.2% of the total energy. The time evolution of the kinetic energy also shows a different trend. For the toroidal case rapidly a kinetic energy appears at the poloidal Alfvén time τ = 1 when the axisymmetric velocity field emerges. Subsequently the helical instability arises and makes the kinetic energy increase around τ = 5. For the cylinder case we have also the instability that develops after τ = 5 but before that there is no increase of the kinetic energy, the energy remains below 10 -5 . This is the energy of the initial perturbation that decays to zero if the threshold of the helical instability is not reached.

COMPARISON WITH EXPERIMENTAL DATA

The numerical results are compared to experimental data of three different RFP devices. The first experimental data comes from the REPUTE experiment, that is described in [START_REF] Taylor | Relaxation and magnetic reconnection in plasmas[END_REF][START_REF] Toyama | Unknown[END_REF], the second set of data is the RFP ZT-40M [START_REF] Taylor | Relaxation and magnetic reconnection in plasmas[END_REF][START_REF] Dimarco | The ZT-40M experiment[END_REF] and the third is from the device RELAX [START_REF] Ikezoe | Extended operational regimes and MHD behavior in a low-aspect-ratio reversed field pinch in RELAX[END_REF]. This last experiment has a low aspect ratio, Γ = 2, close to the aspect ratio used in the present simulations. This device has the particularity to allow Θ to be greater than two, whereas for the majority of RFP experiments the pinch ratio is around 1.5. In RELAX a distinctive deep-reversal configuration of the toroidal magnetic field is found where the flow is characterized by a low level of fluctuations [START_REF] Ikezoe | Extended operational regimes and MHD behavior in a low-aspect-ratio reversed field pinch in RELAX[END_REF][START_REF] Ikezoe | Characterization of equilibria in a low-aspectratio RFP[END_REF]. In Fig. 4.45 the numerical and the experimental results are presented in the Θ-F plane. We also show the theoretical curve found by Taylor [START_REF] Taylor | Relaxation of toroidal plasma and generation of reverse magnetic fields[END_REF] that has been recently critically reviewed in [START_REF] Cappello | The reversed field pinch toward magnetic order: a genuine self-organization[END_REF] with respect to its relevance to the RFP dynamics.

From the figure we can see that the numerical calculations give similar results compared to the RFP experiments in the Θ-F plane: we have the decrease of the F parameter with the increase of the pinch ratio (Figs. 4.45 and 4.46). The reversal of the toroidal magnetic field B T for the considered experiments occurs around the same value, for Θ ≈ 1.5. The numerical calculation that fits the best the experiments is the one performed for a toroidal geometry with Pr = 3 and viscous Lundquist number M = 658. In this figure we note that the geometry and the viscous Lundquist number play an important role in the evolution of the reversal parameter F with Θ. We see that increasing the viscous Lundquist number, the value of Θ, that corresponds to F = 0, decreases. It seems that, using the toroidal geometry and increasing sufficiently the viscous Lundquist number as well as the magnetic Prandtl number, we fit better the experimental data. We do not know at present whether we could get a magnetic toroidal reversal at the torus edge for high enough values of M. In this study the computational resources were not sufficient to carry out higher resolution computations in order to investigate higher values of the Lundquist and Prandtl numbers. It will be possible in a close future to perform such simulations and to check if we can reach the region, in parameter space, where the toroidal magnetic field reverses.

CONCLUDING REMARKS ABOUT THE TOROIDAL AND CYLINDRICAL SIMULATIONS

The toroidal and cylindrical simulations carried out in the present study help to understand the similarities and differences that appear between these two geometries. The helical instability appearing for increasing pinch ratio is quite similar for both cases. Generally, if the viscous Lundquist number is not too important, a single helical mode dominates the system, so that we are close to the quasi-single helicity state found in experiments. For increasing viscous Lundquist number an oscillatory steady state is found in both geometries. A turbulent state appears for the cylindrical case at large M number and pinch ratio. It can be expected that a similar turbulent state will appear in the toroidal configuration if the viscous Lundquist number is sufficiently increased.

The fundamental qualitative difference between the cylinder and the torus is the axisymmetric or axially invariant flow that disappears in the cylindrical configuration. The toroidally invariant poloidal flow is composed of two counterotating vortices. This flow topology has been found experimentally [START_REF] Bonomo | Flow measurements and modelling in helical RFX-mod equilibria[END_REF], and in a straight cylinder this important characteristic of the plasma dynamics is lost.

In summary, despite some qualitative similarities between the dynamics in the two geometries, we stress that quantitatively the differences are quite large. We believe that the curvature of the toroidal geometry needs to be taken into account to correctly reproduce the experimental flow measurements.

V CONCLUSION

The objective of this thesis project was to study the dynamics of a magnetofluid in a toroidal geometry by the means of the incompressible, isothermal, viscoresistive magnetohydrodynamic model with constant transport coefficients. Using this simplified model many characteristics of a real fusion plasma are not taken into account. Hence we study the dynamics of a plasma without many of its complexities. The advantage of this approach is that basic phenomena can be easily identified. Naturally more realistic characteristics could be introduced later on, one at a time, and a separation of the different phenomena that influence a more realistic plasma dynamics could then be highlighted. We believe that, in this way, a better understanding of the dynamical system can be achieved.

First, in Chap. II, a novel numerical method to solve the MHD equations in toroidal geometries is presented, where a forcing term, added to the Navier-Stokes and the magnetic induction equations, allows to take into account Dirichlet boundary conditions in a pseudo-spectral algorithm. The advantage of the method is that the considered geometry can easily be modified. As an example, in this manuscript the MHD equations were solved in different toroidal and cylindrical configurations.

In Chap. III, where we applied this method to investigate the behavior of MHD in toroidal geometry, we showed that the presence of an irrotational toroidal electric and magnetic field in a toroidal configuration automatically produces velocities even in the absence of MHD instabilities. Therefore under the assumptions we used, no static equilibrium is possible. If one would like to describe toroidal plasma dynamics as perturbations around an equilibrium state, this state must therefore be necessarily dynamic. The velocities for high values of the viscosity and resistivity are mainly located in the poloidal plane, but for increasing nonlinearity and for a large safety factor the toroidal velocity component considerably increases and dominates the dynamics. Also it is observed that if an up-down asymmetric cross section of the torus is considered, a net toroidal angular momentum appears, a feature which is also observed in experiments.

For a toroidal magnetofluid, where the imposed helical magnetic field is characterized by a small safety factor (or large pinch ratio), more complex dynamics are found. In this case a transition from an axisymmetric flow field to a helical configuration is observed. This helical topology resembles the quasi-single helicity state found in RFP experiments. For increasing nonlinearity and pinch ratio a transition to a fully turbulent system can be achieved. It is shown, that even though some similar features are present in both, the cylindrical and toroidal geometry, important differences exist between the flow topology observed in them, both qualitatively and quantitatively. Also the comparison with experimental data shows a better agreement if the toroidal geometry is considered.

PERSPECTIVES

The physics of toroidally confined fusion plasmas is very rich and almost every aspect of it can constitute a research field on its own. The MHD model we use can be improved in many ways to approach reality. However, one should be careful not to exaggerate this complexification. Indeed the goal of a description such as MHD is not, in our opinion, to build a numerical fusion reactor but to isolate key physical mechanisms that can help to better understand the plasma dynamics. This understanding can then guide experimentalists in interpreting experimental results in order to improve magnetic fusion configurations. Considering too many features at once can perhaps diminish the gap between numerical and experimental observations, but does not necessarily lead to better understanding. With this in mind we propose here several perspectives, favouring the ones which are most simple to consider and which, possibly, have the biggest impact on the observed plasma behavior.

In the present manuscript, the magnetic Prandtl number was kept close to unity. Experimentally it is difficult to assess the real viscosity value of a plasma, so it seems interesting to perform a parametric study varying the Prandtl number over several orders of magnitude to quantify its influence on the dynamics of the system. This is one of the practically most simple perspectives since it does not need any modification of the numerical code, as long as sufficient numerical resources are available.

The influence of the shape of the domain on the RFP dynamics and in particular on the appearance of quasi-single helicity states deserves definitely further investigation.

The model considered in the present study can further be developed. The first step could be the introduction of a profile for the resistivity. In this way, the dependence of the resistivity on the temperature can be taken into account. A more realistic way to include the influence of non-uniform transport coefficients is to include an evolution equation for the temperature. This, and including the effects of finite compressibility, will necessitate a non negligible effort, notably with respect to the introduction of the penalization term in the energy and in the continuity equations [START_REF] Liu | A Brinkman penalization method for compressible flows in complex geometries[END_REF].

Another perspective, which is not directly related to the field of magnetic fusion, is the use of the numerical code to study less complex plasma experiments that can be used for a more fundamental understanding of the plasma physics. An example is the preliminary study of a plasma dynamo presented in Appendix F.

where α = √ η and e * = e 1 √ η . The subscripts 's' and 'f' stand for the solid and the fluid domain, respectively.

For the penalization with the regularization method, the interpolating polynomial must be defined first. An Hermite interpolation consists in finding a polynomial which fits two separate points where the values of the function and of a chosen number of derivatives are imposed. In our case the Hermite polynomial H(x) must match the values of u and its first derivative at x = 0 and yields H(1) = H ′ (1) = 0 at x = 1. As the exact solution of the unpenalized problem is u(x) = x+1 2 , the expression of the third degree Hermite polynomial can be easily derived:

H(0) = 1 2 , H ′ (0) = 1 2 , H(1) = 0, H ′ (1) = 0 H(x) = 3x 3 2 -5x 2 2 + x 2 + 1 2
In equation (A.1), H(x) is set as u wall and, with the same notations as before, the solutions are The penalization error in the L 2 norm can be computed directly as

       u f (x) = a(x + 1), u s (x) = K 1 e x √ η + K 2 e -x √ η + H(x) + (9x -5)η with a = 14α 3 (1-e 2
||ε|| 2 2 = 0 -1 (u η -u exact ) 2 dx = 0 -1 (a(x + 1) - x + 1 2 ) 2 dx = (1 -2a) 2 12
Then the coefficient a is replaced by the corresponding expressions obtained above in the two cases and the error becomes :

       ||ε|| 2 = | (1-2a) √ 12 | classical penalization : ||ε|| 2 = √ η √ 12
hermite penalization :

||ε|| 2 = 10η √ 12
The convergence results are confirmed by the numerical experiments shown in Fig. A.4. The penalized solutions are plotted for several values of η = [10 -1 , 10 -2 , 10 -3 , 10 -4 , 10 -5 ] and the errors with respect to the exact profiles are computed in the L 2 and L ∞ norms. The order of convergence is indeed improved to η for the regularization method instead of √ η with the classical penalization. The errors are also presented for the first and second derivative (Fig.A.5). It must be noticed that the error for the second derivative in the case of the regularized penalization method still improves when η decreases, whereas it increases for the classical method.

The order of convergence obtained in Chap. II sections 5.1 and 5.2 is higher than the one obtained here. This could result from the dissipative terms in the Navier-Stokes and induction equations, which are not present here. This term smoothens the remaining discontinuities and improves the order of convergence. The principal drawback of this method is that a baseflow must be known to precompute the regularization term. Moreover, it is possible that this baseflow is not continuous if in the numerical solution instabilities appear. As further perspective, this method should be implemented actively so that the regularization term fits the numerical solution as closely as possible at each time-step.

B IMPOSING NON HOMOGENEOUS NEUMANN BOUNDARY CONDITIONS WITH A PENALIZATION METHOD

FORMULATION

A generalization of the volume penalization method to impose non homogeneous Neumann boundary conditions (BC) will be explained in the following.

The advection-diffusion equation of a scalar field θ (t, x) is considered in a domain Ω,

∂ θ ∂t + u j ∂ j θ = α∆θ , (B.1) ∇θ | ∂ Ω = γ wall , θ (t = 0, x) = f 0 (x). (B.2)
Taking the gradient of this equation we can obtain a system of advection-diffusion equations for the different components of the gradient of θ , denoted by γ i = ∂ i θ . The penalization term is added and inhomogeneous Dirichlet conditions are imposed on γ i ,

∂ γ i ∂t + ∂ i (u j ∂ j θ ) = α∆γ i - χ(x) η (γ i -γ i wall ) . (B.3)
We use a Fourier pseudo-spectral method and thus periodic boundary conditions are applied at the edge of the computational domain for the imposed γ i wall . The next step is to integrate this last equation in space to recover the equation for θ . To do so we apply the divergence operator to the system of equations (B.3). Using the relation ∇ • (∇ f ) = ∆ f , the Laplace operator appears in the equation. Thus to recover the equation for θ a Poisson equation has to be solved. In spectral space this is equivalent to multiply the equation (B.3) by the operator

ık i /k 2 , θ (k) = -ık i k 2 γi (k). (B.4)
with k = 0. In consequence, in physical space, the advection-diffusion equation with the penalization term for the gradient, can be written:

∂ θ ∂t + u j ∂ j θ = α∆θ + F -1 ık i k 2 F χ(x) η (γ i -γ i wall ) . (B.5)
This last equation can be solved numerically with a pseudo-spectral method and the solution of the advection-diffusion equation with non homogeneous Neumann boundary conditions can thus be obtained.

To verify this numerical method a one-dimensional test case is presented in the next section considering an instationary diffusion problem with non homogeneous Neumann boundary conditions. A different numerical method to take into account homogeneous boundary conditions using spectral methods is presented in [START_REF] Kadoch | A volume penalization method for incompressible flows and scalar advection-diffusion with moving obstacles[END_REF].

TEST CASE: ONE-DIMENSIONAL INSTATIONARY DIF-FUSION PROBLEM

We consider the instationary diffusion equation for the temperature θ (x,t), with a positive heat flux on the left and right boundaries (see Fig.

∂ θ ∂ x x=a = -q and ∂ θ ∂ x x=b = q. (B.7)
The solution of this system is a combination of a parabolic function in space and a linear time evolution,

θ (x,t) = q(x -a) (x -a) (b -a) -1 + 2q (b -a) t. (B.8)
The initial condition chosen for the simulation is the parabolic function that respects the inhomogeneous boundary conditions,

θ (x, 0) = q(x -a) (x -a) (b -a) - 1 .
(B.9) On the border of the domain the periodicity must be ensured and therefore a zero gradient of θ is imposed at x = 0 and x = L. For this reason we impose a linear evolution of the gradient (γ i wall (x)) in the solid region,

γ i wall (x) x∈[0,a] = -q a x, γ i wall (x) x∈[b,L] = q L -b (L -x). (B.11)
Hence we verify the Neumann boundary conditions at x = a and x = b, as well as the periodicity at the edge of the computational domain. [START_REF] Jause-Labert | Numerical validation of the volume penalization method in three-dimensional pseudo-spectral simulations[END_REF] extended an implicit formulation for the time-integration scheme for Dirichlet conditions, originally proposed by Kolomenskiy and Schneider [START_REF] Schneider | A Fourier spectral method for the Navier-Stokes equations with volume penalization for moving solid obstacles[END_REF] for Burgers' equation, to relax the constraint on the choice of the time-step. Their approach is discussed in some detail in this section. We use the following form of the penalized Navier-Stokes equation:

∂ t u = ξ -∂ x P - χ η u, (C.1)
in this equation u wall = 0, P + u 2 /2 → P and ξ represents the nonlinear term. The viscous term can simply be added in the integral form, u exp(-νk 2 t) → u. In Fourier space we can write ∂ t û = P ⊥ ξ - The semi-implicit penalization method is in this limit analogous to an explicit formulation. The difference is that, if the penalization parameter is small enough, the porosity of the solid walls is given by the value of the time step.

D CALCULATION OF THE IMPOSED POLOIDAL MAGNETIC FIELD

We construct the magnetic field satisfying the following properties: (i) it corresponds to a current density profile ∝ 1/R, (ii) it is parallel to the wall and (iii) it is solenoidal. With respect to our previous investigation [START_REF] Morales | Intrinsic rotation of toroidally confined magnetohydrodynamics[END_REF] the magnetic topology is changed. In fact in that investigation the imposed poloidal magnetic field satisfied the imposed toroidal current density profile j 0 T and the solenoidal constraint, but the normal component did not vanish (as is shown in Here λ is a diffusion coefficient (λ = 1), Ξ is the mask function (it takes the value one in the region where the Poisson equation needs to be solved and zero in the rest of the computational domain) and η is the penalization parameter (η = 5 • 10 -4 ). The size of the domain is (2π) 3 for the asymmetric and (2π × 2π × π) for the symmetric cross section. The value of the Dirichlet boundary condition is χ wall . This equation is evolved in time, reaching a steady state, numerically χ n+1 -χ n < 10 -6 . We then recover with sufficient accuracy the solution of the Poisson equation (D.3) taking into account the Dirichlet boundary condition via the penalization term. The solution of this precomputation will give our basis magnetic field B 0 which will be kept constant during the actual simulation. The resulting poloidal magnetic topology is presented in

E RELATION OF THE DIMENSIONLESS MHD EQUATIONS TO PHYSICAL VALUES FROM THE JET TOKAMAK

The justification of the choice of the parameters we use in our simulations with respect to the experimental values in real life fusion reactors is a difficult exercise, in particular since our assumptions of uniform, scalar transport coefficients are necessarily gross oversimplifications. However, here we present an attempt to estimate the values we would like to use ideally if we would aim to approach physical reality as closely as possible. This estimation is based on the presentation in reference [START_REF] Kamp | Toroidal steady states in visco-resistive magnetohydrodynamics[END_REF].

To determine the dimensionless values that should be introduced in the numerical code we consider the values of the JET experiment [START_REF] Jet-Team | Fusion energy production from a deuterium-tritium plasma in the JET tokamak[END_REF]. In this device there is a toroidal magnetic field of 28 kG (in SI units B JET = 2.8T), a toroidal current I = 3.1 MA, the minor radius is a = 1.5 m and major radius is R 0 = 3 m. Also we have the plasma density 3.6 • 10 19 m -3 . For our numerical simulations the evolution along the big radius of the imposed dimensionless toroidal electric end magnetic fields are: are going to consider its parameters in SI units. The root mean square value of the toroidal magnetic field is taken to be B JET = 2.8T, this is the reference magnetic field.

E * 0 (R) = Ẽ * R 0 R e T , (
Note that the reference value of B in our simulation is chosen to be the toroidal magnetic field at R = R 0 . The value of JET is the surface averaged magnetic field. This introduces a first correction. To find the dimensionless B * value the integration of B * 0 over the surface of the circular cross section must be equal to one, with a = 1.5 m, R 0 = 3 m and A the area of the cross section (considering a circular cross section A = πa 2 ). We calculate this integral using the numerical software Maple. For the circular cross section we find the same value as Kamp and Montgomery [START_REF] Kamp | Toroidal steady states in visco-resistive magnetohydrodynamics[END_REF], which is B * = 0.94. Numerically if the value of the magnetic field is unity it will be in reality B JET = 2.8T. Note that the correction induced by the integral (E.3) is of order unity, which is very small compared to the uncertainties of the transport coefficients.

To find the value of Ẽ * we impose a toroidal loop voltage V T = 1V . To have the corresponding electric field in the cross section V T is divided by the toroidal length,

E section = V t 2πR 0 . (E.4)
To make this value dimensionless we divide E section by C a B JET (E * = E section / C a B JET ).

With C a = B JET / √ ρ µ the Alfvén velocity (C a = 7.2 • 10 6 ms -1 ), we recall that ρ is the density and µ is the permeability of free space. The integration of E * 0 over the cross section must be equal to the dimensionless value of the electric field E * . As for the toroidal magnetic field we write, The characteristic length is L = a/(0.3π) because in the numerical code the minor radius has a length of 0.3π (a = 1.5 m is the real length of the minor radius for JET).

E * = 1 A
The viscous Lundquist number is defined as M = C a L/ν. The high boundary of the value of the kinematic viscosity is ν = 1.25 • 10 11 m 2 s -1 . In fact for a magnetofluid exposed to a strong magnetic field a complicated viscous stress tensor results [START_REF] Braginskii | Transport processes in a plasma[END_REF] with different viscosity coefficient that span about twelve orders of magnitude [START_REF] Montgomery | MHD steady states as a model for confined plasmas[END_REF].

F LARGE SCALE FORCING OF A PLASMA DYNAMO

The generation of planetary magnetic fields through the turbulent movement of liquid metals in the core is an intriguing phenomenon. It is not only an important issue in geophysics, but also explains the self-organization in Reversed Field Pinch plasma devices. This so called dynamo effect is now reproduced in several experimental set-ups such as the VKS experiment in Cadarache [START_REF] Monchaux | Generation of a magnetic field by dynamo action in a turbulent flow of liquid sodium[END_REF], but the detailed physics are still poorly understood. One drawback in the liquid metal experiments is the very low value of the magnetic Prandtl number. Liquid metal experiments require therefore a huge kinetic Reynolds number to obtain only a moderate magnetic Reynolds number. For this reason the construction of a plasma dynamo has been undertaken at the Laboratoire de Physique of the Ecole Normale Supérieure de Lyon. The use of a plasma instead of a liquid metal allows the experiment to attain larger magnetic Reynolds numbers at moderate kinetic Reynolds numbers. The measurement of the dynamics in a plasma experiment, and the forcing of the large scale flow field, is however highly nontrivial. Therefore, in parallel with the experiment, we investigate the plasma dynamics by direct numerical simulations § .

In this work we model the plasma vessel by imposing solid, no-slip boundary conditions at the walls of a cylinder, periodic in the axial direction. An imposed magnetic field B 0 is pointing in the axial direction. An electric field is created in the radial direction as illustrated in Fig.

Figure 1 . 1 :

 11 Figure 1.1: Toroidal geometry surrounded by helical magnetic field lines.

ΩfFigure 2 . 1 :

 21 Figure 2.1: The computational domain Ω contains both the fluid domain Ω f and the solid domain Ω s .

Figure 2 . 2 :

 22 Figure2.2: Leaving a part of the domain unpenalized at the edges of the computational domain allows to use inhomogeneous Dirichlet conditions at the fluid solid interface, without violating the periodicity condition at the edge of the domain (left). Another option is to interpolate the imposed velocity to a zero value with a horizontal tangent at the domain frontier with an interpolating Hermite polynomial. The velocity field in the whole computational region belongs then to the C 1 class (right).

Figure 2 . 3 :

 23 Figure 2.3: Code scaling in the calculator Vargas of IDRIS for hydrodynamic and MHD calculations.

Figure 2 . 4 :

 24 Figure 2.4: Taylor-Couette flow configuration.

2 Figure 2 . 5 :

 225 Figure 2.5: Taylor-Couette: convergence of the relative L 2 error of u θ with the penalization parameter η, convergence order √ η.

Figure 2 . 6 :

 26 Figure 2.6: Taylor-Couette: convergence of the relative L 2 error of u θ with the resolution N which shows a second order convergence.

  2.6).

3 Figure 2 . 7 :

 327 Figure 2.7: Taylor-Couette: convergence of the relative L 2 error with the penalization parameter η, semi-implicit (solid line) and explicit (dashed line) methods.

Figure 2 . 8 :

 28 Figure 2.8: Taylor-Couette: convergence of the relative L 2 error with the resolution N, semiimplicit (solid line) and explicit (dashed line) methods.

Figure 2 . 9 :

 29 Figure 2.9: z-pinch configuration.

N 4 η 1

 41 = 256 -Hermite N = 512 -Hermite N = 256 -No Hermite N = 512 -No Hermite η

Figure 2 .

 2 Figure 2.10: Z-pinch: convergence of the relative L 2 error with the penalization parameter η (left) and with the resolution N (right) for the magnetic field in the z-pinch geometry. Comparison of the results with and without Hermite polynomial interpolation.

  29) and B(x, y, z,t = 0) = β [-2 sin(2y) + sin(z), 2 sin(x) + sin(z), sin(x) + sin(y)] (2.30) with β = 0.8, the initial kinetic energy, E k = 2, and the corresponding magnetic energy, E m = 1.92. The energies are evaluated by E k = 1 2 u 2 and E m = 1 2 B 2 with .. the volume average. The maximum of the current density is calculated by

= λ = 1 •

 1 10 -3 ν = λ = 5 • 10 -3 ν = λ = 1 • 10 -2

7 Figure 2 . 14 :Figure 2 . 15 :

 7214215 Figure 2.14: Evolution of the axial kinetic energy for different Re numbers.

Figure 2 . 16 :

 216 Figure 2.16: Axial velocity u z (color) for Re = 120.

Figure 2 . 17 :

 217 Figure 2.17: Azimuthal vorticity ω θ (color) for Re = 120.

Figure 2 . 18 :

 218 Figure 2.18: Azimuthal velocity u θ (color) for Re = 120.

Figure 2 . 19 :

 219 Figure 2.19: Radial velocity u r (color) for Re = 120.

= 15 Figure 2 . 20 :

 15220 Figure 2.20: Evolution of the axial kinetic energy for different Hartmann numbers.

Figure 2 . 21 :

 221 Figure 2.21: Hydromagnetic flow. Cut in the rz plane of velocity u (left) and magnetic (B) field (right) for Pr = 1, Ha = 7, Re = 100, ζ = 0.5 and Γ = 4.

zFigure 2 . 22 :

 222 Figure 2.22: Helical magnetic field scheme.

Figure 2 .

 2 Figure 2.23: First helical mode, velocity streamlines colored with the axial velocity u z .

Figure 2 .

 2 Figure 2.25: Zoom: average axial current as a function of the average electric field.

Figure 2 .

 2 Figure 2.26: Zoom: total energy dissipation rate as a function of the average electric field.

Figure 2 . 27 :Figure 2 . 28 :

 227228 Figure 2.27: Average axial current as a function of the average electric field.

Figure 3 . 1 :

 31 Figure 3.1: Cross-sections of the toroidal geometries considered in the present work. The toroidal direction is labelled T and the poloidal P.

Figure 3 . 2 :

 32 Figure 3.2: Three-dimensional magnetic field lines colored with the vertical magnetic field (B Z ). For the symmetric (left) and asymmetric (right) cross sections.

Figure 3 . 3 :

 33 Figure 3.3: Streamlines colored with toroidal velocity (u T ) for M = 23.

- 1 • 10 - 3 - 5 • 10 - 4 0 5 • 10 - 4 1

 110351045104 • 10 -3

Figure 3 . 4 :

 34 Figure 3.4: Azimuthally averaged toroidal velocity and poloidal stream function contours (solid line positive, dotted line negative contours).

Figure 3 . 5 :Figure 3 . 6 :

 3536 Figure 3.5: Azimuthally averaged: (a) Toroidal vorticity ω T and poloidal stream function, (b) radial velocity u R and (c) perturbation of the toroidal magnetic field, B ′ T .

2 rFigure 3 . 7 :

 237 Figure 3.7: Poloidal and toroidal square velocity component evolutions at early times, in toroidal Alfvénic time units (τ).

1 • 10 - 7 1 • 10 - 6 1 • 10 - 5 1 • 10 Figure 3 . 8 :

 1071061051038 Figure 3.8: Kinetic energy evolution at large times (left) and oscillatory behavior at early time (right) in toroidal Alfvénic time units, for asymmetric and symmetric geometry.

1 • 10 - 6 1 • 10 - 5 1 • 10 - 4 1 • 10 Figure 3 . 9 :

 1061051041039 Figure 3.9: Fluctuating magnetic energy evolution at large times (left) and oscillatory behavior at early time (right) in toroidal Alfvénic time units, for asymmetric and symmetric geometry.

Figure 3 . 10 :

 310 Figure 3.10: Streamlines colored with toroidal velocity (u T ) for M = 23 (left) and M = 226 (right) for the symmetric torus.

Figure 3 . 11 :

 311 Figure 3.11: Streamlines colored with toroidal velocity (u T ) for M = 23 (left) and M = 226 (right) for the asymmetric torus.

Figure 3 . 12 :Figure 3 . 13 :

 312313 Figure 3.12: The ratio of the mean-square toroidal velocity to the total mean-square u 2 T / u 2 as a function of M.

Figure 3 . 14 :Figure 3 . 15 :

 314315 Figure 3.14: Azimuthally averaged flow visualizations: toroidal velocity u T and poloidal stream function contours (solid line positive, dotted line negative contours) for M = 23 (a), M = 1131 (b) and M = 4524 (c). (d) Toroidal velocity profiles along a vertical cut. The position of these cuts is indicated in (a), (b) and (c) by a dotted vertical line.

1 • 10 - 7 1 • 10 - 6 1 • 10 Figure 3 . 16 :

 10710610316 Figure 3.16: Kinetic energy as a function of M for the asymmetric and symmetric cross sections.

Figure 3 . 17 :

 317 Figure 3.17: Root mean square value of the Lorentz force as a function of the viscous Lundquist number.

Figure 3 . 18 :Figure 3 . 19 :

 318319 Figure 3.18: Azimuthally averaged vector norm of the Lorentz force for M = 23 (a), M = 1131 (b) and M = 4524 (c).

Figure 3 . 20 :Figure 3 . 21 :

 320321 Figure 3.20: Volume-averaged absolute value of the cosine of the angle between the poloidal current density (j pol ) and the poloidal magnetic field (B pol ).

Figure 3 . 22 :

 322 Figure 3.22: RMS value of the non azimuthally symmetric magnetic fluctuations, normalized by the total root mean square perturbed magnetic field.

Figure 3 . 23 :Figure 3 . 24 :

 323324 Figure 3.23: Azimuthally averaged square velocity fluctuations around the azimuthal mean value for M = 23 (a), M = 1131 (b) and M = 4524 (c).

  3.23 and 3.24) show the velocity fluctuations mainly concentrated at the boundaries. These are the areas where the velocity is peaked (see e.g. velocity profiles Fig.3.14 (d)) and where the velocity gradients are important. For the magnetic field the fluctua-tions are spread in a larger region, they are more important at the high and low field side of the torus.

Figure 3 . 25 :

 325 Figure 3.25: Left: kinetic energy evolution. Right: a zoom on the early time instants. Time is given in toroidal Alfvénic time units.

Figure 3 . 26 : 1 • 10 - 6 2 • 10 - 6 3 • 10 - 6 4 • 10 - 6 5 • 10 - 6 6 • 10 - 6 7 • 10 Figure 3 . 27 :

 326110610610610610610610327 Figure 3.26: Left: perturbed magnetic energy evolution. Right: a zoom on the early time instants. Time is given in toroidal Alfvénic time units.

Figure 3 . 28 :Figure 3 . 29 :

 328329 Figure 3.28: Current helicity as a function of the safety factor q.

Figure 3 . 30 :

 330 Figure 3.30: The ratio of the mean-square toroidal velocity to the total mean-square u 2 T / u 2 as a function of q. In the inset we show the average over the domain of the absolute value of the cosine of the angle between the velocity field and magnetic field.

Figure 3 . 31 :

 331 Figure 3.31: Normalized toroidal angular momentum L T /L T rms as a function of q.

- 1 • 10 - 3 - 5 • 10 - 4 0 5 • 10 - 4 1 • 10

 11035104510410

Figure 3 . 32 :

 332 Figure 3.32: Azimuthally averaged flow visualizations: toroidal velocity u T and poloidal stream function contours (solid line positive, dotted line negative contours) for q = 5.7 (a), q = 3.8 (b) and q = 2.9 (c). (d) Toroidal velocity profiles along a vertical cut. The position of these cuts is indicated in (a), (b) and (c) by a dotted vertical line.

Figure 3 . 33 :

 333 Figure 3.33: Toroidal velocity profiles along a horizontal cut, at the center of the domain.

5. 48 • 10 - 3 5. 52 • 10 - 3 5. 56 • 10 - 3 5. 60 • 10 - 3 5. 64 • 10 - 3 5Figure 3 . 34 :

 4810352103561036010364103334 Figure 3.34: Square velocity fluctuations normalized by the total square velocity.

Figure 3 . 35 :

 335 Figure 3.35: Square magnetic fluctuations normalized by the total square perturbed magnetic field.

Figure 3 . 36 :

 336 Figure 3.36: Azimuthally averaged flow visualizations: toroidal velocity u T and poloidal stream function contours (solid line positive, dotted line negative contours) for imposed positive B 0 T (a) and negative B 0 T (b). (c) Toroidal velocity profiles along a vertical cut. The position of these cuts is indicated in (a) and (b) by a dotted vertical line.

Figure 4 . 1 :

 41 Figure 4.1: Wire model (top view): (left) initial position of the current-carrying wire (red) surrounded by the solenoid with poloidal current j P and (right) development of the perturbed poloidal current density, j ′ P (when the wire in the center is kinked) with enhancement of the axial magnetic field in the center (+B z ) and decrease at the edge (-B z ).

Figure 4 . 2 :

 42 Figure 4.2: Toroidal circular geometry.

Θ = 1 . 6 ,

 16 M = 658 Θ = 1.5, M = 658 Θ = 2.0, M = 329

Figure 4 . 3 :

 43 Figure 4.3: Pinch ratio evolution as a function of the poloidal Alfvén time.

Figure 4 . 4 :

 44 Figure 4.4: Reversal parameter evolution as a function of the poloidal Alfvén time.

Θ = 1 . 6 ,

 16 M = 658 Θ = 1.5, M = 658 Θ = 2.0, M = 329

Figure 4 . 5 :

 45 Figure 4.5: Perturbed toroidal component of the magnetic field as a function of the poloidal Alfvén time.

Figure 4 . 6 :

 46 Figure 4.6: Velocity streamlines colored with the toroidal velocity for M = 329, Θ = 0.32 (left) and M = 658, Θ = 1.45 (right).

Figure 4 . 7 :

 47 Figure 4.7: Ratio of the kinetic energy of the dominant toroidal modes over the total kinetic energy for the torus geometry, M = 329 (left) and M = 658 (right) as a function of Θ.

Figure 4 . 8 :

 48 Figure 4.8: Toroidal velocity isocontours +0.01 (blue), -0.01 (orange) for M = 658, Θ = 1.45 (left) and Θ = 1.57 (right).

Figure 4 . 9 :

 49 Figure 4.9: Toroidal perturbed magnetic isocontour colored by the toroidal velocity (top view) for M = 658, Θ = 1.45 (left) and Θ = 1.57 (right).

2 • 10 - 5 4 • 10 - 5 6 • 10 - 5 8 2 nFigure 4 . 10 : 2 • 10 - 5 4 • 10 - 5 6 • 10 2 nFigure 4 . 11 :

 10510510524102105105102411 Figure 4.10: Toroidal spectrum of the square velocity (left) and toroidal spectrum of the square toroidal magnetic field component (right) for M = 658 and Θ = 1.45.

Figure 4 . 12 :

 412 Figure 4.12: Azimuthally averaged toroidal (colour) and poloidal (vectors) components of the velocity field for (a) Θ = 0.32 (here vector norm ×10), (b) Θ = 1.44 and (c) Θ = 2.03, for M = 329.

Figure 4 . 13 :

 413 Figure 4.13: Azimuthally averaged perturbed toroidal magnetic field for Θ = 0.32 (a), Θ = 1.44 (b) and Θ = 2.03 (c), for M = 329 at the saturated state.

ΘFigure 4 . 14 : 1 • 10 - 2 2 • 10 - 2 3 • 10 - 2 4 • 10 Θ = 1 . 6 ,

 41411021021021016 Figure 4.14: Kinetic energy evolution as a function of the poloidal Alfvén time.

Figure 4 . 15 :

 415 Figure 4.15: Evolution of the magnetic energy of the perturbation as a function of the poloidal Alfvén time.

Figure 4 . 16 :

 416 Figure 4.16: Kinetic energy as a function of the pinch ratio.

Figure 4 . 17 :

 417 Figure 4.17: Fluctuating magnetic energy as a function of the pinch ratio.

Figure 4 . 18 :

 418 Figure 4.18: RMS Lorentz force as a function of the pinch ratio.

Figure 4 . 19 :

 419 Figure 4.19: Volume averaged current helicity as a function of the pinch ratio.

Figure 4 . 20 :

 420 Figure 4.20: Azimuthally averaged vector norm of the Lorentz force for Θ = 0.32 (a), Θ = 1.44 (b) and Θ = 2.03 (c), for M = 329 at the saturated state.

Figure 4 . 23 :

 423 Figure 4.23: Temporal standard deviation of the kinetic (left) and perturbed magnetic energy (right) at the saturated state for the torus geometry.

Figure 4 . 25 :

 425 Figure 4.25: Pinch ratio evolution for the cylindrical geometry as a function of the poloidal Alfvén time.

Figure 4 . 26 :

 426 Figure 4.26: Reversal parameter evolution for the cylindrical geometry as a function of the poloidal Alfvén time.

Figure 4 . 27 :

 427 Figure 4.27: Perturbed toroidal component of the magnetic field as a function of the poloidal Alfvén time.

Figure 4 . 28 :

 428 Figure 4.28: Ratio of the kinetic energy of the dominant modes over the total kinetic energy for the cylindrical geometry as a function of Θ, M = 634 (left) and M = 1902 (right).

Figure 4 . 29 :

 429 Figure 4.29: Axial velocity isocontours +0.011 (blue) and -0.011 (orange) for Θ = 2.1 (left) and Θ = 2.8 (right), for both cases M = 634.

Figure 4 . 30 :

 430 Figure 4.30: Axial perturbed magnetic field isocontour colored by the axial velocity for M = 634, Θ = 2.1 (left) and Θ = 2.8 (right).

0 1 • 10 - 5 2 • 10 - 5 3 • 10 - 5 4 • 10 - 5 5 • 10 - 5 6Figure 4 . 31 :

 1105105105105105431 Figure 4.31: Axial spectra of the square velocity field (left) and of the axial component of the perturbed magnetic field (right) for M = 634 and Θ = 2.1.

Figure 4 . 34 :

 434 Figure 4.34: Evolution of the kinetic energy (left) and of the perturbed magnetic energy (right) for the cylindrical geometry as a function of the poloidal Alfvén time, M = 1902 and Pr = 2.

Figure 4 . 35 :

 435 Figure 4.35: Evolution of the kinetic energy (left) and of the perturbed magnetic energy (right) for the cylindrical geometry as a function of the pinch ratio.

Figure 4 . 36 :

 436 Figure 4.36: Root mean square of the Lorentz force as a function of Θ.

Figure 4 . 37 :

 437 Figure 4.37: Current helicity as a function of Θ.

Figure 4 . 38 :

 438 Figure 4.38: Axially averaged vector norm of the Lorentz force for Θ = 0.90 (a), Θ = 1.88 (b) and Θ = 2.70 (c), for M = 1902 at the saturated state.

76

 76

Figure 4 . 39 :

 439 Figure 4.39: Time evolution of the kinetic energy contained in the first axial modes for M = 634 and Pr = 3.

88

 88

38

 38

Figure 4 . 40 :

 440 Figure 4.40: Time evolution of the kinetic energy contained in the first axial modes for M = 1902 and Pr = 2.

Figure 4 .

 4 Figure 4.41: Torus: time and azimuthally averaged toroidal velocity for M = 658, Pr = 3 and Θ = 1.57.

Figure 4 .

 4 Figure 4.42: Cylinder: time and axially averaged axial velocity for M = 1902, Pr = 2 and Θ = 1.88.

Θ

  Torus M = 658 Torus M = 329 Cyl. M = 1902 Cyl. M = 634

Figure 4 . 43 :

 443 Figure 4.43: Ratio of the zero (toroidal or axial) mode kinetic energy over the total kinetic energy for the toroidal and cylindrical geometries as a function of Θ.

0 1 • 10 - 3 2 • 10 - 3 3 • 10 - 3 4 • 10 - 3 5 • 10 - 3 6 • 10 - 3 7 • 10 - 3 8Θ = 1 . 9 ,

 110310310310310310310319 M = 1902 Θ = 2.1, M = 634 Θ = 2.8, M = 634

Figure 4 . 44 :

 444 Figure 4.44: Kinetic energy evolution at the early times for the toroidal (left) and the cylindrical geometry (right).

  Cyl. M = 1902 Cyl. M = 634 Torus M = 658 Torus M = 329

Figure 4 . 45 :

 445 Figure 4.45: Field reversal parameter F as a function of the pinch parameter Θ for different viscous Lundquist numbers, Pr = 3 for the toroidal simulations and for the cylindrical case with M = 634 and Pr = 2 for the cylindrical calculation with M = 1902. Also experimental data of three RFP devices.

90 -

 90

Figure 4 . 46 :

 446 Figure 4.46: Field reversal parameter F as a function of the pinch parameter Θ for different viscous Lundquist numbers, Pr = 3 for the toroidal simulations and for the cylindrical case with M = 634 and Pr = 2 for the cylindrical calculation with M = 1902.

Figure A. 1

 1 Figure A.1 shows the exact profiles of the solution of the Poisson equation in [-1,0[ and the penalization term in [0,1]. The solution of the penalized problem should converge to these profiles when η → 0.

Figure A. 1 :Figure A. 2 :

 12 Figure A.1: Exact solution with the different penalization methods. With the regularization method, the solution is C 1 on the whole domain.

Figure A. 3 :

 3 Figure A.3: First (top) and second derivative (bottom) of the solution of the penalized problem for the classical (left) and regularized penalization method (right).

Figure A. 4 :

 4 Figure A.4: Penalization errors of u η as a function of η for classical (left) and regularized penalization method (right).

Figure A. 5 :

 5 Figure A.5: Penalization errors of u ′ η (top) and u ′′ η (bottom) as a function of η for classical (left) and regularized penalization method (right).

  B.1). The equation to be solved is the following, ∂ θ (x,t) ∂t = ∇ 2 θ (x,t), (B.6) with non homogeneous Neumann boundary conditions at x = a and x = b,

Figure B. 1 :

 1 Figure B.1: One dimensional diffusion problem scheme for positive heat fluxes at the boundaries.

= 1 • 10 - 2 η = 4 • 10 - 3 η = 1 • 10 - 3 η = 4 • 10 - 4 Figure B. 2 :

 11024103110341042 Figure B.2: Time evolution of the L 2 error for the gradient of θ , between the numerical results and the exact solution, for different resolutions, N (left) and penalization parameters, η (right).

  P ⊥ the Riesz projection-operator. The implicit treatment for the penalization term at the first order implies ûn+1in the last term the χ is convoluted with ûn+1 makes it non-trivial to write the implicit formulation. The following formulation is proposed[START_REF] Jause-Labert | Numerical validation of the volume penalization method in three-dimensional pseudo-spectral simulations[END_REF]:ûn+1 = P ⊥ F   u n + ∆t F -1 P ⊥ ξ n 1 + χ n+1 η ∆t   . (C.4) C CORRESPONDING EXPLICIT SCHEME OF THE SEMI-IMPLICIT PENALIZATION METHODWe know that the χ function takes only the values 0 and 1 depending if we are in the fluid or in the solid region respectively, χ(x,t) = 0 for x ∈ Ω f , the fluid domain 1 for x ∈ Ω s , the solid domain.(C.5)An equivalent form of the Eq. (C.4) can be written as followsûn+1 = P ⊥ F u n + ∆t F -1 P ⊥ ξ n 1 -χ n+1 ∆t η + ∆t . (C.6)We recover the Navier-Stokes equation in the fluid domain and the implicit penalized equation (C.4) in the solid region. Using the following relationsP ⊥ P ⊥ â = P ⊥ â P ⊥ û = û, (C.7)the differential form of Eq. (C.6) can be written∂ t û = P ⊥ P ⊥ ξχu η + ∆t -∆t η + ∆t F χ F -1 P ⊥ ξ . (C.8)Here we can identify two extreme cases. If ∆t << η the equation (C.8) converges towards the equation (C.2), we recover the classical penalized Navier-Stokes equation. On the other hand if η << ∆t we recover the following equation:∂ t û = P ⊥ F (1 -χ) F -1 P ⊥ ξχu ∆t . (C.9)This equation is very close to the classical penalized Navier-Stokes equation. Here the time step (∆t) replaces the penalization parameter (η) and the nonlinear term vanishes in the solid region.

1 RFigure D. 1 :

 11 Figure D.1: Poloidal magnetic field lines (χ = RA 0 T = constant) for the different cross sections: (a) symmetric, (b) asymmetric. In (c) we show the field lines for the asymmetric geometry used in Ref. [73].

  Fig. D.1 ((a) and (b)) respectively for the considered symmetric and asymmetric geometries (Fig. 3.1).

Figure E. 1 :

 1 Figure E.1: Toroidal circular geometry.

5 )

 5 Calculating this integral we find Ẽ * = 2.5 • 10 -9 . Then we fix the resistive Lundquist number (S) to have the toroidal current I = 3.1 MA. The value of S is determined using the conductivity σ and considering Ohm's law:I = A σ E section , (E.6)The formula used to find the resistive Lundquist number (S) is:S = C a Lσ 4π • 9 • 10 9 c 2 . (E.7)with L the diameter of the cross section (see Fig.E.1) and c the speed of light. The coefficient 4π • 9 • 10 9 /c 2 allows us to pass from SI units to cgs units. This way we can compare with the value in the literature[START_REF] Kamp | Toroidal steady states in visco-resistive magnetohydrodynamics[END_REF]. Now we can calculate the dimensionless poloidal magnetic field (B * 0 pol ) induced by the dimensionless toroidal electric field (E * 0 ). First we consider the dimensional Ampère's law:∇ × B = µj. (E.8)The second equation needed is the relation between a dimensionless and a dimensional electric field,E * = E C a B JET , (E.9)with E * being the dimensionless field. We replace the current density in Eq. (E.8) by Eq. (E.9) taking into account Ohm's law (j = σ E),∇ × B = µσ B JET C a E * . (E.10)To make dimensionless this last equation we divide it by B JET /L,∇ × B * = C a Lσ µ S E * (E.11)Then we can introduce the dimensionless irrotational electric field profile presented in Eq. (E.1),∇ × B * 0 pol = SE * 0 . (E.12)We consider a constant resistivity hence the current density profile is proportional to the electric field profile. The relation between the dimensionless imposed current density (j 0 ) and electric field (E 0 ) is the following,j * 0 = SE * 0 . (E.13)Finally the dimensionless Ampère's law can be written,∇ × B * 0 pol = j * 0 . (E.14)This equation for B * 0 pol can be solved numerically (see Appendix D) considering the solenoidal constraint, and imposing the normal component of the poloidal magnetic field vanishing at the wall as boundary condition.The resistive Lundquist number is calculated for the circular cross section, S = C a Lσ µ ≈ 1.19 • 10 8 (E.15)

1 . 2 •S 1 . 19 •M 9 . 18 •

 12119918 So we have 1.25 • 10 11 m 2 s -1 > ν > 1.25 • 10 -1 m 2 s -1 [51]. For the circular cross section the viscous Lundquist number is, M ∈ [9.18 • 10 -5 , 9.18 • 10 7 ].(E.16)If we use the arbitrary value ν = 1m 2 s -1 , we find M = 1.51 • 10 7 . Which is at present impossible to take into account numerically.To summarize, the parameters of JET tokamak are presented in Tab. E.1. Note that other characteristic quantities of the spherical tokamak NSTX can be found in the Appendix A of the article by Ferraro and Jardin[START_REF] Ferraro | Calculations of two-fluid magnetohydrodynamic axisymmetric steady-states[END_REF]. 10 -7 kg/m3 or 3.6 • 10 19 m -3 Vacuum Permeability µ 4π • 10 -7 Vs/(Am) Alfvén velocity C a 7.2 • 10 6 m/s B JET / √ ρ µ Reference electric field 2.019 • 10 7 V/m C a B JET Toroidal current I 3.1 • 10 6 A Dimensionless applied electric Ẽ * (found in the formula:E * 0 (R, Z) = Ẽ * R 0 R e T ) density field j * (found in the formula: j * 0 (R, Z) = j * R 0 R e T )0.298 j * = S Ẽ * Plasma conductivity σ allowing to have a current of 3.1 MA 8.27 • 10 6 Siemens / m (circular cross section) σ = I AE section 10 8 (circular cross section) S = C a Lσ µ 10 -5 to 9.18 • 10 7 (circular cross section) M = C a L ν Table E.1: JET tokamak parameters.

F. 1 .

 1 The hereby generated current induces a Lorentz force in the azimuthal direction. It is this Lorentz force which is supposed to induce a large scale motion.For the preliminary results presented here the viscous and resistive Lundquist numbers § This investigation is carried out in collaboration with Francesco Palermo, Fabien S. Godeferd and Nicolas Plihon.

  are M = S = 10, corresponding to unity magnetic Prandtl number. The imposed electric field, magnetic field and current density are given the value E 0 = 5 and B 0 = 0.5 and j 0 = 50, respectively. Preliminary results are shown in Fig. F.1 (right) in which it is shown how the plasma is set into movement in the horizontal plane containing the electrodes. Three-dimensional flow visualizations are shown in Fig. F.2, where it is observed that a dipolar velocity structure is created in the plasma.

Figure F. 1 :

 1 Figure F.1: Studied geometry and imposed fields (left). Two dimensional velocity streamlines colored with the azimuthal velocity in a horizontal plane at the axial coordinate where the electric field is introduced (right).

Figure F. 2 :

 2 Figure F.2: Three dimensional velocity streamlines colored with the azimuthal velocity for early time simulations at low magnetic and kinetic Reynolds numbers. Two views: from a direction orthogonal to the imposed radial electric field (left), and along the imposed radial electric field (right).

Table 2 .

 2 1: Parallel performances of the numerical code.

		100
	Time per iteration per processor (s)	1 10	p -1	HD -128 3 HD -256 3 MHD -256 3 HD -512 3 MHD -512 3 MHD -128 3
		0.1	10	100
			Processors

Table 2 .

 2

	Time step Theory Numerical Error Theory Numerical Error
		ω θ	2	ω θ	2	ω θ	2	ω z	2	ω z	2	ω z	2
	1 • 10 -2 1 • 10 -3 5 • 10 -4 1 • 10 -4	14.68	12.68 13.99 14.13 14.28	14 % 5 % 4 % 3 %	5.78	4.99 5.51 5.57 5.62	14 % 5 % 4 % 3 %

2: Comparison of analytical (see Eq. (2.37)) and numerical calculation of the decay rate for the Ohmic diffusion in an infinite cylinder.

Table 2 .

 2

3: Maximum velocity components comparison and relative errors for Re = 120 in a periodic cylinder.

Table 3 .

 3 1: Corresponding viscous Lundquist number for each safety factor value.

† This chapter is an adapted version of the manuscript[START_REF] Morales | Simulation of confined magnetohydrodynamic flows using a pseudo-spectral method with volume penalization[END_REF], reporting on the investigation carried out in collaboration with Matthieu Leroy and Kai Schneider in Marseille.§ Chap. III is an adaptation of a manuscript[START_REF] Morales | Magnetohydrodynamically generated velocities in confined plasma[END_REF], which was written in collaboration with Kai Schneider and David C. Montgomery.

TWO-DIMENSIONAL VALIDATION

THREE-DIMENSIONAL VALIDATION

† An additional effect limiting the reversal of the toroidal magnetic field can be the periodicity condition imposed by our numerical method and we are currently investigating this. The outcome of the test and possible solution to the problem will be presented elsewhere.

10 -4 10 -3 0 500 1000 1500

10 -4 10 -3 0 500 1000 1500 2000

10 -4 10 -3 0 500 1000 1500 2000 behavior as presented previously: at early times we have strong oscillations followed by a relatively calm stationary state. For the rest of the cases, with large Θ, the oscillations that started in the early times are not damped and a strongly fluctuating steady state develops.

10 -4

10 -3 0 500 1000 1500 2000 2500

10 -4

10 -3 0 500 1000 1500 2000 2500

10 -4 10 -3 0 500 1000 1500 2000

10 -4 10 -3 0 500 1000 1500 2000 Notice also the change in the dominant toroidal mode: after the instability the n = 3 mode holds the highest energy value but for the last oscillating three cases the dominant mode is n = 4.

A THEORETICAL ANALYSIS OF THE HERMITE REGULARIZATION 1 ANALYSIS OF A MODEL PROBLEM

In this section a theoretical analysis of the one-dimensional penalized Poisson equation imposing Dirichlet boundary condition is conducted, and the order of convergence yielded by the Hermite regularization method will be characterized and compared to the classical penalization.

The unpenalized problem in [-1;0] reads:

The penalized problem is solved in the interval [-1;1]:

At the fluid-solid interface we impose the continuity condition

For the classical penalization, we set u wall = 1 2 and thus u wall (1) = 1 2 . The solutions are

* -e 2 At the steady state (when the error between the numerical results and the exact solution is constant in time) we can evaluate the different error norm and compute their convergence as a function of the penalization parameter and the resolution. The