Numerical simulation of the Tayler instability in liquid metals

The electrical current through an incompressible, viscous and resistive liquid conductor produces an azimuthal magnetic field that becomes unstable when the corresponding Hartmann number exceeds a critical value in the order of 20. This Tayler instability, which is not only discussed as a key ingredient of a non-linear stellar dynamo model (Tayler-Spruit dynamo), but also as a limiting factor for the maximum size of large liquid metal batteries, was recently observed experimentally in a column of a liquid metal (Seilmayer et al., Phys. Rev. Lett. 108, 244501, 2012}. On the basis of an integro-differential equation approach, we have developed a fully three-dimensional numerical code, and have utilized it for the simulation of the Tayler instability at typical viscosities and resistivities of liquid metals. The resulting growth rates are in good agreement with the experimental data. We illustrate the capabilities of the code for the detailed simulation of liquid metal battery problems in realistic geometries.


Introduction
Current driven instabilities have been known in plasma physics for many decades [2]. A paradigm for their occurrence is the z-pinch, a cylindrical plasma column with an electrical current in direction of the cylinder axis that produces an azimuthal magnetic field. Recently, Bergerson et al have evidenced the onset of kink-type instabilities in a line-tied screw pinch when the ratio of axial to azimuthal magnetic field (the safety parameter) drops below a certain critical value [3].
For a purely azimuthal magnetic field, the onset of current driven instabilities depends basically on the detailed radial dependence B ϕ (r). For the kink-type (m = 1) instability, the relevant criterion is which had been derived by Vandakurov [4] and Tayler [5]. Strictly speaking, this condition holds only without taking into account the stabilizing role of rotation, viscosity, resistivity, or density stratification. Focusing on the latter, the term Tayler instability had been coined by Spruit [6] to describe a situation in which the azimuthal magnetic field becomes strong enough to act against the stable stratification in a star. This instability is quite remarkable since it could provide the second ingredient, in addition to differential rotation, for an alternative type of stellar dynamos, which is called now the Tayler-Spruit dynamo. Although this sort of non-linear dynamo is presently controversially discussed and might fail to work under realistic conditions [7], the underlying kink instability could have important astrophysical consequences, in particular for the extreme spin down of the cores of white dwarfs [8], for chemical mixing in stars [9], or for the occurrence of helical structures in cosmic jets [10]. Resistivity and viscosity have a similar stabilizing effect on the kink instability as density stratification. By slightly stretching the original semantics [6], we will use here the terminus Tayler instability (or TI) for this viscous and resistive case, too, keeping in mind that in either case the azimuthal magnetic field has to exceed a certain critical value in order to become unstable.
In plasma physics the effect of viscosity and resistivity has been studied for various boundary conditions and profiles of the electrical current [11,12]. Apart from details, it was shown that the onset of TI occurs only when the product of the Lundquist number S = µ 0 σRv A and the so-called Alfvén-Reynolds number A = ν −1 Rv A exceeds a value in the order of 10 3 . Here, µ 0 is the magnetic permeability constant, σ and ν are the conductivity and the viscosity of the fluid, R is the radius of the pinch, and v A := B(µ 0 ρ) −1/2 is the Alfvén speed which is proportional to the magnetic field B (ρ is the density of the fluid).
Going over from plasma physics, with all its peculiarities due to compressibility, anisotropic viscosities and conductivities, as well as complicated boundary conditions, to the paradigmatic case of an incompressible, viscous and resistive cylindrical fluid with a homogeneous current distribution, Rüdiger et al were able to show that the key parameter is the Hartmann number, Ha = BR(σ/ρν) 1/2 , which should exceed a value of approximately 25 for the TI to set in when we take B = B ϕ (R) [13]. This is completely consistent with the former results in plasma physics because Ha 2 = SA. With typical viscosities and conductivities of liquid metals the critical electrical current is in the order of 10 3 A. Using the eutectic alloy GaInSn, we had experimentally confirmed the onset of TI at a critical current of approximately 2.7 kA, as well as the numerically predicted increase of this threshold when insulating inner cylinders were inserted into the liquid metal.
Interestingly, currents in this order would indeed be relevant for large-scale liquid metal batteries which are presently strongly discussed as cheap means for the storage of highly intermittent renewable energies [14,15]. Smaller versions of these devices were already studied in the sixties [16,17]. Basically, liquid metal batteries consist of a selfassembling stratified system of a heavy liquid metalloid (e.g. Bi, Sb) at the bottom, an appropriate molten salt mixture as electrolyte in the middle, and a light alkaline or earth alkaline metal (e.g. Na, Mg) at the top (see figure 1). Choosing Na and Bi as an example, Na will loose one electron during the discharge process, turning into Na + .
This ion diffuses through the electrolyte into the lower Bi layer where it is reduced and alloys with Bi to NaBi. Thus, the lower part of the battery will increase in size during discharge and the volume of the upper part will decrease correspondingly. . The electrolyte works as ion conductor and separates the two liquid metals. A movement of the fluid may wipe the electrolyte and lead to an internal short-circuit and thereby a battery failure (right). The TI will induce one or several convection cells in the conducting liquid metal layers.
While small versions of such batteries have already been tested [15], the occurrence of the TI could represent a serious problem for the integrity of the stratification in larger batteries [18]. This is all the more an issue as the highly resistive electrolyte should be chosen as thin as possible in order to maintain a reasonable cell voltage. In [18] we had proposed a simple trick to suppress the TI in liquid metal batteries by just returning the battery current through a bore in the middle. By the resulting change of the radial dependence of B ϕ (r) it is possible to prevent the (ideal) condition (1) for the onset of the TI.
In spite of such attempts to suppress the TI, it is worthwhile to study in detail the flow structure that would arise from it, and the resulting consequences on the stratification of the three layer system. A first step in this direction is the determination of the final velocity field in the saturated state of the TI. The code utilized up to now for the simulation of the liquid metal experiment was well capable of determining the critical current and the growth rates for the onset of the TI [19], but the resulting scale of the velocity in the saturated state was less secure. This has to do with the fact that liquid metals are characterized by a very small magnetic Prandtl number P m := νµ 0 σ ∼ 10 −6 ...10 −5 . The usual numerical schemes for the simulation of TI, which solve the Navier-Stokes equation for the velocity and the induction equation for the magnetic field, are working typically only for values of P m down to 10 −2 . When applying the code at these much too high P m, and scaling the resulting velocity level by ν/R, one arrives at velocity scales in the order of mm/s. However, an extrapolation over 4 orders of magnitude is somehow risky, and needs definitely a justification by codes that can cope with realistic values of P m.
In this paper, we circumvent the P m limitations of those codes that rely on the differential equation approach. We do this by replacing the solution of the induction equation for the magnetic field by applying the so-called quasistatic approximation [20]. This approximation means that we skip the explicit time dependence of the magnetic field by computing the electrostatic potential by a Poisson equation, and deriving then the electric current density. However, in contrast to many other applications in which this procedure is sufficient, in our case we cannot avoid to compute the induced magnetic field, too. The reason for that is the existence of an externally applied current in the fluid. Computing the Lorentz force term it turns out that the product of the applied current times the induced field is of the same order as the product of the magnetic field (due to the applied current) times the induced current. Here, we compute the induced current density from the induced magnetic field by means of Biot-Savart's law. This way we arrive at an integro-differential equation approach, as it had already been used by Meir and Schmidt [21].
In the following, we will describe the mathematical basis of the integro-differential equation approach in the quasistatic approximation. Then we will present the developed numerical model that utilizes the open source CFD library OpenFOAM R [22], supplemented by an MPI-parallelized implementation of Biot-Savart's law. We will discuss, in particular, the convergence properties of this numerical scheme when applied to a TI problem in cuboid geometry. Based on this, we will study the effect of a varying geometric aspect ratio on the spatial structure of the TI and on the critical current. We will also discuss the effect of spontaneous chiral symmetry breaking as it was recently discussed by Gellert et al [23] and Bonanno et al [24]. Then, we will apply our model to the cylindrical TI experiment [1], for which we will show a surprising agreement between measured and computed growth rates. The paper closes with a discussion of the results, with an outlook towards more detailed simulations of liquid metal batteries, and with prospects for a larger TI-related experiment.

Mathematical Model
The initial point for describing fluid dynamics in a liquid metal is the Navier-Stokes with v denoting the velocity, p the pressure, ρ the density, ν the kinematic viscosity, and f the body force density. For incompressible fluids the continuity equation ∇ · v = 0 has to be taken into account. The trigger of the TI is the Lorentz force, with J meaning the current density. Note that the magnetic field B consists of two parts: the static field B 0 , generated by the applied current I (or the corresponding current density J 0 ), and the induced magnetic field b, generated by the motion of the electrically conducing fluid. Although in our problem b is rather small, it must not be neglected in the expression for the Lorentz force because it is multiplied with the large current density J 0 . The usual way to simulate the magnetic field evolution is by solving the induction equationḂ The solution of this equation requires suitable boundary conditions for the magnetic field, which can be implemented either by solving a Laplace equation in the exterior of the fluid [25,26], or by equivalent boundary element methods [27,28,29].
In the following, we choose a second option, i.e. we compute the magnetic field using Biot-Savart's law which is the inversion of Ampère's law ∇×B = µ 0 J . Here r means the point coordinate where B has to be computed, and r ′ is the integration coordinate, which runs through the whole domain where J exists. This way, the problem is shifted from the determination of the magnetic field B to the determination of the current density J . The approach to use Biot-Savart's law was already proposed and utilized by Meir et al [21], and will be adopted to our case.
The flow chart of the numerical model is shown in figure 2. Assuming the current I through our cylindrical vessel as given (e.g. the battery charging current), we compute the corresponding current density J 0 , and the associated static magnetic field B 0 . For an infinitely long vertical cylinder we would obtain in Cartesian coordinates x, y, z. For real geometries with external leads to the system, B 0 could still be be computed by Biot-Savart's law.
In the main loop of our numerical scheme, the Navier-Stokes equation (2) is solved firstly, followed by a velocity corrector step to ensure continuity (∇· v = 0). Then we have to find the electric current density. Presupposing the magnetic Reynolds number Rm = µ 0 σRv on the basis of the the TI-triggered velocity scale v to be small, we can invoke the quasistatic approximation [20]. This means that we express the electric field by the gradient of an electrostatic potential, E = −∇Φ. Applying the divergence operator to Ohm's law in moving conductors, J = σ(−∇Φ + v × B), and demanding charge conservation, ∇ · J = 0, we arrive at the Poisson equation for the perturbed electric potential ϕ = Φ − J 0 z/σ, where z is the coordinate along the battery axis pointing into the direction of the applied potential difference. By subtracting the electric potential caused by the battery charging current I we can easily set the boundary conditions at the electrodes to ϕ = 0. Since no current is allowed to flow through the insulating rim of the cylindrical vessel, we assume here ∇ϕ = 0.
Having derived the potential in the fluid by solving the Poisson equation, we easily recover the current density induced by the fluid motion as and then the induced magnetic field by using Biot-Savart's law (5). The total current density is In a last step, the Lorentz force, to be implemented into the Navier-Stokes equation, is calculated according to It should be noticed again that the (small) induced magnetic field b cannot be omitted here, since its product with the (large) impressed current J 0 is of the same order as the product of the (small) induced current j with the (large) magnetic field B 0 .
Assuming an implementation of this scheme on a certain grid, one has to ask for the allowed time steps to keep the time evolution of the system stable. For this it is important to consider not only the Courant number based on the hydrodynamic velocity v, but also the Alfvén-Courant number based on the Alfvén velocity v A = B/(µ 0 ρ) 1/2 . The influence of the Alfvén-Courant number on the accuracy of the simulations is discussed in the following section.

Numerical scheme and convergence
In this section we will examine the numerical scheme with respect to its convergence and error characteristics. We will both study variations of the grid spacing as well as of the time steps. We will further consider the case that different grid spacings and time steps are used for the Navier-Stokes equation and for the Biot-Savart law.
In order to utilize the convenience of orthogonal cells for grid refinement, a cuboid with the dimensions of 96 × 96 × 120 mm 3 is used as a reference geometry, see figure 3.
For this case, figure 4 shows a typical evolution of the TI in time, comprising an initial phase were the eigenfield of the TI develops, followed by a long period (over ∼ 600 s) in which this eigenfield grows exponentially, and a final phase in which the TI saturates at a certain velocity level.
In the following we will compare the growth rates, as they can be read off during the exponential growth phase for different grid spacings and time steps, with respect to their relative errors which are derived by dividing the absolute error by the difference of the TI growth rates at 10 kA and 0 kA.

Grid refinement
We start by exploring the effects of different grid sizes of the cubic cells with an edge length between 1 mm and 8 mm. Using the Finite-Volume-Method (FVM) of OpenFOAM R for the simulation, the values from the cell centres have to be interpolated to the cell faces in order to ensure continuity. For this interpolation we employ two different methods, a linear one and a cubic one, and compare the results in figure 5.
As expected, the cubic interpolation provides better results, i.e. the growth rates are typically higher. While it would be numerically very demanding to go below a cell size of 1 mm, we expect the difference between cubic and linear interpolation to be zero in the limit of vanishing cell size. With that assumption we do a Richardson extrapolation to obtain the "true" growth rate of the TI. Taking this value as reference, the relative error of the growth rates is then calculated in dependence on the grid cell size ( figure  6).
Since the solution of the Biot-Savart law is numerically very costly, we have additionally checked the possibility to speed up this part of the simulation. While the two lines in figure 5 result from computations in which the Navier-Stokes equation and the Biot-Savart law are implemented on the same grid, the green circles and the liqht blue squares correspond to cases in which the Biot-Savart law is realized on coarser grids. In contrast to the b-calculation with cells of 8 mm edge length, which does not give good results (see figure 5), a grid cell size of 6 mm already seems to be sufficient. Solving the Navier-Stokes equation on a 3 mm grid, and computing Biot-Savart's law on a 6 mm grid, provides almost the same results as doing both on a 3 mm grid. This is numerically very convenient because Biot-Savart's law needs the square of cell numbers as operations whenever carried out, i.e. a coarse grid makes simulation much faster. The resulting slightly lower growth rates when working with two meshes appear due to the inter-grid mapping. The interpolation flattens the b-field a little bit which reduces the growth rate -but the difference is very small.

Alfvén-Courant number and time steps
As usually, during the time integration of the NSE we have to fulfill the Courant-Friedrichs-Lewy (CFL) condition for the time step. In our special problem, however, information is not only transferred by the flow velocity, but also electromagnetically with the Alfvén velocity v A = B/(µ 0 ρ) 1/2 . We consider this by introducing an "Alfvén-Courant number" and compute the maximum time step on that basis. Due to the high currents at which TI is expected to set in, and the consequently high magnetic field B 0 , this time step has to be very short in spite of the low fluid velocities. Fortunately, it turns out that choosing the Alfvén-Courant number larger than one does not lead to a completely unstable system, as it would be expected for a purely differential equation system. Evidently, the mixed integro-differential equation approach makes the system somewhat more "benign" and less vulnerable to violations of the CFL condition. Of course, the error increases significantly with increasing Alfvén-Courant number as shown in figure 7.  Besides calculating the magnetic field b on a coarser grid than the velocity field, we have also tested to compute it only at every n th timestep. Figure 8 shows that it is well possible to compute the Biot-Savart integrals only every 100 th time step. The benefit in computation time is especially high when solving the NSE and Biot-Savart's law on the same grid, as the latter one takes usually more than 100 times longer than solving all other equations.

Results for cuboid geometry
In this section we will focus on the characteristic eigenmode structure of the TI, on the typical velocities in the saturated regime, and on the critical currents in dependence on various geometry parameters. The simulation are done for a cuboid volume with a base area of 96 × 96 mm 2 , but with varying heights.

Shape, alignment, and chiral symmetry breaking
Basically, the TI appears in form of convection cells that are vertically stacked one over another (see figure 9). The fluid in one cell rotates around some horizontal axis, and the sense of rotation is changing from one cell to the neighboring one. Referring, for the moment, to cylindrical geometry, this cell structure represents a non-axisymmetric flow structure, with an azimuthal wavenumber m = ±1. Since the present simulations are carried out for cuboid geometry, one could imagine, e.g., a preferred alignment of the convection cells with the diagonal of the base area. However, this is generally not the case. In fact, the alignment of the cells depends strongly on the chosen initial conditions. Assuming a small perturbation in the x-direction, the TI will arise around the y-axis ( figure 9). For arbitrary initial velocity, the TI will develop around a random axis.
Most interestingly, a transient helical shape of the TI may result under certain initial conditions (see figure 9, right). This was observed, e.g., with an initial velocity in z-direction and in particular when the height of the simulated geometry did not match an integer multiple of the wave length of the TI. The helical perturbation may develop if, e.g., two convection cells arise simultaneously, but with azimuthal twisted rotation axes. The growth rates of the non-helical and the helical TI state were found to be very similar. Figure 10 shows the velocity and the normalized mean value of the helicity H = v · (∇ × v) of a non-helical (blue) and of a helical (red) state for a 96 × 96 × 360 mm 3 geometry (h = 3.1λ). Reaching a certain velocity level in the saturated regime, the helicity of the helical state reduces significantly and goes finally to zero. Approximetaly at the same time (at ∼2500 s) the saturated state of the formerly non-helical state acquires also some helicity, which goes finally to zero, too. The ultimate convection cells of the two states are not equal (compare figure 9 left and right), and seem to depend on the different initial conditions of the velocity.
This transient appearance of helical TI structures is quite remarkable, since it corresponds to the spontaneous chiral symmetry breaking which was recently observed by Gellert et al [23] and discussed in detail by Bonanno et al [24]. Here we can only evidence its occurrence, while a more comprehensive study of this effect has to be left for future studies. Depending on the initial velocity conditions, the convection cells of the TI align to an arbitrary direction (left: orthogonal, second: diagonal to base cube).
Due to a helical perturbation, a helical TI may grow (3. panel, t = 2 300 s), but almost disappears later (4. panel t = 2 400 s, right t = 2 800 s). The corresponding time graph is shown in figure 10.

Saturation levels of the TI
To investigate the saturation level of the TI we have carried out simulations using different currents for a cuboid geometry of 96×96×120 mm 3 and a grid cell length of 3 mm. When the critical current for the onset of the TI is exceeded, any small disturbance can initiate the instability which then needs some time to establish its eigenmode. Thereafter, this eigenmode grows exponentially, see figure 11. For currents only slightly above the critical one, the instability grows very slowly. In such a case, it may take tens of minutes until the TI reaches its saturated final state. Figure 12 shows the dependence of the final velocities on the applied current. We see that the maximum velocities grow from a few mm/s close to the critical current, to nearly 14 cm/s at 50 kA. Velocities on this scale would clearly be problematic for the stability of the layer stratification in liquid metal batteries.

Critical current in dependence on the aspect ratio
To investigate the dependence of the critical current on the aspect ratio of the fluid volume we consider a cuboid with a base area of 96 × 96 mm 2 which is changed in height from 24 to 480 mm. As a first estimate for the critical wavelength of the TI we can take the value for an infinitely long cylinder of radius R, for which the critical wavenumber k is known to be 2.47 [19]. Embedding this cylinder in the cuboid to be considered here, the resulting wave length would be λ = 122 mm. To obtain the corresponding true spatial wave length for the cuboid geometry, the TI is simulated for a very long cuboid with a height of 1 040 mm. For this case we obtain λ = 115 mm, a value quite close to the corresponding value of the infinitely long cylinder. Figure 13 shows now the critical currents as a function of height. In the limit of very tall cuboids, the critical current converges to a value of approximately 3 kA which is also quite close to the value 2.7 kA for the infinite long cylinder [19]. Shortening the cuboid, a first plateau of the critical current is observed at a height corresponding to the TI wave length of 115 mm. Decreasing the height further, a sharp increase of the critical current is observed, which flattens again to a second plateau when the height comes close to the half wave length (57.5 mm). Shortening the cuboid still further, the critical currents rise extremely steep. With regard to liquid metal batteries this would mean that very flat cells will not be affected by the Tayler instability.
For those heights which are below the half wavelength of the TI (h < λ/2), there is only one convection cell, for slightly taller vessels a second one appears (see figure 13, left). In general, the number of convection cells adapts to the height, and the cell size acquires a maximum if the height is a multiple of the wavelength. If h exceeds 1.5λ, a third cell may appear -but only in case of a helical TI. In case of non-helical TI an even number of convection cells seems to be more stable. Usually, the new convection cells appear at the top and/or bottom. If the height is slightly below an integer multiple of the characteristic wave length, the upper and lower convection cell will appear somehow "compressed" in comparison to the more central ones ( figure 13, right).

Validation of experimental results
In a recent paper [1], we had evidenced the occurrence of the TI in a liquid metal experiment. A cylinder, made of polyoxylmethylene (d = 100 mm, h = 750 mm) was filled with the eutectic GaInSn, being liquid at room temperature. A current of up to 8 kA was injected into the liquid metal through two massive copper electrodes at the top and bottom of the cylinder. In order to avoid any inhomogeneities of the electric current at the interface between Cu and GaInSn, we restrained from inserting ultrasonic transducers for the direct measurement of the velocity perturbations, and relied exclusively on external measurements of the induced magnetic field. The vertical component b z was measured at 11 positions along the vertical axis of the device. With this set-up we were able to identify the growth rates of the short-wavelength perturbations due to the TI, before those gave way to longer wavelength structure that we attributed to the thermal wind resulting from Joule heating in the GaInSn.
For this cylindrical geometry, we simulate the TI and compare the results with the available experimental data. The cylindrical fluid volume is meshed with a grid cell size of about 3 × 3 × 3 mm 3 . As the electrodes are relatively thick and are elongated by long power supply leads, we assume B 0 to be constant over the height of the cylinder.
For an assumed current of 10 kA, the resulting convection cells of the TI are illustrated in figure 14. It should be noticed that the number and shape of these convection cells is quite independent on the initial conditions.
The experimental growth rates are compared with the numerical ones in figure 15. Simulations for the infinitely long cylinder [19] provide the maximum possible growth rates. As we do simulation for the real, i.e. finite cylinder, the resulting growth rates lie slightly below these values, but correspond nicely with the experimental data.
For our geometry the theoretical wave length is λ = 127 mm [19]. This value matches fairly well the wave length obtained by our finite length simulation (121 mm), as inferred from figure 16.  Values are taken at maximal growth rate, i.e. shortly before reaching saturation (black dotted line). The red fit allows to determine the spatial wavelength of λ = 121 mm.

Conclusion and outlook
On the basis of an integro-differential equation approach we have developed a numerical tool that is particularly suited for investigations of current driven magnetohydrodynamic instabilities in liquid metals. Employing the quasistatic approximation, the solution of the induction equation for the magnetic field was substituted by solving a Poisson equation for the induced current density, from which the magnetic field is computed via Biot-Savart's law. This procedure allows to deal with instabilities in conducting fluids characterized by very small magnetic Prandtl numbers. After examining the convergence properties of this numerical scheme, we have simulated the TI in different geometries, including cuboids and cylinders. For the cuboid it was shown in detail how the critical current increases with decreasing aspect ratio, exhibiting two remarkable plateaus at those heights that correspond to one or a half wavelength of the TI mode, respectively. This curve may be directly used for defining a maximum aspect ratio of liquid metal batteries at a given charging current. Depending on the initial conditions our simulations have also shown the transient appearance of helical structures, i.e. the occurrence of chiral symmetry breaking. For cylindrical geometry, we have numerically reproduced the growth rates of the TI that had been measured in a recent liquid metal experiment.
A main result of our simulations is that the saturation level of the TI induced velocity for typical electrical currents are in the order of mm/s until cm/s. This is quite consistent with a formal scaling, by ν/R, of the velocity levels as obtained at much larger magnetic Prandtl numbers [30]. While this good correspondence confirms the interaction parameter Ha 2 /Re to be the governing parameter for the saturation over a wide range of magnetic Prandtl numbers, it makes the physical interpretation of the saturation mechanism a bit more difficult. Actually, for P m close to one, Gellert et al [30] had identified a very plausible saturation mechanism in terms of a radially dependent increase of the TI triggered turbulent resistivity (β effect) that modifies the radial current distribution in such a way as to make the magnetic field distribution just marginally stable. Evidently, this nice picture does not apply any more for small P m, as considered here, since the arising velocity perturbation are much too weak for producing any noticeable turbulent resistivity. In this respect, a simple picture for understanding the saturation of TI at low P m remains elusive.
In a next step we plan to include into the numerical model the effects of thermal convection resulting from the Joule heating due to the strong current, which has turned out important in the experiment. We also develop a multiphase solver which will allow us to simulate three-layer liquid metal batteries in detail.
On the experimental side, in the framework of the DRESDYN project [31] we plan to set-up a large scale liquid sodium experiment for the combined investigation of the magneto-rotational instability, as recently investigated [32], and the TI.