Self-organization of helically forced MHD flow in confined cylindrical geometries

Vortical structure formation plays an important role in magnetohydrodynamic (MHD) flows in a variety of settings and applications, such as the Earthʼs core (Moffatt 1978), metal casting, and fusion devices of the reversed field pinch (RFP) type (Taylor 1974) and of the tokamak type (Wesson 2011). The development of large rotational structures has been addressed before in Bos et al (2008) and Neffaa et al (2008). While some analytic results exist for MHD flows, we must generally rely on numerical simulations in order to understand the behaviour of those fluids in confined geometries. Understanding structures in MHD flows lends insight into the underlying mathematical properties and can help to control MHD flows more efficiently.

In this paper we explore the behaviour of a conducting fluid and the evolution of the flow topology under the influence of external magnetic forcing in two cylindrical geometries. In previous investigations we focused on the structure formation of magnetofluids in two-dimensional (Bos et al 2008, Neffaa et al 2008) and toroidally confined (Morales et al 2012) domains. The case of a cylindrical domain was considered numerically in Shan et al (1991) and Shan and Montgomery (1993b) in the case of non-penetration boundary conditions. This choice was necessary for these simulations because the fields were decomposed into eigenfunctions of the curl operator in order to observe the dynamics of helical structures, and it was argued, via an energy-minimization argument, that the dominant mode of the flow is such a helical basis function. The instability threshold, determined analytically, was confirmed numerically and the evolution of the flow with the magnitude of the forcing was described. Due to the absence of an efficient physical-space to Chandrasekhar–Kendall function-space transform, the nonlinear term was computed directly, with rapidly increasing computational cost. In this paper, we use a pseudospectral method (Canuto et al 1987) combined with a volume penalisation method for imposing boundary conditions in complex geometries (Arquis and Caltagirone 1984, Angot et al 1999). Thus we were able to implement different boundary conditions than used in Shan et al (1991) and to perform simulations in cylindrical geometries with different cross-sections. A related study, which was limited to circular cross-sections and for different parameters focusing on the RFP application can be found in Bonfiglio et al (2013) and Veranda et al (2013). The self-organization of compressible MHD flows into helical structures has also been observed in Horiuchi and Sato (1986) and Zhu (1995), for example, where the mechanism of self-organization is magnetic re-connection.

This paper is organized as follows: in section 2 we give a mathematical description of the problem and present our numerical technique. In section 3, we explore the structure of the velocity field in a circular geometry when forced by an external helical magnetic field, and we show the effect of wrapping number in section 4. Alignment statistics are discussed in section 5. Section 6 shows how the flow topology changes when the cross-section is modified. We give our conclusions and future perspectives in section 7.

2.1. Mathematical description

Let ${\boldsymbol{u}}$ denote the velocity and ${\boldsymbol{B}}$ the magnetic field of the conducting fluid. The time evolution of these quantities using conventional Alfvén normalized variables, is given by

$$\begin{eqnarray}&&\frac{\partial {\boldsymbol{u}} }{\partial t}={\boldsymbol{u}} \times {\boldsymbol{ \omega }} +{\boldsymbol{j}} \times {\boldsymbol{B}} -\nabla P+\nu {{\nabla }^{2}}{\boldsymbol{u}} ,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\frac{\partial {\boldsymbol{B}} }{\partial t}=\nabla \times ({\boldsymbol{u}} \times {\boldsymbol{B}} )+\lambda {{\nabla }^{2}}{\boldsymbol{B}} ,\end{eqnarray} \tag{ 2 }$$

where P is the total pressure, ${\boldsymbol{ \omega }} =\nabla \times {\boldsymbol{u}}$ is the vorticity, and ${\boldsymbol{j}} =\nabla \times {\boldsymbol{B}}$ is the current density. Physical parameters are the kinematic viscosity $\nu =4.5\times {{10}^{-2}}$, and the magnetic diffusivity $\lambda =4.5\times {{10}^{-2}}$. Both the velocity and magnetic field are considered to be incompressible, so $\nabla \cdot {\boldsymbol{u}} =0$ and $\nabla \cdot {\boldsymbol{B}} =0$.

The fluid is confined in a cylinder of length ${{z}_{\ell }}=8$ with its axis oriented along the z-axis. The ends of the cylinder are treated as periodic. The cross-section of the cylinder in the x–y plane is either a circle of radius 1 or an ellipse with major radius 1 along the x-axis and minor radius $1/\sqrt{2}$ along the y-axis corresponding to an eccentricity of $1/\sqrt{2}$. No-slip (and non-penetration) boundary conditions are imposed on the velocity field. A helical forcing is imposed in the boundary conditions for ${\boldsymbol{B}} ,$ where the normal component is set to zero, the z-component is set to B₀, and the poloidal component is set to B_c. This corresponds to the boundary conditions for the magnetic field following a helix along the wall, as illustrated in figure 1. We define the wrapping number q to be the number of rotations achieved by a magnetic streamline following the boundary conditions between z = 0 and $z={{z}_{\ell }}$, and we choose ${{B}_{0}},{{B}_{{\rm c}}}$ pairs which provide integral values of q for the simulations in this paper. It was shown in Shan and Montgomery (1993a, 1993b) that the instability of the current setup is determined by two parameters: the Hartmann number and the pinch-ratio. The Hartmann number was chosen proportional to the axial magnetic field, and the pinch-ratio is proportional to the wrapping number. Our description in terms of B₀ and q is thus equivalent to the previous investigations.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The magnetic field is initialized to a linear profile that matches the boundary conditions. The velocity field is set to an uncorrelated Gaussian solenoidal random field with kinetic energy of the order 10⁻⁷.

2.2. Numerical method

Let us first consider a cylinder with circular cross-section. The kinetic energy of the fluid is initially low, and, after a short time (a few Alfvénic time units, depending on the boundary forcing) the kinetic energy grows exponentially with no change in topology, which is consistent with the linear instability of helical modes as described in Shan et al (1991), before reaching a new steady state. In certain cases (low wrapping numbers / boundary forcing magnitude) this does not occur, and the kinetic energy is essentially zero for the duration of the simulation. Example kinetic energy and kinetic energy dissipation curves are shown in figure 3.

Zoom In Zoom Out Reset image size

We set B₀ = 4.5 and vary B_c so that integral values of wrapping numbers q are achieved. For q = 1 no flow is observed. At q = 2 (${{B}_{{\rm c}}}\approx 7.06$) we observe a helical structure in the velocity field consisting of two pairs of helices, as shown in figure 4. With these forcing parameters, the velocity modes have the same wrapping number as the boundary forcing of the magnetic field, which is also observed with wrapping number 3 (${{B}_{{\rm c}}}\approx 10.6$). At wrapping number 4, the structure begins to break down, and is replaced with a helical dipole at q = 5 (${{B}_{{\rm c}}}\approx 17.6$), with wavenumber 2 in the axial direction, as shown in figure 5. This low-order mode, which resembles the minimal dissipation mode as in Montgomery et al (1989), endures even at q = 20 (see figure 6), at which point the velocity field in the z-direction is disordered and exhibits a full spectrum, as shown in figure 7.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the previous section, we increased forcing by increasing the wrapping number to explore the effect on the topology of the flow. In order to understand the individual roles of forcing and wrapping number, we also perform a number of simulations with wrapping number $q\in \{1,2,3,4\}$ and various forcing amplitudes.

For wrapping number 1, we observe the formation of structures first at forcing amplitude ${{B}_{\parallel }}=\sqrt{B_{0}^{2}+B_{{\rm c}}^{2}}=20$, consisting of two helical structures of opposite-signed axial velocity and axial wavenumber 1, matching the forcing wrapping number. This structure is also present for ${{B}_{\parallel }}\in \{30,40,80,100\}$, and exhibits decreasing regularity as the forcing parameter increased. For wrapping number 2, self-organization occurs at lower forcing values, first appearing at ${{B}_{\parallel }}=10$, at which point we observe three helical pairs, as shown in figure 8. The structure of the flow becomes increasingly complex as ${{B}_{\parallel }}$ increases, and an axial-wavenumber 1 mode appears at ${{B}_{\parallel }}\gt 40$, which persists up to ${{B}_{\parallel }}=100$. At wrapping number 3, we observed axial-wavenumber 1 helical structures at ${{B}_{\parallel }}\in \{30,40,60,80,100\}$, with two helical pairs at low forcing degenerating into less ordered (but still helical) structures as forcing magnitude increased. Finally, simulations with q = 4 show an axial wavenumber 2 helical pair at ${{B}_{\parallel }}=10$, which persists up to ${{B}_{\parallel }}=100$. These results are summarized in figure 9, which may be compared to figure 1 in Shan and Montgomery (1993b), which shows the linearly unstable modes as a function of pinch ratio $\Theta ={{B}_{{\rm c}}}/{{B}_{0}}$ and Hartmann number $H={{B}_{0}}/\sqrt{\eta \lambda }$ in a parameter region different than given here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In figure 10 we show the probability density functions (pdfs) of the cosine of the angles between $({\boldsymbol{u}} ,{\boldsymbol{ \omega }} )$, $({\boldsymbol{j}} ,{\boldsymbol{B}} )$ and $({\boldsymbol{u}} ,{\boldsymbol{B}} )$. These angles are chosen since they locally determine the strength of the Navier–Stokes nonlinearity, Lorentz force and magnetic field induction. Two cases are compared: a laminar one in which the velocity and magnetic field are stationary, and a turbulent multi-mode state. It is shown that already the laminar state shows non-trivial alignment properties. When the flow becomes turbulent the dynamics relax to a state in which in particular the magnetic field aligns more strongly with the velocity field compared to the laminar state. The ${\boldsymbol{u}} -{\boldsymbol{ \omega }}$ alignment becomes weaker and the ${\boldsymbol{j}} -{\boldsymbol{B}}$ alignment does not dramatically increase or decrease in strength.

Zoom In Zoom Out Reset image size

Such rapid local relaxation processes, in particular observed here for the ${\boldsymbol{u}} -{\boldsymbol{B}}$ alignment have been observed previously in hydrodynamic (Pelz et al 1985) and MHD (Servidio et al 2008) turbulence. In particular in the latter work similar alignments as the ones in the present work are investigated in the case of isotropic MHD turbulence. It is observed in their results that the ${\boldsymbol{j}} -{\boldsymbol{B}}$ alignment becomes very strong. Furthermore, due to isotropy all their pdfs showed alignment–antialignment symmetry. This symmetry is lost in the pdfs of ${\boldsymbol{u}} -{\boldsymbol{ \omega }}$ and ${\boldsymbol{j}} -{\boldsymbol{B}}$ of our simulations due to the imposed helical character of the magnetic field. We do not think that the alignment in the turbulent case can be explained in an easy way. However, it was suggested by Kraichnan and Panda (1988) that the ${\boldsymbol{u}} -{\boldsymbol{ \omega }}$ alignment observed in isotropic turbulence might be one representation of a more generic property of turbulent systems to suppress their nonlinearity. The rapid relaxation observed in the present work can also be thought to be an expression of this property.

Let us now consider a cylinder with elliptical cross section defined by ${{x}^{2}}+2{{y}^{2}}=1$. As in the previous section, we keep the wrapping number constant and modify the forcing magnitude. The onset of kinetic energy growth occurs at much stronger forcing than in the circular case, with q = 1, ${{B}_{\parallel }}=60$ our first detection. At this point, the velocity field forms five pairs of opposite-sign helices, as shown in figure 11. As forcing is increased, 6 pairs are observed at ${{B}_{\parallel }}=80$, and 5 pairs again at ${{B}_{\parallel }}=100$, 110, 120, and 130. As forcing amplitude is increased, the velocity field becomes increasingly concentrated at the border.

Zoom In Zoom Out Reset image size

In this article we studied the self-organization of the flow in a conducting fluid confined in a cylinder. The flow is forced by a helical magnetic field imposed via boundary conditions. This phenomenon was first explored by Shan et al (1991) using a purely spectral method for circular cross-sections with non-penetration boundary conditions. The simulations presented here used a large number of grid points, different parameter regimes, and were performed in cylinders with both circular and elliptical cross-sections. The penalisation technique can also be used in more general geometries. In cylinders of circular cross-section, various helical flows develop depending on the value of the wrapping number and forcing amplitude, tending toward a unique low-order mode at large values of the forcing parameters. This mode is present even when the fluid starts to exhibit turbulent dynamics. For the elliptical cross-section, the situation is different: a ring of vortices with poloidal mode number five or six forms and encloses a weak flow area. These vortices follow the topology of the imposed magnetic field. The magnitude of the forcing required to trigger this instability is much larger than for circular cases. For further work, it would be interesting to study the flow for larger ranges of Hartman number and pinch ratio, and to investigate the influence of the eccentricity of the ellipse on the topology of the flows obtained. Indeed, these first results in cylinders with non-circular cross-section reveal the significant influence of the shape on the modes which are excited. This observation might have important implications for the control of RFP fusion reactors, where the shape and frequency of the helically extended modes determine the quality of the plasma confinement. Even though in realistic RFP simulations the strong temperature gradient will have an important influence on the dynamics, already the isothermal visco-resistive MHD description is shown to reproduce a number of key features of the dynamics in Bonfiglio et al (2013).

This work was supported by the contract SiCoMHD (ANR-Blanc 2011-045), with simulations performed on IDRIS under project 91664 and made use of the HPC resources of Aix-Marseille Université financed by the project Equip@Meso (ANR-10-EQPX-29-01) of the program 'Investissements dʼAvenir' supervised by the ANR.

Let $\Omega ={{\Omega }_{{\rm f}}}\cup {{\Omega }_{{\rm s}}}$ be the computational, fluid, and solid domains, respectively, with $\left| {{\Omega }_{{\rm f}}}\cap {{\Omega }_{{\rm s}}} \right|=0$. Given a Dirichlet boundary condition on $\partial {{\Omega }_{{\rm f}}}$ represented by $\{{{{\boldsymbol{v}} }_{{\rm bc}}}({\boldsymbol{x}} )|{\boldsymbol{x}} \in \partial {{\Omega }_{{\rm f}}}\}$ such that

$$\begin{eqnarray}&&{{\int }_{\partial {{\Omega }_{{\rm f}}}}}{{{\boldsymbol{v}} }_{{\rm bc}}}\cdot {\boldsymbol{n}} \;{\rm d}s=0,\end{eqnarray} \tag{ A.1 }$$

where ${\boldsymbol{n}}$ is the normal vector on $\partial {{\Omega }_{{\rm f}}}$. We wish to extend ${{{\boldsymbol{v}} }_{{\rm bc}}}$ to a sufficiently regular function ${\boldsymbol{v}}$ on the entire computational domain $\Omega$ such that $\nabla \cdot {\boldsymbol{v}} =0$. In order to find ${\boldsymbol{v}}$ numerically, we introduce a pseudo-time $\tau$ and time-step the equation

$$\begin{eqnarray}&&\frac{\partial {\boldsymbol{v}} }{\partial \tau }=\kappa {{\nabla }^{2}}{\boldsymbol{v}} -\frac{1}{{{\eta }_{\tau }}}{{\chi }_{\partial {{\Omega }_{{\rm f}}}}}({\boldsymbol{v}} -{{{\boldsymbol{v}} }_{{\rm bc}}}),\end{eqnarray} \tag{ A.2 }$$

where $\kappa$ is a diffusion parameter, ${{\chi }_{\partial {{\Omega }_{{\rm f}}}}}$ is the characteristic function for $\partial {{\Omega }_{{\rm f}}},$ and we call ${{\eta }_{\tau }}$ the pseudo-time penalisation parameter. For the simulations in this paper, we chose $\kappa =100$, and we used an adaptive 1st/2nd order embedded Runge–Kutta time-integrator and a pseudospectral method in space, projecting ${\boldsymbol{v}}$ onto the solenoidal manifold via a Helmholtz decomposition after each time-step. The time-stepping in $\tau$ is stopped when

$$\begin{eqnarray}&&\left\|{{{\boldsymbol{v}} }_{{\rm bc}}}-{\boldsymbol{v}} \right \|{{}_{\infty ,\partial {{\Omega }_{{\rm f}}}}}\lt 0.2\times \sqrt{\eta },\end{eqnarray} \tag{ A.3 }$$

where η is the penalisation parameter used to advance the fluid. This is to say that the maximum error between the penalisation field ${\boldsymbol{v}}$ and the desired value ${{{\boldsymbol{v}} }_{{\rm bc}}}$ on the boundary is less than one fifth the error expected from the penalisation term used in time-stepping the fluid. Since the boundary conditions obey the compatibility condition (equation A.1), the resulting field is consistent with a solenoidal field and should not be affected by projection onto the solenoidal manifold. Moreover, the stopping condition (inequality A.3) is satisfied only when the resulting penalisation field closely matches the boundary conditions.

The stop condition (A.3) implies that the error in the boundary conditions is sub-leading-order, and the use of a Helmholtz decomposition provides a penalisation field which is divergence-free up to machine precision. The choice of $\kappa$ determines the regularity of the field and the number of time-steps required for condition (A.3) to be achieved.