Shallow water modeling of rolling pad instability in liquid metal batteries

TL;DR

This paper develops a shallow water model to analyze magnetohydrodynamic interface instabilities in liquid metal batteries, revealing complex dynamics similar to rolling pad instability and interactions between interfacial waves.

Contribution

A new two-dimensional nonlinear shallow water model for liquid metal battery interface instability is introduced and numerically studied.

Findings

Instability resembles rolling pad instability in certain conditions.

Complex dynamics emerge when electrolyte-metal density differences are similar.

Interaction between symmetric and antisymmetric interfacial waves observed.

Abstract

Magnetohydrodynamically induced interface instability in liquid metal batteries is analyzed. The batteries are represented by a simplified system in the form of a rectangular cell, in which strong vertical electric current flows through three horizontal layers: the layer of a heavy metal at the bottom, the layer of a light metal at the top, and the layer of electrolyte in the middle. A new two-dimensional nonlinear model based on the conservative shallow water approximation is derived and utilized in a numerical study. It is found that in the case of small density difference between the electrolyte and one of the metals, the instability closely resembles the rolling pad instability observed earlier in the aluminum reduction cells. When the two electrolyte-metal density differences are comparable, the dynamics of unstable systems is more complex and characterized by interaction between two nearly symmetric or antisymmetric interfacial waves.

Figures

Tables

Table 1. Table 1: Verification test of the instability parameter ( 82 ). Exponential growth rates γ 𝛾 \gamma and periods T 𝑇 T of growing perturbations at viscosity ν A = ν B = ν E = 5 × 10 − 7 superscript 𝜈 𝐴 superscript 𝜈 𝐵 superscript 𝜈 𝐸 5 superscript 10 7 \nu^{A}=\nu^{B}=\nu^{E}=5\times 10^{-7} m 2 /s and the analogous coefficients γ ~ ~ 𝛾 \tilde{\gamma} , T ~ ~ 𝑇 \tilde{T} at ν A = ν B = ν E = 5 × 10 − 11 superscript 𝜈 𝐴 superscript 𝜈 𝐵 superscript 𝜈 𝐸 5 superscript 10 11 \nu^{A}=\nu^{B}=\nu^{E}=5\times 10^{-11} m 2 /s are shown. The base case is the simulation with Π = 5.734 Π 5.734 \Pi=5.734 illustrated in Figs.  2 and 3 . In each of the other 5 simulations, two battery parameters are changed in comparison to the base case so that the value of Π Π \Pi remains the same. Current density J 0 subscript 𝐽 0 J_{0} is in A/m 2 , magnetic field B z subscript 𝐵 𝑧 B_{z} is in T, lengths are in m, and density is in kg/m 3 .

Run	Base	$J_0=5\times {10}^3$	$J_0=5\times {10}^3$	$J_0=2\times {10}^4$	$J_0=2\times {10}^4$	$J_0=2\times {10}^4$
	case	$L_x=1.061$	$B_0=0.002$	$H_0^E=0.01$	$H_0^A=H_0^B=0.2$	$\Delta {\rho }^A=200$
$\gamma$	0.379	0.357	0.378	0.417	0.404	0.384
$T$	0.947	0.946	0.947	0.945	0.947	0.952
$\tilde{\gamma }$	0.439	0.439	0.439	0.440	0.438	0.425
$\tilde{T}$	0.944	0.943	0.944	0.948	0.946	0.952

Equations

$$H^A,H^B,H^E\ll L.$$

$$H^A,H^B,H^E\ll \lambda \sim L.$$

$$\delta \equiv \frac{H}{L}\sim \frac{H}{\lambda }\ll 1.$$

$${\rho }^A<{\rho }^E<{\rho }^B.$$

$${\sigma }^E\ll {\sigma }^A\sim {\sigma }^B.$$

$${\mathbf{J}}_0=J_0{\mathbf{e}}_z,J_0=const.$$

$${\mathbf{B}}_0=B_0{\mathbf{e}}_z,B_0=const.$$

$$z={\zeta }_0^A=\frac{H_0^E}{2},z={\zeta }_0^B=-\frac{H_0^E}{2},$$

$$\eta (x,y,t)=H^E(x,y,t)-H_0^E={\zeta }^A(x,y,t)-{\zeta }^B(x,y,t)-H_0^E$$

$$\left(J_0+j_z^E(x,y,t)\right){\mathbf{e}}_z$$

$$J_0+j_z^E={\sigma }^E\frac{{\Phi }^A-{\Phi }^B}{H^E}\text{ (perturbed state)},J_0={\sigma }^E\frac{{\Phi }_0^A-{\Phi }_0^B}{H_0^E}\text{ (base state)},$$

$$⟨J_0+j_z^E⟩\equiv {\int }_0^{L_y}{\int }_0^{L_x}\left(J_0+j_z^E\right)dxdy=⟨J_0⟩.$$

$$j_z^E={\sigma }^E\left(\frac{C(t){\Phi }_0^A}{H^E}-\frac{{\Phi }_0^A}{H_0^E}\right)=J_0\left(\frac{C(t)H_0^E}{H^E}-1\right),$$

$$C(t)=\frac{⟨J_0⟩}{⟨J_0{H}_0^E/H^E⟩}.$$

$$j_z^E\approx -J_0\frac{\eta }{H_0^E}.$$

$$\nabla \cdot {\mathbf{j}}^{\mathbf{A}}=0$$

$${j_z^E|}_{z={\zeta }^A}={j_z^A|}_{z={\zeta }^A}-{j_x^A|}_{z={\zeta }^A}\frac{\partial {\zeta }^A}{\partial x}-{j_y^A|}_{z={\zeta }^A}\frac{\partial {\zeta }^A}{\partial y}$$

$$\frac{\partial J_x^A}{\partial x}+\frac{\partial J_y^A}{\partial y}=j_z^E,$$

$$J_x^A={\int }_{{\zeta }^A}^{z_t}{j}_x^{A}dz,J_y^A={\int }_{{\zeta }^A}^{z_t}{j}_y^{A}dz$$

$${\nabla }_{\perp }\times \mathbf{J}=0$$

$${\mathbf{J}}_{\perp }^A={\nabla }_{\perp }{\Psi }^A,{\nabla }_{\perp }^2{\Psi }^A=j_z^E,$$

$${\nabla }_{\perp }\times {\mathbf{J}}_{\perp }^A={j_x^A|}_{z={\zeta }^A}\frac{\partial {\zeta }^A}{\partial y}-{j_y^A|}_{z={\zeta }^A}\frac{\partial {\zeta }^A}{\partial x}.$$

$${\mathbf{J}}_{\perp }^B=-{\mathbf{J}}_{\perp }^A.$$

$${\mathbf{F}}_{\perp }^A={\mathbf{J}}_{\perp }^A\times \left(B_0{\mathbf{e}}_z\right),{\mathbf{F}}_{\perp }^A=-{\mathbf{F}}_{\perp }^B.$$

$$J_z^E=H^E{j}_z^E,J_z^A=\frac{H^A}{2}{j}_z^E,J_z^B=\frac{H^B}{2}{j}_z^E,$$

$${\mathbf{f}}^A$$

$${\mathbf{f}}^E$$

$${\mathbf{f}}^B$$

$$\frac{\partial \zeta }{\partial t}+u\frac{\partial \zeta }{\partial x}+v\frac{\partial \zeta }{\partial y}=w,$$

$$w\sim \delta u\sim \delta v.$$

$$\frac{\partial p}{\partial z}=-\rho g+f_z.$$

Full text

1 Introduction

The liquid metal battery is a promising conceptual device for stationary energy storage Kim:2013 . Recently developed small-scale (about few cm) laboratory prototypes have demonstrated technical feasibility of the concept and its advantages of higher efficiency and longer operational life in comparison to traditional solid-electrode batteries (see e.g. Bradwell:2012 ; Kim:2013calcium ; Wang:2014 ; Ouchi:2016 ; Xu:2016 ). Scaling up the concept to large commercially attractive devices has not, however, yet been achieved. As we discuss below, the hydrodynamic effects, in particular that of electromagnetically modified interfacial waves, are expected to be non-negligible factors of this transition.

In this paper we consider one version of the battery, which, on the level of simplification sufficient for out purposes, can be viewed as the system illustrated in Fig. 1a. It is a cuboid or cylindrical cell filled with three liquid layers: the layer of a heavy metal (e.g. Bi, Sb, Zn, or PbSb) at the bottom, the layer of a light metal (e.g. Na, Li, Ca, or Mg) at the top, and the thin layer of molten-salt electrolyte sandwiched in the middle. The electrolyte is selected so that it is immiscible with the metals and conductive to positive ions of the light metal. It has the density intermediate between those of the light and heavy metals, so the entire system is stably stratified by density. The sidewalls are electrically insulated, while the horizontal top and bottom walls are conducting and serve as current collectors. The system is maintained at a temperature above the melting points of all the three materials (200 to 700*∘*C depending on the materials used).

The energy stored in the battery is the difference between the Gibbs free energies of the light metal in its pure state and the state of alloy with the heavy metal. The battery is charged when the light metal is electrochemically reduced from the alloy and discharged when the alloy is formed. The reactions occur in liquid state in the bottom layer in the presence of strong ($\sim$ 1 A/cm2) electric currents flowing in the vertical direction.

The illustration in Fig. 1a lets a fluid dynamicist to recognize immediately that the hydrodynamic instabilities can be a major factor of the system’s operation, especially when the transition from small laboratory prototypes to large commercial devices is made. The effect of the instabilities may be positive, e.g. when the resulting flow in the bottom layer enhances mixing of reactants. It can also be negative, when the deformation of the interfaces between the layers becomes so strong that it causes rupture of the electrolyte layer, which means short circuit between the metals and disruption of the battery’s operation.

Several mechanisms of instability have been identified and analyzed so far including the Tayler instability Stefani:2011 ; Seilmayer:2012 ; Weber:2014 ; Herreman:2015 , the thermal convection caused by the volumetric Joule heating of the electrolyte Shen:2016 ; Xiang:2017 or by bottom heating Kelley:2014 , and the rolling pad instability Zikanov:2015 ; Weber:2016 ; Bojarevics:2017 ; Horstmann:2017 . We should also mention the electrovortex effect, which may become a significant factor under certain conditions Weber:2015 . The analysis is far from complete, but it can already be said that all these mechanisms are likely to be active in large-scale batteries. At the same time, the preliminary estimates made in Herreman:2015 for the Tayler instability and in Shen:2016 for the Joule-heating convection show significant damping effect of the stable density stratification among the layers. It is possible, although not yet proven, that in a wide range of operational parameters, the instabilities are present but do not lead to a disruption. Rather, the instabilities lead to saturated states with weak flows and small amplitudes of interface deformation.

In this paper, we study the rolling pad instability. It is conceptually similar to the instability observed and well studied in the Hall-Héroult aluminum reduction cells (see e.g. Davidson:2016 ). The reduction cell is a shallow and large rectangular bath filled with molten aluminum at the bottom and molten salt electrolyte with alumina dissolved in it at the top. Nearly vertical currents of density $\sim$ 1 A/cm2 pass through the two layers and cause the desired electrochemical reduction of aluminum from its oxide.

The reduction cell is different from the liquid metal battery in many respects: two layers instead of three, turbulent mixing of melts driven by gas bubbles and Lorentz forces, the top wall consisting of many separate carbon anodes, etc. The one important similarity is that in both the systems, the liquid domain consists of horizontal layers with vastly different electric conductivities. The conductivity of molten-salt electrolytes is about four orders of magnitude lower than that of liquid metals. This means that even a small variation of the local thickness of the electrolyte causes large variation of local electric resistance and, thus, significantly changes the distribution of the electric currents within the system. Lorentz forces are created by the current perturbations interacting with the magnetic field inevitably present within the cell. In unstable systems, the forces modify the flows of the melts in such a way that the interface deformation is enhanced.

In the aluminum production, the instability is viewed as a major concern. It develops in a reduction cell when the thickness of the electrolyte is reduced below a certain cell-specific threshold and leads to sloshing waves that grow and, eventually, cause short circuit between the aluminum and the anode. The thickness has to be kept large enough, which results in energy losses to the excessive Joule heating of the electrolyte. The situation improved in the 1980s and 1990s when it was understood that the instability was caused by the interaction between the horizontal components of the current perturbations and the vertical component of the magnetic field generated by the currents flowing in the neighboring cells and electric supply lines.

A detailed discussion of the instability mechanism can be found, e.g. in Bojarevics:1994 ; Sneyd:1994 ; Davidson:1998 ; Davidson:2016 . Briefly, as illustrated in Fig. 1b, in which we should ignore the metal layer A for to consider a reduction cell, the vertical current perturbations $\bm{j}_{z}$ form as a result of the variation of the local thickness of the electrolyte caused by a deformation of interface. The current perturbations close within the highly conducting metal layer resulting in a horizontal component $\bm{j}_{\bot}$. The interaction between this component and the base vertical magnetic field $\bm{B}_{0}$ create the horizontal Lorentz force $\bm{f}$ that drives the metal in the horizontal direction perpendicular to the plane of the interface deformation. The result is a large-scale interfacial wave structure rotating around the cell - a so-called rolling pad.

The mechanism by which the rotating wave grows in amplitude, and, thus, the instability develops is easiest to understand using the linearized shallow water model. In this case, the problem is reduced to a single wave equation (see e.g. Davidson:1998 ). Application of the Fourier expansion shows that the Lorentz forces modify the natural gravitational standing waves on the interface and introduce coupling between them. In the unstable situation, the eigenvalues of one or several couples of such waves merge to form complex-conjugate pairs, each corresponding to a rotating and exponentially growing wave. The shape of the horizontal cross-section of the cell determines the natural gravitational modes and, thus, strongly affects the threshold of the instability. In particular, cells of square or circular cross-section are degenerate in the sense that they have multiple gravitational modes with equal frequencies. In the inviscid limit, such cells are unstable at arbitrarily weak electromagnetic effect. For rectangular cells, the strength of the effect must exceeds a certain not very large threshold. In the real cells, the threshold is increased by the viscosity and interfacial tension, and the situation is further complicated by the background melt flows, gas bubbles, block structure of the anode, and other factors, but the principal instability mechanism remains the same.

Fig. 1b illustrates the evident fact that a similar instability can develop in a liquid metal battery. The presence of the second metal layer, in which perturbation currents $\bm{j}_{\bot}$ of opposite orientation develop (we discuss this in detail in section 2.2), makes the picture more complicated. The first attempt to analyze the instability was made in Zikanov:2015 . A mechanical model based on the approach developed earlier for the aluminum reduction cells in Davidson:1998 was used. In the model, the sloshing motions of the metal layers are imitated by pendulum-like oscillations of two solid metal slabs, which are independently suspended and separated from each other by an electrolyte. The oscillations are modified and coupled to each other by the same physical mechanism as in the battery, i.e. via the Lorentz forces associated with the current perturbations caused by the local changes of the electrolyte thickness. The model is a drastic simplification of a real battery, but provides what can be viewed as a low-mode, nondissipative, linearized analogy of some aspects of the instability.

The main result of Zikanov:2015 is the clear indication that the rolling pad instability similar to the instability in the reduction cells should be expected in the batteries. The model also suggests existence of an instability caused by the interaction between the current perturbations and the azimuthal magnetic field induced by the base current.

The existence of the rolling pad instability in the liquid metal batteries was recently confirmed in the three-dimensional Weber:2016 and two-dimensional shallow water Bojarevics:2017 numerical simulations. Batteries of cylindrical Weber:2016 or rectangular Bojarevics:2017 shapes were considered. The typical fluid densities of a Mg-Sb battery were chosen, for which the density jump and, thus, the stabilizing effect of the buoyancy force is much larger at the lower than at the upper interface. The situation becomes particularly close to that of a reduction cell in the sense that the instability leads to noticeable deformation of only one interface. As another indication of the closeness, it has been found in Weber:2016 that the onset of the instability is best determined by the generalization of the stability criterion developed for the reduction cells in Sele:1977 and Urata:1985 .

This paper reports the results of a new computational study. We present a detailed derivation of the nonlinear shallow water model of a liquid metal battery. The model is then applied in a numerical study of the rolling pad instability. The goals are to (i) further analyze the instability appearing at a single interface in the situation of strongly different density jumps and (ii) extend the analysis to the more complex case when the density jumps are comparable to each other, and so both the interfaces are significantly deformed.

As we discuss in detail in the model derivation in section 2, its validity is limited to shallow cells, in which the horizontal dimensions are much larger than the height. This limitation does not appear especially strict when we consider the main technological requirements to a large commercial battery. It must contain large quantities of both metals, have good vertical mixing in the bottom layer. The electrolyte layer must be thin and have large area in the horizontal plane, so as to maximize the reaction rate and minimize the Joule heat losses. This suggests the optimal shape of the battery in the form of a horizontally large (several m) and shallow (perhaps tens of cm) cell. The additional benefit of such a shape can be avoidance of disruptively strong Tayler and thermal convection instabilities.

2 Physical model and numerical method

A battery in the form of a shallow cell shown schematically in Fig. 1a is considered. The cell is filled with three liquid layers: metal B at the bottom, metal A at the top, and electrolyte E in the middle. The sidewalls are vertical and electrically insulated. The top and bottom walls are horizontal and serve as the electrodes, between which strong electric current of density $\bm{J}$ flows through the cell.

In the following discussion, variables with the superscripts A, B, or E are associated with particular layers. Variables without superscripts either mark universal variables common for all layers (such as the horizontal length scale $L$) or, for the sake of brevity, stand for the respective variables in each layer (e.g. $H$ standing for the layer thicknesses $H^{A}$, $H^{B}$, and $H^{E}$).

2.1 Assumptions made in the model

The model describes the behavior related to the interaction between the magnetic fields and the perturbations of electric current arising in response to the interface deformation. The physical effects leading to the other instabilities mentioned in section 1 are ignored. The justification for this approach is based on the following two arguments. The small-scale, possibly turbulent motions of melts caused by thermal convection and interfacial forces increase the dissipative effects, but are unlikely to substantially change the physical mechanisms of the large-scale rolling pad instability. The validity of this argument is confirmed by the results obtained for the aluminum reduction cells (see, e.g. Zikanov:2000 ; Moreau:1984 ) and further supported in the following discussion. Furthermore, as we discuss below, the Tayler instability and the electrovortex flow are either absent or weak in the simplified shallow geometry we consider.

The model is based on the following simplifying assumptions.

The cell is shallow in the sense that the thicknesses of all three layers are much smaller than the typical horizontal size of the cell:

[TABLE]

This is the strongest of our assumptions, and the only one that seriously limits the applicability of the model. We further assume that the typical horizontal wavelength $\lambda$ of the instability is large:

[TABLE]

The error of approximation of the model can be expressed in terms of the shallowness parameter

[TABLE]

As we will see in the following discussion, it is partially of the second order (e.g. for the inviscid part of the momentum transport) and partially of the first order (e.g. for the effects of viscosity or some electromagnetic effects). 2. 2.

The metals are immiscible with the electrolyte, so sharp deformable interfaces exist between the layers. 3. 3.

Each layer is assumed to contain a liquid of constant physical properties. This is a simplification, especially in the case of the bottom layer, where a mixture of metal B with the alloy of A and B in liquid or, possibly, solid intermetallic state is to be found unless the battery is fully charged. 4. 4.

The system is stably stratified by density:

[TABLE] 5. 5.

The electric conductivities of the liquids satisfy

[TABLE]

In real batteries, $\sigma^{E}$ is about four orders of magnitude smaller than $\sigma^{A}$ and $\sigma^{B}$. 6. 6.

The base electric current flowing through a cell with unperturbed horizontal interfaces is purely vertical and uniform, so its density is

[TABLE]

This means that the electrovortex effect is excluded from consideration. 7. 7.

A simplified model of the magnetic field is used. We limit the analysis to the effect of the vertical magnetic field that is generated in a battery by the currents in external circuits. The field is approximated as

[TABLE]

The other components, in particular the magnetic field generated by the base current $\bm{J}_{0}$ within the cell and the perturbation field $\bm{b}$ induced by the current perturbations are not considered. This does not allow us to analyze full three-dimensional dynamics of the battery and excludes from consideration the Tayler instability, the effect of $\bm{b}$ on the rolling pad instability, and the possible second type instability suggested by the mechanical model Zikanov:2015 . The rationale of the approach is as follows. Firstly, we leave the simulations of the full three-dimensional dynamics of the system to future studies. The focus of our work is on the rolling pad instability during its low-amplitude stages when it can, as an approximation, be separated from the other magnetohydrodynamic effects. A similar approach was used for analyzing the classical rolling pad instability in the reduction cells Bojarevics:1994 ; Sneyd:1994 ; Davidson:1998 ; Sele:1977 ; Urata:1985 ; Sun:2004 ; Zikanov:2000 and in the recent studies of liquid metal batteries Weber:2016 ; Bojarevics:2017 ; Horstmann:2017 . Secondly, the Tayler instability is not anticipated in our system, since it is associated with tall cells (as shown, e.g. in Rudiger:2007 , the typical axial wavelength of the instability is larger than the horizontal dimension of the cell). Thirdly, the general effect of the magnetic field perturbations $\bm{b}$ on the rolling pad instability, while present, is likely to be insignificant. This is indicated by the results obtained in Sun:2004 for the reduction cells. Finally, a proper analysis of the instability of the second type predicted in Zikanov:2015 would require accurate evaluation of the magnetic fields induced by both the base and perturbation currents, which can only be done in the framework of the three-dimensional model that includes the interior of the cell and the adjacent conductors. This analysis is left to future studies, as well. 8. 8.

As explained in detail in section 2.2, only the electric current perturbations caused by the deformation of the interfaces are considered. This is justified since, as demonstrated in the three-dimensional simulations Weber:2016 and confirmed by our results presented below, in an unstable system these perturbations are typically much stronger than the electric currents induced by the melt velocities. One consequence of this assumption is that the current perturbations can be represented as gradients of electric potential functions. Full Ohm’s law may need to be applied in a three-dimensional nonlinear analysis of melt flows in a battery with full three-dimensional magnetic field. 9. 9.

The melts are assumed to be at a constant temperature. The thermal convection flows are neglected despite the fact that, as discussed in Shen:2016 , they are practically unavoidable and likely to be turbulent. This is justified by the estimates (see e.g. Shen:2016 ) showing that the stable density stratification between the layers is sufficiently strong to prevent significant interface deformation by the convection in all, but very large batteries. Our assumption is that the convection results in small-scale, possibly turbulent flows but does not critically interfere with the large-scale rolling pad instability. 10. 10.

We now describe the base (unperturbed) state, which is designated in the following discussion by subscript 0. The simplifying assumptions made above imply that the base state has purely vertical current (6) and flat interfaces (see Fig. 1a):

[TABLE]

where $H_{0}^{E}$ is the unperturbed thickness of the electrolyte layer. The Lorentz force $\bm{J}_{0}\times\bm{B}_{0}=0$, and the melt velocities are all zero. This form of the base state is a simplification. In the real battery, $\bm{J}_{0}$ is not uniform and $\bm{B}_{0}$ is a three-dimensional field created by the currents within and without the cell. This and other features of the battery system, for example the thermal convection, inevitably lead to a base state with deformed interfaces and possibly turbulent melt flows. The analysis of the effect of a nontrivial base state on the instability is left to future studies. We note that the model derived below can be applied to systems with background flows and interface deformations after minor modifications (see Sun:2004 ; Zikanov:2000 for examples of such analysis in the case of aluminum reduction cells).

As a final comment, we note that the model is nonlinear and can be applied to simulation of finite-amplitude as well as small-amplitude perturbations.

2.2 Electric current perturbations and Lorentz forces

In the model, we consider the Lorentz forces that arise as a result of the interface deformations, more specifically, as a result of the local change in the thickness of the electrolyte layer

[TABLE]

(see Fig. 1a). The following derivation is a generalization of the derivation used for the aluminum reduction cells (see e.g. Sun:2004 ). The starting point is the approximation based on the relation (5). Due to the much higher electric conductivities of the metals, the distribution of currents in the electrolyte can be accurately approximated as corresponding to equipotential interfaces $z=\zeta^{A}(x,y,t)$ and $z=\zeta^{B}(x,y,t)$. Combined with the shallowness of the electrolyte layer, this allows us to approximate the perturbed current in this layer as vertical (following the shortest path through a resistive medium) and, due to their zero divergence, $z$-independent:

[TABLE]

and to express it via the Ohm’s law as

[TABLE]

where $\Phi^{A}$ and $\Phi^{B}$ are the values of the electric potential at the top and bottom interfaces.

We now set $\Phi^{B}=\Phi^{B}_{0}=0$ and notice that the total electric resistance of the electrolyte changes by the quantity $\sim\left(\eta/H^{E}_{0}\right)^{2}$. Either the net current through the cell or the voltage drop has to change when $\eta(x,y,t)$ is not infinitesimal. Following the approach of Zikanov:2000 , we assume that the voltage drop adjusts as $\Phi^{A}=C(t)\Phi^{E}_{0}$ while the total current remains the same (this corresponds, e.g., to the situation when the considered cell is one of many connected in parallel):

[TABLE]

This allows us to determine the current perturbations as

[TABLE]

where the non-dimensional adjustment coefficient is

[TABLE]

At infinitesimal perturbations $\eta\ll H_{0}^{E}$, the coefficient becomes $C\approx 1$, and the expression (13) can be approximated using the first-order Taylor expansion in $\eta/H_{0}^{E}$ as

[TABLE]

As the next step, we assume that the perturbation currents completely close within the metal layers A and B. For the currents in the layer A, we integrate

[TABLE]

from $z=\zeta^{A}$ to the top boundary $z=z_{t}$. Using the Leibnitz integration rule, the condition that the vertical perturbation current at $z=z_{t}$ is zero, and the condition of continuity of normal current at the interface

[TABLE]

we obtain

[TABLE]

where

[TABLE]

are the vertically integrated horizontal components of the perturbation current in the layer A.

The solution for $\bm{J}^{A}$ can be found in a straightforward manner if we assume that

[TABLE]

and express (18) as

[TABLE]

where $\bot$ marks the horizontal $(x,y)$ plane.

To verify (20), we consider that the three-dimensional current perturbations are curl-free (see assumption 8 in section 2.1). Integration of $\nabla\times\bm{j}^{A}=0$ from $z=\zeta^{A}$ to $z=z_{t}$ gives

[TABLE]

The right-hand side of (22) is of the second order in terms of the perturbation amplitude, so it can be neglected when linear stability to infinitesimal perturbations is analyzed. The situation is more complex for perturbations of finite amplitude. The individual terms in the right-hand and left-hand sides of (22) are of the same order in terms of the shallowness parameter (3). This follows from the estimates $j^{A}\sim HJ^{A}$, $\zeta^{A}\sim H$, and $\partial/\partial x\sim\partial/\partial y\sim\lambda^{-1}$. A closer consideration, however, reveals that the entire right-hand side of (22) is likely to be negligibly small. It is a cross-product of two two-dimensional vector fields: $\left.\bm{j}^{A}\right|_{z=\zeta^{A}}$ and $\nabla_{\bot}\zeta^{A}$. The currents $\left.\bm{j}^{A}\right|_{z=\zeta^{A}}$ are flowing in the direction parallel to the local steepest gradient of the perturbation current $j_{z}^{E}$ injected into the metal. This direction is, in turn, aligned with the direction of $\nabla_{\bot}\eta$. This follows from (13) and is immediately seen in the linearized version (15). If the density differences between the metals and the electrolyte satisfy $\rho^{E}-\rho^{A}\ll\rho^{B}-\rho^{E}$, so only the upper interface is significantly deformed (see section 3.1), $\nabla_{\bot}\eta\approx\nabla_{\bot}\zeta^{A}$. The situation is less certain when the density jumps and, thus, the deformations of the two interfaces are comparable (see section 3.2). For this case we can argue that the interfacial waves are nearly anti-symmetrically coupled (see Horstmann:2017 and our results in section 3.2), so $\nabla_{\bot}\eta$ and $\nabla_{\bot}\zeta^{A}$ are nearly parallel.

In summary, we see that the vectors $\left.\bm{j}^{A}\right|_{z=\zeta^{A}}$ and $\nabla_{\bot}\zeta^{A}$ are perfectly or nearly parallel in the flow. We can neglect the right-hand side of (22) and consider the field of integrated currents $\bm{J}^{A}$ as curl-free. While highly plausible and certainly correct in many cases, this approximation may become inaccurate in some situations, in which case the model becomes only valid for infinitesimal perturbations.

For the vertically integrated perturbation currents in the layer B, we can derive similar relations with $-j_{z}^{E}$ in the right-hand side of the Poisson equation or simply observe that the conservation of electric charge requires

[TABLE]

In our case of purely vertical magnetic field, the Lorentz forces are zero in the layer E and purely horizontal in the layers A and B. The vertically integrated forces are

[TABLE]

In the general case of a three-dimensional magnetic field $\bm{B}$ we would need to use the general three-dimensional expression $\bm{F}=\bm{J}\times\bm{B}$. In addition to the horizontal currents (21), (23), the field of the perturbation currents $\bm{J}$ would include the vertically integrated vertical components

[TABLE]

where the last two expressions are approximations with the error $\sim\delta^{2}$.

2.3 Shallow water approximation

The derivation of the two-dimensional model follows the principal steps used for shallow flows in hydrology. There are differences, in particular the presence of Lorentz forces and immiscible layers of distinct physical properties in our case. Similarly to the approach applied in Sun:2004 to the aluminum reduction cells, the shallow water equations in the conservative de St. Venant form are derived and used. The main advantage of this approximation over the simpler non-conservative shallow water model is that its error in representing the inviscid flow’s dynamics is of the second rather than first order in terms of the shallowness parameter (3). As discussed in Sun:2004 , the two models produce nearly identical solutions in the case of weak interface deformations, but diverge noticeably for the nonlinear states, where the deformations are strong.

The derivation starts with the full three-dimensional momentum equations for flows in each layer, with the liquids treated as incompressible, inviscid and having constant and uniform physical properties. The effect of viscosity will be introduced at a later stage. For the sake of generality and possible future use, the equations are derived for the case of a general three-dimensional Lorentz force that would exist in a cell with a general three-dimensional magnetic field. With the purely vertical currents in the electrolyte and three-dimensional currents in the metals, the general force field has the form

[TABLE]

The reduction to the actually considered in this paper case of a purely vertical magnetic field, in which the force is two-dimensional, will be done later.

The boundary conditions are those of zero velocity at all walls. At the interfaces $z=\zeta^{A}(x,y,t)$ and $z=\zeta^{B}(x,y,t)$, continuity of pressure and the kinematic condition

[TABLE]

where $u$, $v$, and $w$ are the velocity components, are required.

The incompressibility implies the estimate of the vertical velocity

[TABLE]

Applying it to the $z$-momentum equation and dropping all the terms $\sim\delta^{2}$ we obtain, in each layer:

[TABLE]

Next, we introduce the pressure distribution at the mid-plane of the electrolyte layer $z=0$

[TABLE]

and integrate (31) in $z$. The resulting pressure distributions satisfying the requirement of continuity at the interfaces are:

[TABLE]

They are substituted into the horizontal momentum equations written in conservation form for each layer. Together with the incompressibility conditions, the equations are:

[TABLE]

The approaches of the simple non-conservative shallow water model and the conservative de St. Venant model diverge at this point. In the former, we would just assume that $w=0$ and all the other variables are $z$-independent. This would be equivalent to neglecting all but the zero-order terms in the $\delta$-expansions of the flow fields. Instead, we follow the routine of the second order in the shallowness parameter (3) based on the integration of the equations vertically across each layer.

The integration requires the use of the Leibniz rule for the $x$-, $y$-, and $t$-derivatives, the zero-velocity conditions at $z=z_{b}$ and $z=z_{t}$, and the interface conditions (29). We use the velocity fluxes

[TABLE]

and the integrated horizontal components of the Lorentz forces

[TABLE]

We also need the vertically integrated terms of the horizontal gradients of the pressure fields (33)-(35), of which only the Lorentz force terms are non-trivial. For them, we apply the Leibniz rule and the approximation to the second order in $\delta$ to obtain

[TABLE]

where

[TABLE]

In the integration of the nonlinear terms of the momentum equations, we utilize the second order approximations, such as

[TABLE]

The final momentum equations are

[TABLE]

where $\bm{\tau}^{A}$, $\bm{\tau}^{B}$, and $\bm{\tau}^{E}$ are the viscous friction terms, which we introduce here for the first time and discuss in detail below, and the material derivative in the left-hand side of each equation has the $x$- and $y$-components

[TABLE]

The incompressibility condition is also integrated using the Leibniz rule and the interface and wall boundary conditions. This leads to

[TABLE]

In our model, (50) and (51) are applied to determine the evolution of interfaces. The equation (52) is replaced by the combination of (50)-(52) expressing the global conservation of mass

[TABLE]

The approximation of the viscous friction terms follows the approach commonly applied in hydrology. We retain the horizontal components of the Laplacian and approximate the friction at the horizontal walls and interfaces by finite differences. The friction terms in (45)–(47) are:

[TABLE]

This friction model is decisively simple. In particular, the coefficients $\mu$ are the molecular dynamic viscosity coefficients of the respective liquids. No attempt is made to account for the possible small-scale turbulence by using eddy viscosities or applying the models derived for turbulent open channel flows. The main justification of this approach is that, as expected for large-wavelength instabilities and confirmed by our results in section 3, the effect of the viscous friction terms on the evolution of the large-scale modes of the rolling pad instability is limited to minor quantitative changes.

In the computational analysis presented in section 3, the equations are simplified for the case of a purely vertical magnetic field by setting all the vertical force components $F_{z}$ and the forces in the electrolyte $\bm{F}_{\bot}^{E}$ to zero and using the expressions (24) for the remaining horizontal components. The perturbation electric currents are determined as described in section 2.2.

2.4 Equations in non-dimensional form

As a physical system, even the simplified battery considered in this paper is described by a large number of parameters: horizontal dimensions, thicknesses and liquid properties of the three layers, base electric current and magnetic field. Converting the problem into non-dimensional form reduces this number, but still leaves many non-dimensional groups. Furthermore, since the battery is a real device, each physical parameter should be allowed to vary only within some practically meaningful limits. In such a situation, it is convenient to perform the analysis using dimensional variables. At the same time, as we will see in section 3, the non-dimensional form of the problem helps to understand the physics of the instability. Our discussion takes the hybrid approach. The simulation cases are identified using dimensional parameters, while the discussion is conducted primarily in terms of non-dimensional quantities. The non-dimensionalization is presented in this section.

We take the larger horizontal dimension of the cell, say, $L=L_{x}$ in a rectangular cell, as the length scale, the unperturbed current $J_{0}$ as the scale for the density of electric current, the constant vertical component $B_{0}$ as the scale of the magnetic field, and the physical properties of one liquid (we take the electrolyte) $\rho^{E}$, $\nu^{E}=\mu_{E}/\rho_{E}$ as the scales of density and kinematic viscosity. As the velocity scale, we use the typical velocity of a large-scale gravitational wave on the upper interface:

[TABLE]

This choice, also practiced for the aluminum reduction cells (see, e.g. Davidson:1998 ), is consistent with the physical nature of the instability. The typical scales of time and pressure are $L/U_{0}$ and $\rho^{E}U_{0}^{2}$, respectively. The typical scales for the Lorentz force and the electric potential $\Psi^{A}$ are $J_{0}B_{0}L$ and $J_{0}L^{2}$, respectively.

Substituting the scaled variables into the equations of sections 2.3 and 2.2 and performing the non-dimensionalization procedure, we obtain the final non-dimensional system. It includes the momentum equations written for each of the three layers, the equations of the interface deformations, and the incompressibility condition:

[TABLE]

where all the variables are now non-dimensional and the material derivatives in the left-hand sides of (58)-(60) are expressed as in (48)-(49). The non-dimensional friction terms are

[TABLE]

The Lorentz forces are expressed in terms of the non-dimensional variables as

[TABLE]

where $A$ is the non-dimensional area of the horizontal cross-section of the cell. For the system analyzed in this paper, the non-dimensional magnetic field is $\bm{B}_{0}=\bm{e}_{z}$, so only the horizontal components of the Lorentz force are non-zero, and the model for the current perturbations reduces to (67), (70), (71).

The boundary conditions include the no-slip and perfect electric insulation conditions at the walls:

[TABLE]

The non-dimensional parameters of the problem are the Froude number

[TABLE]

the electromagnetic parameter

[TABLE]

the Reynolds number

[TABLE]

and the ratios of densities and kinematic viscosities

[TABLE]

We should also add the parameters defining the cell’s geometry: the non-dimensional unperturbed layer thicknesses $\tilde{H}_{0}^{A}$, $\tilde{H}_{0}^{B}$, $\tilde{H}_{0}^{E}$, and the parameter describing the shape in the horizontal plane, for example, the aspect ratio

[TABLE]

for a rectangular cell.

2.5 Numerical method

The problem is solved numerically using a finite-difference scheme similar to the scheme applied in Zikanov:2000 and Sun:2004 to the shallow water equations for aluminum reduction cells. The time discretization is based on the explicit time-splitting (projection) scheme of the third order. Each time step includes

Solution of the Poisson equation (70) and computation of the electric current perturbations and Lorentz forces, 2. 2.

Advancement of the time-evolution equations (58)-(62) according to the stiffly-stable scheme of Karniadakis:1991 :

[TABLE]

where $f$ stands for a velocity flux or interface deformation and $q$ for the respective right-hand side of the evolution equation. For the momentum equations, the pressure gradient term is not included into the right-hand side, and $f_{*}$ is an intermediate velocity flux. For the interface equations, $f_{*}$ is the value at the next time step. 3. 3.

Solution of the pressure equation

[TABLE]

with zero normal gradients of $p_{0}$ at the sidewalls as the boundary condition. 4. 4.

Correction of the velocity fluxes, in each layer, according to

[TABLE]

and enforcement of the no-slip boundary conditions.

The spatial discretization is of the second order of approximation and uses central differences implemented on a staggered grid. The Poisson equation (70) and the non-separable elliptic equation (80) are solved using the multi-grid algorithms of the library Mudpack Mudpack .

A grid sensitivity study was carried out. It has been found that the large-scale character of the instability allows us to obtain accurate results on moderately fine grids. For example, for a battery of rectangular cross-section with the aspect ratio $\Gamma=2.0$, the grid of $128\times 64$ points is certainly sufficient, while the grid of $64\times 32$ points provides the results, which are generally correct with only minor quantitative inaccuracy.

The analysis consists of multiple simulations, each reproducing the flow’s evolution for a certain set of parameters. Each simulation begins with zero melt velocities and the flat interfaces (8) perturbed randomly with the amplitude of $10^{-5}H^{E}_{0}$. The simulation continues until the maximum deformation of one of the interface exceeds $0.5H^{E}_{0}$ or the growing perturbations reach finite-amplitude saturation, or, in the stable case, for not less than 3000 non-dimensional time units $L_{x}/U_{0}$, with $U_{0}$ given by (57).

3 Results

The analysis consists of two parts different from each other by the relation between the density jumps across the interfaces $\Delta\rho^{A}=\rho^{E}-\rho^{A}$ and $\Delta\rho^{B}=\rho^{B}-\rho^{E}$. The simpler case of $\Delta\rho^{A}\ll\Delta\rho^{B}$ or $\Delta\rho^{B}\ll\Delta\rho^{A}$, in which the significant deformation only occurs at the interface with the smaller jump is considered in section 3.1. The main goal is to continue the analysis of Weber:2016 ; Bojarevics:2017 toward a more accurate evaluation of the extent, to which the predictions developed for the aluminum reduction cells apply to the batteries. The more complex case when $\Delta\rho^{A}\sim\Delta\rho^{B}$, so both the interfaces are significantly deformed, is discussed in section 3.2.

3.1 Single interface deformation

As a typical example of the instability observed at $\Delta\rho^{A}\ll\Delta\rho^{B}$, Figs. 2 and 3 show the flow in a rectangular cell with the aspect ratio $\Gamma=2.0$, $L_{x}=0.75$ m, $\rho^{A}=1000$ kg/m3, $\rho^{B}=8000$ kg/m3, $\rho^{E}=1100$ kg/m3, $J_{0}=10^{4}$ A/m2, $B_{0}$=0.001 T, $H_{0}^{A}=H_{0}^{B}=0.1$ m, $H_{0}^{E}=0.005$ m, and $\nu^{A}=\nu^{B}=\nu^{E}=5\times 10^{-7}$ m2/s.

After an initial adjustment, the time signals of the local and integral characteristics show combined periodicity and exponential growth. Analysis of the interface deformations such as those in Figs. 3a,b shows that these features are caused by, respectively, rotation and growth of an interfacial wave. At the positive (directed upward) base current $J_{0}$, the rotation is counterclockwise if viewed from above. A negative base current changes the sense of rotation, but does not affect the flow in any other way.

The solution in Figs. 2 and 3 is qualitatively similar to the solutions obtained in Weber:2016 ; Bojarevics:2017 in terms of the spatial shape, dynamics, and typical length and time scales of the growing perturbations. The amplitude of the deformation is much smaller at the lower than at the upper interface, approximately in proportion to $\Delta\rho^{A}/\Delta\rho^{B}$.

The instability can be characterized by two coefficients computed during the phase of exponential growth: the time period of oscillations $T$ and the growth rate $\gamma$. Both can be estimated using the time signals shown in Fig. 2. For example, the point signal of the interface deformation in Fig. 2a can be used to determine $T$, while an exponential fit of the appropriate segments of the curves in Fig. 2b,c provides $\gamma$. For the specific flow in Fig. 2, the values are $T=0.947$ and $\gamma=0.379$.

The unstable wave illustrated in Figs. 2 and 3 grows until the maximum interface deformation exceeds $0.5H_{0}^{E}$, at which point the simulations are stopped. Such a growth is typically observed in our simulations in strongly unstable systems. In systems with positive but small growth rate $\gamma$, the perturbations may saturate nonlinearly at some finite amplitude of the interface deformation, typically six to three times smaller than $H_{0}^{E}$.

The evident similarity with the instability in the aluminum reduction cells indicates that the observed instability may be controlled by the same non-dimensional parameters: the horizontal aspect ratio $\Gamma$, the instability parameter, which we discuss shortly, and the Reynolds number, which appears since viscosity is included in our model.

An expression for the parameter controlling the instability is suggested by the studies of the reduction cells. The simplified models, in which viscosity and other complicated effects are neglected, and modified wave equations for low-amplitude large-scale perturbations are derived, all produce principally the same results (see Sele:1977 ; Urata:1985 ; Sneyd:1994 ; Bojarevics:1994 ; Davidson:1998 ). The parameter is directly proportional to the product $J_{0}B_{0}$ and the square of the horizontal dimension of the cell, and inversely proportional to the product of the thicknesses of the metal and electrolyte layers, density difference $\Delta\rho$, and gravity acceleration. In our notation, for the instability developing at the upper interface, this becomes

[TABLE]

A similar parameter with $\Delta\rho^{B}$ instead of $\Delta\rho^{A}$ should be used at $\Delta\rho^{B}\ll\Delta\rho^{A}$ when the instability develops at the lower interface.

The expression (82) is equivalent to the parameter $\beta$ derived for the reduction cells in the pioneering work of Sele:1977 and to the parameter, in terms of which the results of Davidson:1998 were recently recast in Davidson:2016 . It can be rewritten in terms of the non-dimensional parameters introduced in section 2.4 as

[TABLE]

We have carried out extensive testing in order to verify that $\Pi$ is a suitable control parameter for the single-interface instability in the battery. Firstly, we verified that a significant change of $\Pi$ leads to a significant change of the wave characteristics $\gamma$ and $T$. Furthermore, tests were performed, in which the battery parameters, such as $L_{x}$, $B_{0}$, $J_{0}$, $H_{0}^{A}$, or $H_{0}^{E}$ were varied substantially in such a way that the value of $\Pi$ remained the same. The results of one of such tests conducted at $\Gamma=2.0$ are shown in table 1. Two sets of simulations are exhibited: at the normal values of viscosities and at the values reduced by four orders of magnitude, which practically removed the effect of viscosity on the flow.

We see in table 1 that at negligible viscosity the characteristics of the growing wave $\tilde{\gamma}$ and $\tilde{T}$ are determined by $\Pi$ with very good accuracy. The computed values are within 1%. An exception is the simulation with $\Delta\rho^{A}$ increased to 200 kg/m3, in which $\tilde{\gamma}$ and $\tilde{T}$ change by about 3%. This effect will be further discussed in section 3.2.

The variation is stronger at the realistic $\nu=5\times 10^{-7}$ m2/s. This can be attributed to the non-negligible effect of the variation of the Reynolds number.

We have also performed simulations with $\Delta\rho^{B}\ll\Delta\rho^{A}$, i.e. when the significant deformations are expected only on the lower interface. The results are qualitatively the same as when $\Delta\rho^{A}\ll\Delta\rho^{B}$. The only significant differences are the opposite sense of rotation (clockwise at $J_{0}>0$ and counterclockwise otherwise), which is evidently caused by the opposite direction of the horizontal current perturbations in the bottom metal layer (see (23) and the illustration in Fig. 1b), and the change of $T$ and $\gamma$ attributed to the change of the typical time scale.

Parametric studies were performed to explore the effect of $\Pi$ and the aspect ratio $\Gamma$ on the instability. The realistic values $\nu^{A}=\nu^{B}=\nu^{E}=5\times 10^{-7}$ m2/s were used. In the study, the same values $\rho^{A}=1000$ kg/m3, $\rho^{E}=1100$ kg/m3, $\rho^{B}=8000$ kg/m3, and $H_{0}^{A}=H_{0}^{B}=0.1$ m were used, while the other parameters were varied in the ranges $J_{0}=5\times 10^{3}$, $10^{4}$, A/m2 $B_{0}=10^{-3}$, $2\times 10^{-3}$ T, $0.5\leq L_{x}\leq 1.5$ m, and $0.005\leq H_{0}^{E}\leq 0.1$ m so as to change the value of $\Pi$. The results are presented in Fig. 4. We find that, as one would expect by analogy with the reduction cells, the behavior is strongly influenced by $\Gamma$. In particular, the instability threshold $\Pi_{cr}$ such that the battery is stable at all $\Pi<\Pi_{cr}$ but unstable at $\Pi>\Pi_{cr}$, is a non-monotonic function of $\Gamma$. At $\Gamma=1.0$, $\Pi_{cr}$ is very close to zero, although it cannot be exactly zero due to the stabilizing effect of viscosity and numerical dissipation. In rectangular cells, $\Pi_{cr}$ can be small, as, for example, $\Pi_{cr}(3.0)\approx 0.15$ or large as $\Pi_{cr}(2.0)\approx 3.7$ or $\Pi_{cr}(4.0)\approx 3.5$. Even small variations of $\Gamma$ may change $\Pi_{cr}$ quite significantly. An example is the difference between $\Pi_{cr}(3.0)\approx 0.15$ and $\Pi_{cr}(3.33)\approx 0.67$. This and the strong effect of $\Gamma$ in general can be explained in the same way as, e.g. in Davidson:1998 for the aluminum reduction cells. The aspect ratio determines the set of available natural gravitational wave modes and, so, strength of the electromagnetic effect needed to transform two of them into a pair with complex-conjugate eigenvalues.

As an illustration of the aspect ratio effect, Fig. 5 shows the interface deformations computed for the unstable flows at $\Gamma=3.0$ (Fig. 5a) and $\Gamma=3.33$ (Fig. 5a). Almost all the battery parameters are the same in the two cases: $J_{0}=10^{4}$ A/m2, $B_{0}=0.001$ T, $L_{x}=1.5$ m, $\rho^{A}=1000$ kg/m3, $\rho^{E}=1100$ kg/m3, $\rho^{B}=8000$ kg/m3, $H_{0}^{A}=H_{0}^{B}=0.1$ m. The only differences are in the values of $L_{y}$ (0.5 m or 0.45 m) and the values of $H_{0}^{E}$ (0.02548 m or 0.02317 m) that are selected in such a way that the instability parameter is $\Pi=3.0$ in both cases. The time signals of the two unstable solutions show strongly different characteristics: $T=0.690$, $\gamma=0.110$ at $\Gamma=3.0$ and $T=0.335$, $\gamma=0.073$ at $\Gamma=3.33$. The growing interface deformations shown in Fig. 5 demonstrate that the explanation suggested above is valid. The electromagnetic coupling, which causes the instability, involves different gravitational wave modes in the two cases, with much larger $x$-wavelength at $\Gamma=3.0$ than at $\Gamma=3.33$. A more specific description can be obtained with the eigenvalue analysis.

The main conclusion of this section is that the single-interface rolling pad instability occurring when one density difference is much smaller than the other is very similar to the instability observed in the aluminum reduction cells. The instability is caused by the same mechanism, develops in essentially the same form of a rotating interfacial wave, and is determined by three non-dimensional parameters: the instability parameter (82) (or its analog based on $\Delta\rho^{B}$, $H_{0}^{B}$ instead of $\Delta\rho^{A}$, $H_{0}^{A}$ if the instability occurs on the lower interface), the Reynolds number Re, and the horizontal aspect ratio $\Gamma$.

3.2 Double interface deformation

The situation when $\Delta\rho^{A}$ and $\Delta\rho^{B}$ are of comparable magnitudes, and, thus, the instability may significantly deform both the interfaces, is considered in this section. The analysis is performed as a parametric study, in which the metal densities are kept constant at $\rho^{A}=1000$ kg/m3 and $\rho^{B}=8000$ kg/m3, while the density of the electrolyte $\rho^{E}$ varies in the range between $1100$ and $6000$ kg/m3. This can be compared with the material densities in the currently considered battery schemes (see Kim:2013 ): 500 to 1500 kg/m3 for $\rho^{A}$, 6000 to 10000 kg/m3 for $\rho^{B}$, and 1000 to 3000 kg/m3 for $\rho^{E}$. In the study, we keep constant metal layer thicknesses $H^{A}=H^{B}=0.1$ m and viscosities $\nu^{A}=\nu^{B}=\nu^{E}=5\times 10^{-7}$ m2/s. The other parameters vary as $5\times 10^{3}\leq J_{0}\leq 4\times 10^{4}$ A/m2, $10^{-3}\leq B_{0}\leq 3\times 10^{-2}$ T, $1\leq L_{x}\leq 2$ m, $0.005\leq H^{E}\leq 0.04$ m. The behavior at the aspect ratios $\Gamma=2.0$ and 3.0 is explored.

As a typical example, Figs. 6-8 show the instability in a system with $\rho^{E}=3000$ kg/m3, $J_{0}=2\times 10^{4}$ A/m2, $B_{0}=0.002$ T, $L_{x}=2.0$ m, $L_{y}=1.0$ m, $H_{0}^{A}=H_{0}^{B}=0.1$ m, $H_{0}^{E}=0.005$ m. The results of the parametric studies are summarized in Fig. 9.

We see in Figs. 6-8 that the main features of the instability remain the same as in the case of the single-interface deformation. The time signals combine periodic oscillations with exponential growth. The spatial structure in Fig. 7 shows large-scale interfacial waves quite similar to the wave in Fig. 3. A conclusion can be made that the instability is still of the rolling pad type.

At the same time, the interaction between the two waves growing at the lower and upper interfaces results in more complex and diverse dynamics. The waves are always coupled to each other in the sense that they have the same oscillation period and sense of rotation (counterclockwise at $\Delta\rho^{A}<\Delta\rho^{B}$ and clockwise at $\Delta\rho^{A}>\Delta\rho^{B}$ if $J_{0}>0$). Two types of coupling are observed. At substantially different $\Delta\rho^{A}$ and $\Delta\rho^{B}$ (e.g. at $\rho^{E}\leq 3000$ kg/m3 or $\rho^{E}\geq 6000$ kg/m3 at $\Gamma=2.0$), the waves at the upper and lower interfaces are nearly antisymmetric (the sign of $\delta\zeta^{A}$ is opposite to the sign of $\delta\zeta^{B}$). The antisymmetry is not exact, but with a small phase shift. For example, the time signals of the single-point oscillations in Fig. 6a show the time shift $\Delta t\approx 0.042$, which can be compared with the oscillation period of 1.114. The situation changes at close $\Delta\rho^{A}$ and $\Delta\rho^{B}$. In such cases, the two waves become nearly symmetric (with the signs of $\delta\zeta^{A}$ and $\delta\zeta^{B}$ being mostly the same) and, again, with a small phase shift. This type of behavior is illustrated in Fig. 10.

The observed change of the type of wave coupling is consistent with the results of the recent three-dimensional potential flow analysis of interfacial waves in a three-layer cylindrical cell without electromagnetic forces Horstmann:2017 . It has been shown in this work that the type of coupling between the waves at the two interfaces changes from antisymmetric to symmetric when the amplitudes of the waves become comparable. While not a proof due to the evident differences between the systems, the consistency can be considered as an indication that the change of coupling observed in our battery is likely to be related to purely hydrodynamic interaction between the waves, rather than to electromagnetic forces.

The results of the parametric study presented in Fig. 9 show that the effect of the relation between $\Delta\rho^{A}$ and $\Delta\rho^{B}$ on stability is not very significant in the strongly unstable system with $\Gamma=3.0$. The values of the instability growth rate $\gamma$ decrease as the two density differences become closer to each other, but the instability continues to occur at small values of $\Pi$.

On the contrary, for the system with $\Gamma=2.0$, increase of $\Delta\rho^{A}$ has stabilizing effect, which is not only substantial, but also clearly exceeding the effect of the density stratification already incorporated into the expression (82) of the instability parameter. The critical values $\Pi_{cr}$, such that the system is always stable at $\Pi<\Pi_{cr}$ but unstable at $\Pi>\Pi_{cr}$, can be estimated as $\Pi_{cr}=3.7$, 4.0, 4.5, 6.2, and 21.5 at $\Delta\rho^{A}=100$, 500, 1000, 2000, and 3000 kg/m3, respectively. This growth is in contradiction with the prediction of the mechanical model Zikanov:2015 , where the presence of the second pendulum was found to be always destabilizing. We attribute the contradiction to the oversimplified character of the model used in Zikanov:2015 .

The stabilizing effect of the second layer can be given a qualitative explanation based on the fact that the horizontal perturbation currents in the top and bottom layers are always flowing in the opposite directions (see (23)). This implies the opposite directions of the horizontal Lorentz forces (24) and, thus, that forces in one layer always oppose the motion of a coupled wave system. Whether this happens in the bottom or top layer is determined by the relation between $\Delta\rho^{A}$, $\Delta\rho^{B}$, $H_{0}^{A}$, $H_{0}^{B}$ and, possibly, other parameters.

4 Concluding remarks

We have presented a new shallow water model and the results of its application to electromagnetically coupled waves in a simplified liquid metal battery of rectangular shape. The results primarily concern the linear stages of one particular instability in a simplified situation when the base state has flat interfaces and zero flow, but the model itself has much broader applicability. It is well suited for analysis of more general situations: with background melt flows and interface deformations, spatially complex or even time-dependent distributions of base electric current and magnetic field, etc.

The instability caused by the interaction between the externally generated vertical magnetic field and horizontal electric current perturbations associated with the interface deformations is demonstrated and analyzed. The growing perturbations have the form of rotating large-scale interfacial waves. The instability mechanism is similar to the mechanism of the rolling pad instability known for the aluminum reduction cells.

In the case when the density jump at one interface is much smaller than at the other, only one interface is significantly deformed in the course of the perturbation growth, and the similarity to the behavior of reduction cells is quite close. In particular, the instability is controlled by the same non-dimensional groups: the horizontal aspect ratio (78), the instability parameter (82), and the Reynolds number (76). Critically, the instability parameter (82) is proportional to the square of the horizontal size of the cell.

In the case when the density jumps at the two interfaces are comparable, both the interfaces are significantly deformed, and the behavior of the system is more complex and quite different from that of a reduction cell. The two interfacial waves can be coupled either symmetrically or antisymmetrically with a small phase shift between the waves. The presence of the second deformable interface may have a stabilizing effect.

The model used in this study is idealized and cannot produce a complete and quantitatively accurate picture of a real system. Further work is advised, preferably in the way of experiments performed with large battery cells and three-dimensional numerical models based on accurate representation of the geometry and the complex distributions of the base current and magnetic field. Eigenvalue analysis of the simplified 2D and 3D system would also be interesting.

Even at this point, however, the available results can be utilized to make a preliminary assessment of the potential effect of the rolling pad instability on battery design. We have seen above that the unstable perturbations tend to grow strongly resulting in high amplitudes of interface deformation or even rupture of the electrolyte layer. The situation is dissimilar to that predicted for the thermal convection or Tayler instability in the sense that reliable uninterrupted operation of a battery is only possible if the instability is avoided. Our conclusion is that the instability is a limiting factor for the versions of the battery design with small density difference between the metal and the electrolyte. One such version is the Mg-Sb battery Kim:2013 , in which $\rho^{A}\approx 1577$ kg/m3, $\rho^{E}\approx 1715$ kg/m3, $\rho^{B}\approx 6270$ kg/m3. Taking, as an example, the relatively stable geometry of a rectangular cell with $\Gamma=2.0$ and using the criterion $\Pi_{cr}=3.7$ derived in section 3.1 we find that at typical $J_{0}=10^{4}$ A/m2, $B_{0}=0.001$ T, a battery with $H_{0}^{A}=H_{0}^{B}=0.1$ m, $H_{0}^{E}=0.005$ m becomes unstable at the horizontal size $L_{x}$ exceeding 0.7077 m. This is larger than in the currently built laboratory prototypes or commercial concepts, but certainly in the range anticipated for future large-scale energy storage facilities Bradwell:2012 .

We should also mention that, similarly to the practice of the aluminum smelting industry, the instability can be avoided by rearranging the cell’s wiring so as to reduce the vertical component of the magnetic field within the cell. Computational simulation tools built on the basis of our model can be utilized to predict the need for such a rearrangement and to assist in its optimization.

The instability appears less important for other battery concepts, in which $\Delta\rho^{A}$ is substantially larger. The results presented in section 3.2 suggest strong stabilization effect. For example, increasing $\Delta\rho^{A}$ in the example above to realistically possible 2000 kg/m3 we find $L_{x}=3.5$ m, which does not appear to impose a serious limit on the technology.

Acknowledgements

The author is thankful to Norbert Weber, Tom Weier, Valdis Bojarevics, Janis Priede, and Gerrit Horstmann for interesting and stimulating discussions. Financial support was provided by the US NSF (Grant CBET 1435269).