MHD stability of large scale liquid metal batteries

MHD stability of large scale liquid metal batteries

A. Tucs, V. Bojarevics Email address for correspondence: V.Bojarevics@gre.ac.uk    K. Pericleous
Abstract

The aim of this paper is to develop a stability theory and a numerical model for the three density-stratified electrically conductive liquid layers. Using regular perturbation methods to reduce the full 3d problem to the shallow layer model, the coupled wave and electric current equations are derived. The problem set-up allows the weakly non-linear velocity field action and an arbitrary vertical magnetic field. Further linearisation of the coupled equations is used for the linear stability analysis in the case of uniform vertical magnetic field. New analytical stability criteria accounting for the viscous damping are derived for particular cases of practical interest and compared to the numerical solutions for variety of materials used in the batteries. The new criteria are equally applicable to the aluminium electrolysis cell MHD stability estimates.

Key words: Interfacial waves, MHD interaction, linear stability, liquid metal battery, shallow layer approximation, aluminium electrolysis cell

 

1 Introduction

Liquid Metal Batteries (LMBs) have many important characteristics for efficient practical use in combination with renewable energy sources on a national energy grid scale in the future. A relatively high voltage efficiency at high current densities of this storage technology are due to liquid-liquid electrode-electrolyte interfaces that enable high speed charge transfer, high total current capability, low ohmic losses, as well as rapid mass transport of reactants and products to and from the electrode-electrolyte interfaces by means of liquid-state diffusion (Kim et al. 2013).

The liquid state of the main components necessitates consideration of the fluid dynamics in LMBs. A number of recent publications are devoted to the problem, e.g. a possibility of the Tayler instability (Weber et al. 2014; Herreman et al. 2015), the thermal convection (Shen & Zikanov 2016), observation of vortical flow in a LMB model (Kelley & Sadoway 2014), and a simplified model of sloshing in a three layer system (Zikanov 2015). The main motivation of these investigations is to prevent the possibility of direct contact between the molten metal anode and cathode that may occur due to electro-magnetically driven destabilizing interface motion. On the other hand, the controlled mixing enhances mass transport improving the cell performance, preventing the accumulation of intermetallic compounds at the electrode-electrolyte interface.

LMBs are thought to be easily scalable on the cell level due to their simple construction using the natural density stratification of the liquid layers. Large cells of several cubic meters total volume have a potential to operate at very high power value (Bojarevics & Tucs 2017). High current densities coupled to the magnetic field (created by the currents in the cell, the supply bars and the neighbour cells) lead to significant electromagnetic forces. Such forces in stratified liquid layers with large surface areas may cause a long wave interfacial instability as it is well known in the case of Hall-Heroult cells (HHC), as first described by Sele (1977). The manifestation of this instability in LMB is the subject of this paper.

In a typical HHC the electric current, of total magnitude kA, enters the cell from the carbon anodes, passes through the liquid electrolyte and aluminium layer, and exits via the carbon cathode blocks at the bottom of the cell. The liquid layers are relatively shallow, cm in depth vs m in horizontal dimension. The small relative depth of the layers and the small difference of the liquid densities facilitates the instability development.

The ratio of electrical conductivities of the cell materials is another significant parameter. The liquid metal is a better conductor ( S/m) than the carbon ( S/m), while the electrolyte is about two orders of magnitude less conductive ( S/m). The significantly lower conductivity of the electrolyte means that this layer is responsible for the majority of electrical losses in the cell. Joule heating is necessary to heat the cell and to keep the metal liquid, however the total voltage drop must be as low as possible in order to achieve a better electrical efficiency. A small perturbation of the interface between liquid layers may cause a substantial redistribution of the current in the cell.

First attempts to explain the interfacial instabilities were made by Sele (1977), Urata (1985), Sneyd (1985) and Moreau & Ziegler (1986). A more involved understanding of the physical mechanism was provided by Sneyd & Wang (1994), Bojarevics & Romerio (1994), and Davidson & Lindsay (1998). The mechanism is based on the standing gravity wave modification due to the electric current redistribution. The electric current density in the electrolyte increases above the wave crests, resulting in a high density horizontal current in the shallow liquid metal layer. In the presence of a vertical magnetic field the electromagnetic force excites another standing wave mode orthogonal to the initial perturbation. The new wave mode is coupled to the original mode, and the oscillation frequency is shifted. The frequency shift increases with the rise of the magnetic field until at a critical value, when the two wave frequencies coincide, an exponential growth of the amplitude indicates the onset of instability. In general, the above process is described by the following set of equations:

(1.0)

where is a vector which represents the amplitudes of the original gravitational modes , is the matrix of the gravitational frequencies, is the interaction matrix, is the dimensionless parameter characterizing the electromagnetic forces. The mode coupling is included in , where each column represents the Lorentz force (Fourier decomposed) in response to the gravitational wave modes. These coupled equations represent an eigenvalue problem for the square of the new complex frequencies (). The matrix is real anti-symmetric, and in a general case the eigenvalues are shifted increasing the magnetic field (Bojarevics & Romerio 1994). Onset of the instability starts at a critical value of at which the exponentially growing part of the complex eigenvalue appears. See also the paper by Antille & von Kaenel (2002) where the numerical eigenvalue solution is analysed when approaching the instability threshold. A key point noted in the papers, is that the dominant contribution to the perturbed Lorentz force arises from the interaction between a horizontal current in the aluminium layer and the vertical component of the background magnetic field.

Davidson & Lindsay (1998) derived a simple mechanical analogue which captures the basic features of the instability. The liquid aluminium layer is represented by a compound pendulum that consists of a large flat aluminium plate attached to a top surface by a light, rigid strut. The strut is pivoted at its top end so that the plate is free to swing along two horizontal axes and . The fluid system of infinite motion freedom is reduced to only two degrees of freedom. Zikanov (2015) constructed a similar mechanical model for instability description in the LMB taking into account an additional top liquid metal layer. The metal layers of the battery are represented by solid metal slabs rigidly attached to weightless rigid struts pivoted at the top. The free oscillations of the slabs imitate the sloshing motion of the liquid layers. The slabs are separated from each other by a layer of a poorly conducting electrolyte. Two destabilization mechanisms were considered: 1) interaction of a purely vertical magnetic field and horizontal currents, similar to HHC, 2) interaction between the current perturbations and the azimuthal self-magnetic field from the total vertical current. The first mechanism will occur in real batteries if a sufficiently strong vertical magnetic field is present due to the presence of external current supply. The batteries of a square or a circular horizontal cross section will be always unstable if even a small field is present. The second mechanism appears to be more challenging since the azimuthal magnetic field, unlike the vertical magnetic field, cannot be reduced via optimization of the current supply lines unless they cross the liquid layer (Weber et al. 2014). The existence of the second instability type was predicted by Munger & Vincent (2008) for HHC case, yet needs more clarification for the LMB case. The approach developed by Davidson & Lindsay (1998) and Zikanov (2015) is purely mechanical. However, the principal physical mechanism could be valid, due to the fact that sloshing motions generated in the shallow liquid layers are inherently large scale, and so their qualitative behaviour can be approximately described using the coupled pair of long wave modes approach.

More realistic fluid dynamic description can be achieved starting from the full set of Navier-Stokes equations by means of the shallow layer approximation and a systematic derivation of a set of coupled wave equations governing the three fluid layers. The hydrodynamic coupling is realised by pressure continuity at the common interfaces. The continuity of the electric potential and the supplied electric current will introduce the electromagnetic coupling of waves. In the following sections the linear stability of coupled modes will be investigated in the presence of a purely vertical magnetic field, accounting for the continuous electric current in the LMB model. The role of dissipation rate will be analysed using both analytical tools and numerical solutions. Analytical criteria for the cell stability will be established using approximations suitable for a practical cell design, including the solid bottom and top friction effects on the shallow layers.

2 Interfacial dynamics

Hydrodynamics of the three density stratified electrically conductive liquid layer system, schematically represented in figure 1, in the presence of electro-magnetic fields, are described by the following equations

(2.0)
(2.0)

where the indices correspond to the coordinates , the velocity components are given as , the summation over repeated indices is implied, represents density, - effective viscosity, - the gravitational acceleration, - pressure and the vector of Lorentz force is computed in each of the 3 layers, is the current density and the magnetic field. In this paper the horizontal dimensions of the cell are assumed to be much larger compared to the vertical depth, so that the description can be based on a systematically derived shallow layer approximation. The velocity components in each layer can be represented as an expansion in a small aspect ratio parameter , where is a typical depth, for instance the unperturbed metal layer, and is the characteristic horizontal dimension (width of the cell):

(2.0)
(2.0)

where a stretched vertical coordinate is introduced. The expansion starts with the -order due to (2). If all three components of the electromagnetic force density are of the same order of magnitude: , and the horizontal pressure gradient components are of the same order as the corresponding force components: then the vertical component of the gradient . According to these estimates, the leading horizontal components of (2) are

Figure 1: 3 layer liquid metal model under consideration.
(2.0)

The vertical component of the equation (2) gives the leading order terms as:

(2.0)

The hydrostatic pressure in the liquid layers adjacent to the interface , see the figure 1, can be expressed by

(2.0)

where the index stands for the layer number and is the reference pressure at the moving interface . Similarly the same pressure can be referenced to the interface , the corresponding pressure given by

(2.0)

where in this particular case . The respective horizontal gradients of the pressure required in the horizontal momentum equation (2) are:

(2.0)
(2.0)

The next step is to introduce the depth averaging within each layer. The depth averaging for horizontal velocity components is performed in the following way:

(2.0)

where is the layer number (no summation over ) and is the local variable depth, see figure 1. The depth averaging can be applied to the continuity equation (2):

(2.0)

where . The vertical velocity at the is given by the kinematic condition, stating that the interface moves with the local velocity:

(2.0)

Substituting (2.0) into (2.0) leads to

(2.0)

The last equation can be linearised if an additional approximation of a small wave amplitude is introduced: for the layer thickness or equivalently for the interface position, where the additional small parameter is introduced. is a typical wave amplitude and , are the unperturbed values, are the interfacial perturbations. For each particular layer the depth average horizontal velocity divergence can be expressed as:

(2.0)
(2.0)
(2.0)

The volume conservation requires:

(2.0)

In order to estimate the leading order of terms in the depth averaged momentum equations (2) dimensionless variables of order are introduced using the following scaling: - for the coordinates ; for a typical layer thickness; - for the wave velocity; - for the time; for the pressure (). For typical geometries considered in this paper 0.2 m/8 m = , whereas, 0.005 m/0.2 m = .

The horizontal pressure gradient from the expressions (2) and (2), neglecting terms of the and higher order, can be substituted in the depth averaged, nondimensionalized horizontal momentum equation (2). For the layers adjacent to the lower interface the respective momentum equations are

(2.0)
(2.0)

where the depth averaged force is defined similarly to (2). For the upper interface the respective equations are:

(2.0)
(2.0)

The equations (2.0) and (2.0) formally give the connection between the reference pressures and defined on the two interfaces, however being valid in the same fluid layer (electrolyte). The alternative representations are required for the wave equation derivation. After the integration over depth the dissipative terms in (2) are replaced by empirical expressions used for the shallow layer approximation (Moreau & Evans 1984; Rodi 1987) using a linear in velocity friction law with the coefficients . The electromagnetic interaction parameter (the ratio of electromagnetic force to the gravity force perturbation) is introduced as , where is the total electric current, is a typical magnitude of magnetic field. The corresponding magnitude of can be estimated, using typical values for A, T, m (width of cell), kg m (liquid magnesium for the top metal), m s, : . The electromagnetic term is of the same order of magnitude as the leading terms, while the nonlinear wave motion terms are of lower order () and will be neglected later in the linear theory, but retained for a numerical solution.

The wave equations for the coupled interfaces can be derived following the procedure described in Bojarevics (1992):

  1. take time derivative of the non-dimensional linearised equations (2.0),(2.0):

    (2.0)
    (2.0)
  2. substitute (2.0), (2.0) into the horizontal divergence of (2.0),(2.0),

  3. take the difference of the resulting equations.

This procedure eliminates the common unknown pressure on the interface :

(2.0)

The corresponding boundary conditions for the normal velocity at the side walls can be obtained by taking the difference between (2.0) and (2.0) to eliminate the common pressure at the interface :

(2.0)

In a similar manner the wave equation for the upper interface can be obtained:

(2.0)

and the corresponding boundary conditions are

(2.0)

The new constants introduced in the above equations are defined as:

(2.0)
(2.0)
(2.0)

Note, that in (2) and (2) summation over repeated indices is limited to the two horizontal dimensions. As it can be seen from (2) and (2), both interfaces can not be considered independently due to the presence of coupling terms. In the following we will assume that the interfacial friction is negligible. The set of the equations in the electrically nonconductive limit and in the absence of viscous dissipation is in correspondence with the set of the equations obtained by Robino et al. (2001) derived for the dynamics of internal solitary waves in stratified 3 layer ocean. The aluminium electrolysis cell MHD wave model can be recovered if in (2).

3 Electric current flow

For energy storage and supply the LMBs must operate in two regimes: charge and discharge, resulting in the current flowing (upwards or downwards). In this paper only the charging process is considered due to the physical symmetry of both operational regimes. The current flow in the layered structure is illustrated in figure 1. Similarly to Davidson & Lindsay (1998), we assume that the characteristic time-scale for the wave motion is much larger than the diffusion time of the magnetic field to satisfy the low magnetic Reynolds number approximation, leading to , where is the magnetic permeability, the electrical conductivity, is a typical velocity. Using typical values for the liquid metal H m, S m, m s, m the estimated value of . A similar estimate can be obtained for the wave motion magnetically induced electric current ration to the basic current density: .

In the low approximation, when the flow effect is neglected, the electric current can be expressed as

(3.0)

and is described by a set of coupled Laplace equations for the electric potential :

(3.0)

where corresponds to the layer number. The continuity conditions for the electric potential and the normal current component at the interfaces are

(3.0)
(3.0)

The normal derivatives at the deformed interfaces are defined as (assuming the summation over the repeated index only):

(3.0)

With (3.0) the current continuity (3.0) at the interfaces and can be written explicitly in the nondimensional form in order to estimate the leading order terms:

(3.0)
(3.0)

where the four orders of magnitude difference in the electrical conductivities permit to define , and the stretched . These definitions allow us to compare numerically the electrical conductivities in the poorly conducting electrolyte relative to the well conducting liquid metals, and the effect of the small depth () of the layers. The side walls of the domain are considered to be electrically insulating:

(3.0)

In this paper we assume, that the applied current distributions at the top and the bottom are uniform, equal and the interfacial perturbation do not influence the electric current distributions in the collectors:

(3.0)

where . In principle, could be used, however requiring an external circuit solution.

The set of Laplace equations (3.0) can be rewritten in a nondimensional form

(3.0)

The shallow layer approximation requires that the potential is expanded in terms of the parameter :

(3.0)

where the expansion terms are expressed in a similar manner as in Bojarevics & Romerio (1994):

(3.0)
(3.0)
(3.0)

where , , , , , are the coordinate and dependent functions that correspond to the unperturbed interfaces. The functions , , , , , are , and time dependent, corresponding to the perturbed interfaces. Taking into account the previously described boundary conditions and neglecting the higher order terms, the unknown coefficients can be determined as shown in the Appendix A by equalising the similar order of magnitude terms.

Finally, the resulting set of the equations governing the electric current distribution in the system is obtained by introducing the perturbed potentials in both metal layers: , . For application in the LMB the dimensional equations for the electric potential perturbations are linearly correlated to the respective interface perturbations:

(3.0)
(3.0)

where

(3.0)
(3.0)

where is the unperturbed layer thickness. The current distribution in the electrolyte is almost purely vertical due to fact that (similarly to Zikanov (2017)). According to (3.0) and (3.0), the current flow is perturbed by the electrolyte thickness perturbations. Finally, the corresponding dimensional current components can be expressed as

(3.0)
(3.0)
(3.0)

In the following section it will be shown that some of these perturbations may become electromagnetically coupled due to the presence of magnetic field and may lead to an instability resulting in a short circuit state at the extreme case.

4 Linear stability analysis

4.1 Coupled 3 layer problem

The wave equations (2)-(2.0) can be linearised neglecting the -order nonlinear terms and assuming a given magnetic field distribution, which can be expanded in terms of the small parameter : . This expansion follows the same principle as the electric current representation . At this stage we assume that the leading part of the magnetic field is caused by external sources (connectors, supply lines, neighbouring batteries etc.), therefore in the liquid zone. This allows us to neglect the EM force components resulting from the horizontal magnetic field interaction with the vertical unperturbed current: , due to the zero contribution to the horizontal divergence term in (2), (2.0): . The leading horizontal force components contain the horizontal electric current and the vertical magnetic field: , , noting that is -order lower than the horizontal perturbation current. This confirms to the assumptions made in the Introduction and implied in the previous studies on MHD stability of HHC, that the magnetic field is purely-vertical , and it is caused by external sources.

The set of the wave equations for this case has the following form, after assuming that the friction at the electrolyte top and bottom is negligible in comparison to the friction at the solid top and bottom:

(4.0)
(4.0)

with the boundary conditions at the side-walls:

(4.0)
(4.0)

The electric potential distribution is governed by the set of equations (3.0), (3.0). Note, that the coefficients on the left hand side of (3.0) and (3.0) contain only the constant parts of the layer thickness.

The problem can be rewritten in a weak form by means of integrating the equations on the horizontal interface against a set of regular test functions (see the integral representation in Appendix B). The solution can be constructed in a Sobolev space , so that and are satisfying the corresponding equations for all test-functions and that belong to . The following set of functions is introduced:

(4.0)

where is the normalization coefficient. The elements of form orthogonal basis in and the corresponding physical unknowns can be expressed in a similar form as the series:

(4.0)
(4.0)
(4.0)
(4.0)

where , are the spectral wave amplitudes and the perturbed potentials in Fourier space, whereas . Note that the boundary conditions are satisfied in the weak sense only when using the functions (4.0)-(4.0). Taking into account the orthogonality properties of the cosine functions the set of wave equations including the boundary conditions can be rewritten (see the Appendix B) in the spectral coefficient space for and mode interactions:

(4.0)
(4.0)

where the new coefficients are defined as:

(4.0)
(4.0)
(4.0)
(4.0)

The corresponding uncoupled shallow layer gravity wave frequencies are

(4.0)
(4.0)

The selection of the magnetic field modes in (4.0</