2.1 Mathematical setting

We apply the potential flow analysis to an ideal cylinder, as depicted in figure 2. The cylindrical geometry was chosen because it is a generic case. A cylindrical cell is always unstable ( $\unicode[STIX]{x1D6FD}_{c}=0$ ) in the inviscid limit (Lukyanov et al. Reference Lukyanov, El and Molokov2001). Furthermore, there is no need to consider the influence of a horizontal aspect ratio on stability. Anyway, the theory described below is not restricted to cylinders but can be easily transferred to arbitrarily shaped vertical side walls. We show this on the example of a rectangular box in appendix A.

The cylindrical tank with the radius $R$ was defined to contain three immiscible phases $i=1,2,3$ of densities $\unicode[STIX]{x1D70C}_{i}$ and heights $h_{i}$ , where $\unicode[STIX]{x1D70C}_{1}<\unicode[STIX]{x1D70C}_{2}<\unicode[STIX]{x1D70C}_{3}$ has to be fulfilled to realize a stratified system. Two interfaces with interfacial tensions $\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}$ and $\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}$ are located at $z=\unicode[STIX]{x1D702}_{1}(x,y,t)$ and $z=\unicode[STIX]{x1D702}_{2}(x,y,t)$ . At this point we neglect interfacial tension effects at the side walls restricting the analysis to contact angles of $90^{\circ }$ and free-slip contact lines. The coordinate origin of the system is placed in the centre of the upper unperturbed interface $\unicode[STIX]{x1D702}_{1}$ .

Figure 2.

Sketch of the idealized cylindrical three-layer system determined by the radius $R$ , three phases $i=1,2,3$ of densities $\unicode[STIX]{x1D70C}_{i}$ and heights $h_{i}$ as well as two interfaces $\unicode[STIX]{x1D702}_{1}$ and $\unicode[STIX]{x1D702}_{2}$ with corresponding interfacial tensions $\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}$ and $\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}$ .

To facilitate potential flow theory, a few simplifications have to be made. First, an inviscid flow has to be assumed, by which damping effects at the tank walls and on the interfaces are neglected. Actually, viscosity affects the natural frequencies of sloshing interfaces only very slightly (Ibrahim Reference Ibrahim2005) and can be neglected without problems in comparison to the pressure exchange forces basically determining the coupling dynamics. Secondly, the flow is supposed to be irrotational. That seems to be an counterintuitive assumption at first glance since the MPRI is driven by rotational Lorentz forces. By using this approximation, we completely neglect the rotational part of the flow that causes and drives the instability. Nevertheless, as mentioned before, the induced wave dynamics in the unstable regime can be well approximated by the irrotational part of the flow only. In § 3.2.1 it is shown that rotational Lorentz forces primarily influence the shape of the waves, but not their frequencies and amplitude ratios. Finally, the wave amplitudes are restricted to be sufficiently small to facilitate linear wave theory and all three fluids are treated as incompressible.

2.2 Potential flow analysis

Under the above simplifications the velocity fields $\boldsymbol{u}_{i}$ of all three phases $i=1,2,3$ can be expressed unambiguously by the gradient of scalar potentials $\unicode[STIX]{x1D719}_{i}$ as

(2.1) $$\begin{eqnarray}\boldsymbol{u}_{i}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}_{i}.\end{eqnarray}$$

For fulfilling momentum and mass conservation, Bernoulli’s and Laplace’s equation have to be solved in all phases

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{i}}{\unicode[STIX]{x2202}t}}+{\displaystyle \frac{P_{i}}{\unicode[STIX]{x1D70C}_{i}}}+gz=c_{i}(t), & \displaystyle\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D719}_{i}=0, & \displaystyle\end{eqnarray}$$

where $P_{i}$ denotes the pressure, $g$ the standard acceleration due to gravity and $c_{i}(t)$ are arbitrary time-dependent integration constants. In addition, using cylindrical coordinates ( $r,\unicode[STIX]{x1D711},z$ ), the following sets of dynamic boundary conditions must be satisfied:

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \text{top wall:}\quad {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z}}=0|_{z=h_{1}}, & \displaystyle\end{eqnarray}$$

(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \text{bottom wall:}\quad {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}z}}=0|_{z=-(h_{2}+h_{3})}, & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \displaystyle \text{side wall:}\quad {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}r}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}r}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}r}}=0|_{r=R}, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle \text{upper interface:}\quad {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z}}\bigg|_{z=\unicode[STIX]{x1D702}_{1}}, & \displaystyle\end{eqnarray}$$

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle \text{lower interface:}\quad {\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}z}}\bigg|_{z=\unicode[STIX]{x1D702}_{2}}. & \displaystyle\end{eqnarray}$$

No-outflow boundary conditions were formulated at the side wall $r=R$ (2.6), at the top $z=h_{1}$  (2.4) and the bottom wall $z=-(h_{2}+h_{3})$  (2.5), while conditions (2.7) and (2.8) ensure interface preservation. By assuming a harmonic time dependency $\unicode[STIX]{x1D719}_{i}\sim \cos (\unicode[STIX]{x1D714}t)$ , where $\unicode[STIX]{x1D714}$ marks the angular frequency, potential solutions for the Laplace equations (2.3) in cylindrical coordinates fulfilling all above conditions can be found as superpositions of infinite wave modes $m\in \mathbb{N}_{0}$ and $n\in \mathbb{N}_{1}$ :

(2.9a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{1}(r,\unicode[STIX]{x1D711},z,t) & = & \displaystyle -\mathop{\sum }_{m=0}^{\infty }\mathop{\sum }_{n=1}^{\infty }\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}{\displaystyle \frac{\cosh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}(z-h_{1})\right)}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}(h_{1})\right)}}\text{J}_{m}\left(\unicode[STIX]{x1D716}_{mn}{\displaystyle \frac{r}{R}}\right)\cos (m\unicode[STIX]{x1D711}-\unicode[STIX]{x1D714}t),\end{eqnarray}$$

(2.9b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{2}(r,\unicode[STIX]{x1D711},z,t) & = & \displaystyle \mathop{\sum }_{m=0}^{\infty }\mathop{\sum }_{n=1}^{\infty }\left[\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}\left(\text{e}^{{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}}\cosh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}z\right)-\text{e}^{-{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}z}\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)\right)\right.\nonumber\\ \displaystyle & & \displaystyle \left.-\,\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}\cosh \vphantom{\frac{\mathop{\sum }_{\int }}{\mathop{\sum }_{\int }}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}z\right)\right]\cdot {\displaystyle \frac{1}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)}}\text{J}_{m}\left(\unicode[STIX]{x1D716}_{mn}{\displaystyle \frac{r}{R}}\right)\cos (m\unicode[STIX]{x1D711}-\unicode[STIX]{x1D714}t),\qquad\end{eqnarray}$$

(2.9c ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D719}_{3}(r,\unicode[STIX]{x1D711},z,t) & = & \displaystyle \mathop{\sum }_{m=0}^{\infty }\mathop{\sum }_{n=1}^{\infty }\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}{\displaystyle \frac{\cosh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}(z+h_{2}+h_{3})\right)}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}(h_{3})\right)}}\text{J}_{m}\left(\unicode[STIX]{x1D716}_{mn}{\displaystyle \frac{r}{R}}\right)\cos (m\unicode[STIX]{x1D711}-\unicode[STIX]{x1D714}t).\qquad\end{eqnarray}$$

$\text{J}_{m}$

$m$

$\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}$

$\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}$

$m$

$m=0$

$\unicode[STIX]{x1D716}_{mn}$

$n$

$m$

(2.10) $$\begin{eqnarray}\text{J}_{m}^{\prime }(\unicode[STIX]{x1D716}_{mn})=0,\end{eqnarray}$$

to fulfil the radial no-outflow boundary condition. The roots can be easily determined numerically or found in, e.g. Abramowitz & Stegun (Reference Abramowitz and Stegun1972). Exploiting the fact that the vertical velocity on the interfaces is identical to its linearized time variation

(2.11a,b ) $$\begin{eqnarray}\left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}_{1}}{\unicode[STIX]{x2202}t}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z}}\right|_{z=0}\quad \text{and}\quad \left.{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}_{2}}{\unicode[STIX]{x2202}t}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}z}}\right|_{z=-h_{2}},\end{eqnarray}$$

the interfacial shapes can be easily derived from the scalar flow potentials as:

(2.12a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{1}(r,\unicode[STIX]{x1D711},t) & = & \displaystyle \mathop{\sum }_{m=0}^{\infty }\mathop{\sum }_{n=1}^{\infty }\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R\unicode[STIX]{x1D714}}}\text{J}_{m}\left(\unicode[STIX]{x1D716}_{mn}{\displaystyle \frac{r}{R}}\right)\cos (m\unicode[STIX]{x1D711}-\unicode[STIX]{x1D714}t),\end{eqnarray}$$

(2.12b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D702}_{2}(r,\unicode[STIX]{x1D711},t) & = & \displaystyle \mathop{\sum }_{m=0}^{\infty }\mathop{\sum }_{n=1}^{\infty }\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R\unicode[STIX]{x1D714}}}\text{J}_{m}\left(\unicode[STIX]{x1D716}_{mn}{\displaystyle \frac{r}{R}}\right)\cos (m\unicode[STIX]{x1D711}-\unicode[STIX]{x1D714}t).\end{eqnarray}$$

$\text{Mg}\Vert \text{Sb}$

$m=n=1$

$\unicode[STIX]{x1D716}_{11}\approx 1.841$

$\unicode[STIX]{x0394}P$

$\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}$

$\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}$

(2.13a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}P|_{\unicode[STIX]{x1D702}_{1}} & = & \displaystyle P_{2}|_{\unicode[STIX]{x1D702}_{1}}-P_{1}|_{\unicode[STIX]{x1D702}_{1}}=-\unicode[STIX]{x1D70C}_{2}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}t}}\bigg|_{\unicode[STIX]{x1D702}_{1}}-\unicode[STIX]{x1D70C}_{2}g\unicode[STIX]{x1D702}_{1}+\unicode[STIX]{x1D70C}_{1}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}t}}\bigg|_{\unicode[STIX]{x1D702}_{1}}+\unicode[STIX]{x1D70C}_{1}g\unicode[STIX]{x1D702}_{1},\end{eqnarray}$$

(2.13b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}P|_{\unicode[STIX]{x1D702}_{2}} & = & \displaystyle P_{3}|_{\unicode[STIX]{x1D702}_{2}}-P_{2}|_{\unicode[STIX]{x1D702}_{2}}=-\unicode[STIX]{x1D70C}_{3}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}t}}\bigg|_{\unicode[STIX]{x1D702}_{2}}-\unicode[STIX]{x1D70C}_{3}g\unicode[STIX]{x1D702}_{2}+\unicode[STIX]{x1D70C}_{2}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}t}}\bigg|_{\unicode[STIX]{x1D702}_{2}}+\unicode[STIX]{x1D70C}_{2}g\unicode[STIX]{x1D702}_{2}.\end{eqnarray}$$

$c_{i}(t)$

$\unicode[STIX]{x1D719}_{i}$

$R_{1}$

$R_{2}$

(2.14a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}P|_{\unicode[STIX]{x1D702}_{1}} & = & \displaystyle \unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}\left({\displaystyle \frac{1}{R_{1}}}+{\displaystyle \frac{1}{R_{2}}}\right)\bigg|_{\unicode[STIX]{x1D702}_{1}},\end{eqnarray}$$

(2.14b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}P|_{\unicode[STIX]{x1D702}_{2}} & = & \displaystyle \unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\left({\displaystyle \frac{1}{R_{1}}}+{\displaystyle \frac{1}{R_{2}}}\right)\bigg|_{\unicode[STIX]{x1D702}_{2}}.\end{eqnarray}$$

(2.15a,b ) $$\begin{eqnarray}{\displaystyle \frac{1}{R_{1}}}\bigg|_{\unicode[STIX]{x1D702}_{j}}\approx {\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D702}_{j}}{\unicode[STIX]{x2202}x^{2}}},\quad {\displaystyle \frac{1}{R_{2}}}\bigg|_{\unicode[STIX]{x1D702}_{j}}\approx {\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D702}_{j}}{\unicode[STIX]{x2202}y^{2}}}\quad \text{with }j=1,2.\end{eqnarray}$$

By now, the dynamic boundary condition can be expressed in dependence on the flow potentials $\unicode[STIX]{x1D719}_{i}$ by calculating the first time derivatives of (2.13a ), exploiting the preservation conditions (2.11) and eliminating the pressure drops. One gets

(2.16a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \left.\unicode[STIX]{x1D70C}_{1}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}t^{2}}}+\unicode[STIX]{x1D70C}_{1}g{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z}}=\unicode[STIX]{x1D70C}_{2}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}t^{2}}}+\unicode[STIX]{x1D70C}_{2}g{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z}}-\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}\left({\displaystyle \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}}+{\displaystyle \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}}\right){\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z}}\right|_{z=\unicode[STIX]{x1D702}_{1}},\quad\end{eqnarray}$$

(2.16b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \left.\unicode[STIX]{x1D70C}_{2}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}t^{2}}}+\unicode[STIX]{x1D70C}_{2}g{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z}}=\unicode[STIX]{x1D70C}_{3}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}t^{2}}}+\unicode[STIX]{x1D70C}_{3}g{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}z}}-\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\left({\displaystyle \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}x^{2}}}+{\displaystyle \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}y^{2}}}\right){\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}z}}\right|_{z=\unicode[STIX]{x1D702}_{2}}.\end{eqnarray}$$

(2.17) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{i}}{\unicode[STIX]{x2202}z^{2}}}=-\left({\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{i}}{\unicode[STIX]{x2202}x^{2}}}+{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{i}}{\unicode[STIX]{x2202}y^{2}}}\right)\end{eqnarray}$$

so that the kinematic boundary conditions finally read

(2.18a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \left.\unicode[STIX]{x1D70C}_{1}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}t^{2}}}+\unicode[STIX]{x1D70C}_{1}g{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{1}}{\unicode[STIX]{x2202}z}}=\unicode[STIX]{x1D70C}_{2}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}t^{2}}}+\unicode[STIX]{x1D70C}_{2}g{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}{\displaystyle \frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z^{3}}}\right|_{z=\unicode[STIX]{x1D702}_{1}},\end{eqnarray}$$

(2.18b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \left.\unicode[STIX]{x1D70C}_{2}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}t^{2}}}+\unicode[STIX]{x1D70C}_{2}g{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{2}}{\unicode[STIX]{x2202}z}}=\unicode[STIX]{x1D70C}_{3}{\displaystyle \frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}t^{2}}}+\unicode[STIX]{x1D70C}_{3}g{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}z}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}{\displaystyle \frac{\unicode[STIX]{x2202}^{3}\unicode[STIX]{x1D719}_{3}}{\unicode[STIX]{x2202}z^{3}}}\right|_{z=\unicode[STIX]{x1D702}_{2}}.\end{eqnarray}$$

$x$

$y$

$\unicode[STIX]{x1D714}(\unicode[STIX]{x1D716}_{mn})$

(2.19a ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{1|2mn}^{2} & = & \displaystyle {\displaystyle \frac{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}}{\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\left(\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)-{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}}{\displaystyle \frac{1}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)}}\right)}},\end{eqnarray}$$

(2.19b ) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D714}_{2|3mn}^{2} & = & \displaystyle {\displaystyle \frac{(\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}}{\unicode[STIX]{x1D70C}_{3}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{3}\right)+\unicode[STIX]{x1D70C}_{2}\left(\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)-{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}{\displaystyle \frac{1}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)}}\right)}}.\end{eqnarray}$$

$\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}/\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}$

$\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}/\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}$

(2.20) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}a\unicode[STIX]{x1D714}^{4}+b\unicode[STIX]{x1D714}^{2}+c=0,\quad \text{with}\\ a=\left(\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)+\unicode[STIX]{x1D70C}_{3}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{3}\right)\right)\left(\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)+\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{1}\right)\right)\\ \quad \quad -{\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}^{2}}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)^{2}}},\\ b=-\left((\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}\right)\left(\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)+\unicode[STIX]{x1D70C}_{3}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{3}\right)\right)\\ \quad \quad -\left((\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}\right)\left(\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)+\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{1}\right)\right),\\ c=\left((\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}\right)\left((\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}\right),\end{array}\right\} & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

which provides two independent analytical solutions

(2.21) $$\begin{eqnarray}\unicode[STIX]{x1D714}_{\pm }^{2}={\displaystyle \frac{-b\pm \sqrt{b^{2}-4ac}}{2a}}.\end{eqnarray}$$

Consequently, in contrast to two-layer ARCs, LMBs involve always two eigenfrequencies $\unicode[STIX]{x1D714}_{+}$ and $\unicode[STIX]{x1D714}_{-}$ . It can be seen readily from (2.21) that $\unicode[STIX]{x1D714}_{+}$ is always larger than $\unicode[STIX]{x1D714}_{-}$ . Hence, in analogy to the coupled harmonic oscillator, $\unicode[STIX]{x1D714}_{+}$ is called the fast mode and $\unicode[STIX]{x1D714}_{-}$ the slow mode, respectively. Note that Weber et al. (Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier2017b ) gave the dispersion relation (2.21) for the specific case of $\text{Mg}\Vert \text{Sb}$ cells.

Figure 3.

Schematic of the symmetric fast mode $\unicode[STIX]{x1D714}_{+}$ (red) and the antisymmetric slow mode $\unicode[STIX]{x1D714}_{-}$ (blue) of the wavenumber $\unicode[STIX]{x1D716}_{11}$ in cross-section.

Moreover, by knowing the frequency solutions from (2.21), the amplitude ratios can be analytically determined by transposing (2.19a ) or (2.19b ) giving

(2.22) $$\begin{eqnarray}\displaystyle {\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}_{\pm } & = & \displaystyle {\displaystyle \frac{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)}{\unicode[STIX]{x1D70C}_{2}}}\left[\vphantom{\frac{\mathop{\sum }_{\int }}{\mathop{\sum }_{\int }^{2}}}\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)+\unicode[STIX]{x1D70C}_{3}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{3}\right)\right.\nonumber\\ \displaystyle & & \displaystyle -\,\left.{\displaystyle \frac{(\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}}{\unicode[STIX]{x1D714}_{\pm }^{2}}}\right].\end{eqnarray}$$

Corresponding to the two frequency solutions, we also find two different amplitude ratios, whereby it can readily be shown from (2.22) that for the plus mode the amplitude ratio is always positive and for the minus mode always negative

(2.23a,b ) $$\begin{eqnarray}{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}_{+}>0,\quad {\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}_{-}<0.\end{eqnarray}$$

From that directly follows that in the plus mode the interfaces propagate symmetrically (in phase) and in the minus mode antisymmetrically (with a $180^{\circ }$ phase shift). Exemplarily, both coupling states are schematically sketched in figure 3.

To conclude, we have to distinguish between two essential coupling modes, the one fast and symmetric and the other slow and antisymmetric. Both modes may occur in LMBs having different consequences for the operation stability. This result is well known from the coupled harmonic oscillator and was to be reproduced since we employed linear wave theory. However, the formulas (2.21) and (2.22) now allow us to study the natural frequencies and the coupling dynamics in dependence on nine geometrical cell parameters and material properties $\unicode[STIX]{x1D70C}_{1},\unicode[STIX]{x1D70C}_{2},\unicode[STIX]{x1D70C}_{3},\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}},\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}},h_{1},h_{2},h_{3},R$ as presented in the next section.

2.3 Analysis of the coupling dynamics

In the following we analyse the strength of the interfacial wave coupling to predict in which parameter regimes the battery can be considered as decoupled and two-layer analysis to be sufficient to understand the MPRI, and in which regimes the wave coupling has a large share in the instability mechanism.

At fist we compare the expressions (2.19a ) and (2.19b ) of the three-layer dispersion relation with the corresponding two-layer dispersion relations $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{1}}$ and $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{2}}$ of both interfaces $\unicode[STIX]{x1D702}_{1}$ and $\unicode[STIX]{x1D702}_{2}$ that can be analogously derived and are given by

(2.24a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{1}}^{2}={\displaystyle \frac{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}}{\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)}},\end{eqnarray}$$

(2.24b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{2}}^{2}={\displaystyle \frac{(\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}}{\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)+\unicode[STIX]{x1D70C}_{3}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{3}\right)}}.\end{eqnarray}$$

(2.25a,b ) $$\begin{eqnarray}{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}}{\displaystyle \frac{1}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)}}\quad \text{or}\quad {\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}{\displaystyle \frac{1}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)}},\end{eqnarray}$$

respectively. These terms can be exploited to make first estimations in which cases the battery can be considered as decoupled. The terms are vanishing for two different limits: on the one hand if the salt layer height $h_{2}$ is becoming large in comparison to the radius $R$ . In that limit both interfaces are placed too far apart from each other to interact and are both propagating in their own two-layer eigenfrequency $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{1}}$ and $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{2}}$ . Effectively, this limit is not practically relevant for LMBs since the salt layer should always be as thin as possible to realize maximum efficiency. More important is the second limit, where the amplitude ratio is becoming sufficiently small or large

(2.26a,b ) $$\begin{eqnarray}{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}\gg 1\quad \text{or}\quad {\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}\ll 1.\end{eqnarray}$$

In these limits one of the frequencies decreases asymptotically to zero and the other one reaches the two-layer limit. For which parameters this will be the case is principally determined by (2.22). However, in order to understand the key mechanisms of the interfacial coupling some deeper analyses are required.

Figure 4.

(a) Three-layer dispersion relations due to (2.21) and two-layer dispersion relations due to (2.24a ) and (2.24b ) of the $\text{Mg}\Vert \text{Sb}$ cell (table 1) for different salt-layer heights $h_{2}=0.15,0.2,0.3,0.5$ and $5\,\text{ cm}$ . The plus mode $\unicode[STIX]{x1D714}_{+}$ is red and the minus mode $\unicode[STIX]{x1D714}_{-}$ blue coloured. The solid black line shows the two-layer frequency of the upper interface $\unicode[STIX]{x1D702}_{1}$ and the dashed black line the two-layer frequency of the lower interface $\unicode[STIX]{x1D702}_{2}$ both for $h_{2}=5\,\text{ cm}$ . (b) Absolute value of the amplitude ratios due to (2.22) in dependence on $\unicode[STIX]{x1D716}_{mn}$ for $h_{2}=0.1,1.1,2.2,3.1,4.1,5.1,6.1$ and $7.1\,\text{ cm}$ .

Table 1.

Cell parameters of the $\text{Mg}\Vert \text{Sb}$ cell from Weber et al. (Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier2017b ) used to calculate the wave frequencies and amplitude ratios of figure 4.

In the first instance, to give the reader a better understanding, figure 4 shows exemplarily the frequencies (a) due to (2.21) as well as the amplitude ratios (b) due to (2.22). Both are given for varying electrolyte-layer heights in dependency of the wavenumber $\unicode[STIX]{x1D716}_{mn}$ for the same $\text{Mg}\Vert \text{Sb}$ cell numerically investigated by Weber et al. (Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier2017b ). The cell parameters and material properties used are shown in table 1. Figure 4(a) clarifies that both frequency modes are progressively approaching the two-layer limits for increasing salt-layer heights. Another remarkable feature is the frequency transition appearing at $\unicode[STIX]{x1D716}_{mn}\approx 17$ , where the two-layer frequencies are crossing each other. The transition is emphasized in the zoomed window of figure 4(a) showing the two-layer frequencies $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{1}}$ and $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{2}}$ as well as the decoupled three-layer frequencies $\unicode[STIX]{x1D714}_{+}$ and $\unicode[STIX]{x1D714}_{-}$ for large salt layers. It can be seen that for small wavenumbers the fast mode $\unicode[STIX]{x1D714}_{+}$ is following the two-layer eigenfrequency $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{2}}$ of the lower interface while the slow mode $\unicode[STIX]{x1D714}_{-}$ oscillates with the two-layer frequency $\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{1}}$ of the upper interface. From $\unicode[STIX]{x1D716}_{mn}\approx 17$ onwards both modes behave in precisely the opposite manner. Hence, both modes cannot be uniquely assigned to the natural frequencies of uncoupled interfaces. Such a behaviour has not been observed yet in coupled internal and surface waves (Issenmann et al. Reference Issenmann, Laroche and Falcon2016) and leads to some new physical properties. Figure 4(b) shows the absolute values of the amplitude ratio clarifying the strength of coupling in dependence on $\unicode[STIX]{x1D716}_{mn}$ and $h_{2}$ . If we for instance focus on the smallest wavenumber $\unicode[STIX]{x1D716}_{11}$ representing the metal pad roll, one can see that for the $\text{Mg}\Vert \text{Sb}$ cell the amplitude ratios of both modes are not symmetric around $|\tilde{\unicode[STIX]{x1D702}}_{1}^{11}/\tilde{\unicode[STIX]{x1D702}}_{2}^{11}|=1$ . In this case, oscillations of the upper interface will influence the lower one only very slightly. The slow mode $\unicode[STIX]{x1D714}_{-}$ will be excited at the upper interface since it is close to its eigenfrequency. Further, the amplitude ratio is very large with $|\tilde{\unicode[STIX]{x1D702}}_{1}^{11}/\tilde{\unicode[STIX]{x1D702}}_{2}^{11}|_{-}\approx 38$ such that $\tilde{\unicode[STIX]{x1D702}}_{2}^{11}$ becomes negligibly small for reasonable values of $\tilde{\unicode[STIX]{x1D702}}_{1}^{11}$ . In contrast, perturbations of the lower interface with an eigenfrequency close to $\unicode[STIX]{x1D714}_{+}$ will substantially affect the upper one because there the amplitudes are strongly coupled $|\tilde{\unicode[STIX]{x1D702}}_{1}^{11}/\tilde{\unicode[STIX]{x1D702}}_{2}^{11}|_{+}\approx 0.7$ . The above mentioned transition can be seen as well in the amplitude ratio curves at $\unicode[STIX]{x1D716}_{mn}\approx 17$ . There, all curves intersect with $|\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}/\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}|_{\pm }\approx 1$ and both waves are perfectly coupled by both modes. That clarifies that there is not only an exchange of the frequencies, but also of the modes within the interfaces.

Equations (2.21) together with (2.22) are principally sufficient to completely describe the gravitational–capillary coupling dynamics. However, it is worthwhile also to analytically predict the intersection point and to understand the physical nature of this transition. While the presented transition as a function of the wavenumber is not of practical importance for LMBs, we have observed this transition also in dependence on many further cell parameters significantly changing for the various possible liquid metal and electrolyte combinations. We discovered that the influence of all cell parameters can be unified by only two dimensionless parameters. We can prove that by eliminating the frequencies in (2.19a ) and (2.19b ). This leads to a quadratic equation for the amplitude ratio:

(2.27) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}l@{}}a\left({\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}\right)^{2}+b{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}-c=0,\quad \text{with}\\ a={\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)\left(\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)\right)}},\\ b=\left({\displaystyle \frac{(\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}}{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}}}-{\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)+\unicode[STIX]{x1D70C}_{3}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{3}\right)}{\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)}}\right),\\ c={\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}}{\sinh \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)\left(\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)\right)}}\\ \qquad \times {\displaystyle \frac{(\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}}{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g{\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}\right)^{3}}}.\end{array}\right\}\qquad & & \displaystyle\end{eqnarray}$$

As clarified in figure 4(b), the transition always appears at the intersection point of the modes, where the absolute values of the amplitude ratios have to equalize

(2.28) $$\begin{eqnarray}\left|{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}\right|_{+}\stackrel{!}{=}\left|{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}\right|_{-}.\end{eqnarray}$$

That is automatically fulfilled if (2.27) has only one solution, which in turn is guaranteed only if $b$ vanishes. This is again fulfilled only for

(2.29) $$\begin{eqnarray}{\mathcal{A}}\stackrel{!}{=}{\mathcal{B}},\end{eqnarray}$$

where the two dimensionless numbers ${\mathcal{A}}$ and ${\mathcal{B}}$ are defined with

(2.30) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{A}}:={\displaystyle \frac{(\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})gR^{2}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\unicode[STIX]{x1D716}_{mn}^{2}}{(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})gR^{2}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}\unicode[STIX]{x1D716}_{mn}^{2}}}, & \displaystyle\end{eqnarray}$$

(2.31) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{B}}:={\displaystyle \frac{\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)+\unicode[STIX]{x1D70C}_{3}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{3}\right)}{\unicode[STIX]{x1D70C}_{1}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{1}\right)+\unicode[STIX]{x1D70C}_{2}\coth \left({\displaystyle \frac{\unicode[STIX]{x1D716}_{mn}}{R}}h_{2}\right)}}. & \displaystyle\end{eqnarray}$$

Hence, the parameters ${\mathcal{A}}$ and ${\mathcal{B}}$ completely determine the transition. ${\mathcal{A}}$ may be interpreted as the ratio of the restoring gravity and capillary forces acting on both waves, whereas ${\mathcal{B}}$ describes the inertia distribution in the battery system. However, it can be shown that in practical cases the transition process is mainly captured by ${\mathcal{A}}$ alone. For the desired thin electrolyte-layer limit $h_{2}\rightarrow 0$ , ${\mathcal{B}}$ converges asymptotically to 1 and (2.27) simplifies to

(2.32) $$\begin{eqnarray}\left({\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}\right)^{2}+({\mathcal{A}}-1){\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}-{\mathcal{A}}=0.\end{eqnarray}$$

In this limit the amplitude ratio has become independent of ${\mathcal{B}}$ and offers two analytical solutions

(2.33a,b ) $$\begin{eqnarray}{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}_{+}=1\quad \text{and}\quad {\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}_{-}=-{\mathcal{A}}.\end{eqnarray}$$

From this we note that for thin salt layers the interfaces are perfectly coupled for all parameters when they are excited by the fast mode $\unicode[STIX]{x1D714}_{+}$ . In contrast, the amplitude ratio of the slow mode $\unicode[STIX]{x1D714}_{-}$ is directly determined by ${\mathcal{A}}$ .

Finally, ${\mathcal{A}}$ and ${\mathcal{B}}$ can be further exploited to physically understand the mode transition. At first we consider the limits of pure gravity and pure capillary waves. In these limits, ${\mathcal{A}}$ simplifies to

(2.34a,b ) $$\begin{eqnarray}{\mathcal{A}}_{gravity}={\displaystyle \frac{\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2}}{\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1}}}\quad \text{or}\quad {\mathcal{A}}_{capillary}={\displaystyle \frac{\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}}{\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}}},\end{eqnarray}$$

respectively. Hence, in large-scale LMBs the wave coupling is determined only by the ratio of the density differences. ${\mathcal{A}}_{gravity}$ is the most important control parameter that should always be considered for predicting resonance. Further, the wavenumber-dependent mode transition, as shown in figure 4, can be explained considering both limits: for the $\text{Mg}\Vert \text{Sb}$ cell we find ${\mathcal{A}}_{gravity}>1$ but ${\mathcal{A}}_{capillary}<1$ . Consequently, this mode transition arose due to the change from gravity to capillary waves turning the sign of ${\mathcal{A}}$ .

The coupling mechanism can be further understood by noting that ${\mathcal{A}}/{\mathcal{B}}$ can be expressed by the ratio of the two-layer frequencies

(2.35) $$\begin{eqnarray}{\displaystyle \frac{{\mathcal{A}}}{{\mathcal{B}}}}={\displaystyle \frac{\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{2}}^{2}}{\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{1}}^{2}}}.\end{eqnarray}$$

Figure 5.

Absolute values of the amplitude ratios due to (2.22) for the metal pad roll mode $\unicode[STIX]{x1D716}_{11}$ in dependence on ${\mathcal{A}}/{\mathcal{B}}$ with varying salt-layer thicknesses $h_{2}=0.1,1,2,3,4,5,6$ and $7~\text{cm}$ exemplarily plotted using the generic $\text{Mg}\Vert \text{Sb}$ cell parameters from Weber et al. (Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier2017b ) (table 1). To diversify ${\mathcal{A}}$ and ${\mathcal{B}}$ , $\unicode[STIX]{x1D70C}_{3}$ was changed from $1716$ up to $5000~\text{kg}~\text{m}^{-3}$ .

Thus, as one would intuitively expect, the highest coupling indeed arises where the two-layer eigenfrequencies coincide and the interfaces step into resonance. The ratio ${\mathcal{A}}/{\mathcal{B}}$ quantifies the difference of the natural frequencies of both interfaces and therefore determines the strength of the coupling. For the important thin salt-layer limit $h_{2}\rightarrow 0$ , we further obtain

(2.36) $$\begin{eqnarray}\left|{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}\right|_{-}={\mathcal{A}}={\displaystyle \frac{\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{2}}^{2}}{\unicode[STIX]{x1D714}_{\unicode[STIX]{x1D702}_{1}}^{2}}},\end{eqnarray}$$

showing that the amplitude ratio corresponding to the slow mode $\unicode[STIX]{x1D714}_{-}$ is directly determined by the ratio of the two-layer frequencies.

All these results are exemplified in figure 5 showing the absolute values of the amplitude ratios in dependence on ${\mathcal{A}}/{\mathcal{B}}$ for different salt-layer thicknesses of the $\text{Mg}\Vert \text{Sb}$ cell. All curves are indeed crossing at ${\mathcal{A}}={\mathcal{B}}$ . The two inner lines for the smallest chosen $h_{2}$ correspond to the solutions (2.33) bounding the maximum coupling state. For large ${\mathcal{A}}/{\mathcal{B}}\gg 1$ the antisymmetric mode $\unicode[STIX]{x1D714}_{-}$ dominates in the upper interface $\unicode[STIX]{x1D702}_{1}$ , while for low values ${\mathcal{A}}/{\mathcal{B}}\ll 1$ this mode is more prominent in the lower one $\unicode[STIX]{x1D702}_{2}$ . The coupling state of the symmetric mode $\unicode[STIX]{x1D714}_{+}$ in contrast is not limited by ${\mathcal{A}}/{\mathcal{B}}$ since for sufficiently small $h_{2}$ the interfaces are always perfectly coupled by that mode.

2.4 Coupling criterion for LMBs

These findings now can be exploited to determine for which parameter regimes coupled dynamics must be considered and in which cases two-layer stability analysis, as practiced e.g. in the context of ARCs, is sufficient to understand interfacial instabilities in LMBs. At first it is desirable to predict which coupling mode will be excited by the arising Lorentz forces. As presented in the introduction, in two-layer systems the MPRI can be appropriately described by the Sele criterion. Since interfacial tension was respected for the coupling theory, it also has to be included into the Sele criterion. Following Gerbeau et al. (Reference Gerbeau, Le Bris and Lelièvre2006), $\unicode[STIX]{x1D6FD}$ can be regarded as the ratio of the electromagnetic force due to interface displacements to the gravity force exerted on the interface due to the same displacements. By considering the sum of gravity and capillarity as the acting restoring force, we get modified Sele criteria for both interfaces

(2.37) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D702}_{1}}={\displaystyle \frac{I_{0}b_{z}}{\left[(\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})g+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{11}}{R}}\right)^{2}\right]h_{1}h_{2}}}>\unicode[STIX]{x1D6FD}_{crit}, & \displaystyle\end{eqnarray}$$

(2.38) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D702}_{2}}={\displaystyle \frac{I_{0}b_{z}}{\left[(\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})g+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\left({\displaystyle \frac{\unicode[STIX]{x1D716}_{11}}{R}}\right)^{2}\right]h_{2}h_{3}}}>\unicode[STIX]{x1D6FD}_{crit}. & \displaystyle\end{eqnarray}$$

From these criteria we can easily identify which interface is becoming unstable first. For example, for large-size LMBs it is always the interface with the lower density difference. The ratio $\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D702}_{2}}/\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D702}_{1}}$ can be directly expressed in dependence of ${\mathcal{A}}$

(2.39) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D702}_{2}}}{\unicode[STIX]{x1D6FD}_{\unicode[STIX]{x1D702}_{1}}}}={\displaystyle \frac{((\unicode[STIX]{x1D70C}_{2}-\unicode[STIX]{x1D70C}_{1})gR^{2}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}\unicode[STIX]{x1D716}_{11}^{2})h_{1}}{((\unicode[STIX]{x1D70C}_{3}-\unicode[STIX]{x1D70C}_{2})gR^{2}+\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}\unicode[STIX]{x1D716}_{11}^{2})h_{3}}}={\displaystyle \frac{1}{{\mathcal{A}}}}{\displaystyle \frac{h_{1}}{h_{3}}}\end{eqnarray}$$

showing that for large ${\mathcal{A}}$ the upper interface will be more unstable while for small ${\mathcal{A}}$ the lower interface becomes more susceptible to electromagnetic excitation. By considering that at the first moment the natural two-layer frequencies will be excited until the third layer is interacting after some response time, it can be clarified that only the slow mode $\unicode[STIX]{x1D714}_{-}$ will be excited in LMBs for small and large values of ${\mathcal{A}}$ because for large ${\mathcal{A}}$ the natural two-layer frequency of the first excited upper interface is slower and close to $\unicode[STIX]{x1D714}_{-}$ (see (2.19) or figure 4 a) while for small ${\mathcal{A}}$ the lower interface becomes slower and even close to $\unicode[STIX]{x1D714}_{-}$ . In both cases only the antisymmetric mode $\unicode[STIX]{x1D714}_{-}$ develops such that LMBs can principally decouple. On that basis a coupling criterion can be deduced. Note, that for the limit $h_{2}\rightarrow 0$ the coupling is always on its highest level. When neglecting $\unicode[STIX]{x1D714}_{+}$ , we find from (2.33) that

(2.40a,b ) $$\begin{eqnarray}\left|{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}\right|\geqslant {\mathcal{A}}\quad \text{for }{\mathcal{A}}\gg 1\quad \text{and}\quad \left|{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}}{\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}}}\right|\leqslant {\mathcal{A}}\quad \text{for }{\mathcal{A}}\ll 1.\end{eqnarray}$$

By defining a decoupling threshold such that the interfaces are assumed to be decoupled for $|\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}/\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}|>\unicode[STIX]{x1D712}$ or $|\tilde{\unicode[STIX]{x1D702}}_{1}^{mn}/\tilde{\unicode[STIX]{x1D702}}_{2}^{mn}|<1/\unicode[STIX]{x1D712}$ , we suggest the following criterion to evaluate the importance of wave coupling in LMBs:

(2.41a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D702}_{2}\approx 0\quad \text{ if}\left\{\begin{array}{@{}ll@{}}{\mathcal{A}}>\unicode[STIX]{x1D712}{\displaystyle \frac{h_{3}}{h_{1}}} & \text{for }h_{3}\geqslant h_{1}\\ {\mathcal{A}}>\unicode[STIX]{x1D712} & \text{for }h_{3}<h_{1},\end{array}\right.\end{eqnarray}$$

(2.41b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D702}_{1}\approx 0\text{ if}\left\{\begin{array}{@{}ll@{}}{\mathcal{A}}<{\displaystyle \frac{1}{\unicode[STIX]{x1D712}}}{\displaystyle \frac{h_{3}}{h_{1}}} & \text{for }h_{3}\leqslant h_{1}\\ {\mathcal{A}}<{\displaystyle \frac{1}{\unicode[STIX]{x1D712}}} & \text{for }h_{3}>h_{1}.\end{array}\right.\end{eqnarray}$$

Table 2.

Coupling parameter ${\mathcal{A}}_{gravity}$ calculated for different possible working material combinations. The densities are reported at working temperature $T$ , see Lyon (Reference Lyon1954), Agruss, Karas & Decker (Reference Agruss, Karas and Decker1962), Weaver, Smith & Willmann (Reference Weaver, Smith and Willmann1962), Karas & Mangus (Reference Karas and Mangus1963), Cairns et al. (Reference Cairns, Crouthamel, Fischer, Foster, Hesson, Johnson, Shimotake and Tevebaugh1967), Herédy et al. (Reference Herédy, Iverson, Ulrich and Recht1967), Cairns & Shimotake (Reference Cairns and Shimotake1969a ,Reference Cairns and Shimotake b ), Shimotake, Rogers & Cairns (Reference Shimotake, Rogers and Cairns1969), Swinkels (Reference Swinkels, Braunstein, Mamantov and Smith1971), Janz et al. (Reference Janz, Tomkins, Allen, Downey, Gardner, Krebs and Singer1975, Reference Janz, Allen, Bansal, Murphy and Tomkins1979), Chum & Osteryoung (Reference Chum and Osteryoung1980, Reference Chum and Osteryoung1981), Zinkle (Reference Zinkle1998), Smithells, Gale & Totemeier (Reference Smithells, Gale and Totemeier2004), Sobolev (Reference Sobolev2007), International Atomic Energy Agency (2008), Sobolev (Reference Sobolev2010), Aqra & Ayyad (Reference Aqra and Ayyad2011), Bradwell (Reference Bradwell2011), Bradwell et al. (Reference Bradwell, Kim, Sirk and Sadoway2012), Spatocco, Burke & Sadoway (Reference Spatocco, Burke and Sadoway2014), Wang et al. (Reference Wang, Jiang, Chung, Ouchi, Burke, Boysen, Bradwell, Kim, Muecke and Sadoway2014).

a Aluminium-refining cell (Pearson & Phillips Reference Pearson and Phillips1957)

These conditions guarantee both that only one interface is excited due to (2.39) and the amplitude ratio is becoming large or small enough due to (2.40). $\unicode[STIX]{x1D712}$ has to be chosen arbitrarily depending on the expediency. In the context of the metal pad roll dynamics we numerically found that $\unicode[STIX]{x1D712}=10$ is suitable as it will be discussed in the second part. In the case of the $\text{Mg}\Vert \text{Sb}$ battery numerically investigated with and without interfacial tension by Bojarevics & Tucs (Reference Bojarevics, Tucs and Ratvik2017) and Weber et al. (Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier2017b ) we find high coupling parameters ${\mathcal{A}}\approx 28$ and ${\mathcal{A}}_{gravity}\approx 33$ both clearly exceeding the coupling threshold. Indeed, both studies obtained wave motion essentially in the upper interface propagating with the frequency of the slow mode $\unicode[STIX]{x1D714}_{-}$ . To determine the practical importance of wave coupling in large-size LMBs, ${\mathcal{A}}_{\text{gravity}}$ was calculated for further metals. Table 2 shows the densities at operation temperature and the corresponding coupling parameter ${\mathcal{A}}_{gravity}$ for most common metal and electrolyte combinations concerned as possible working materials for LMBs in the literature. It can be seen that the so far studied $\text{Mg}\Vert \text{Sb}$ cells coincidentally typifies the extreme case of lowest wave coupling. Most batteries reveal values of ${\mathcal{A}}_{gravity}$ between 3 and 6 where motion is present in both interfaces. Hence, coupling dynamics is indeed of high importance for most LMBs and may not be generally neglected as has previously been assumed on the basis of the obtained dynamics in $\text{Mg}\Vert \text{Sb}$ cells.

Interestingly, the density ratios in three-layer aluminium-refining cells are such that the interfaces in those cells are strongly coupled as well (see table 2). The following discussion applies therefore equally to aluminium-refining cells. However, the refining process – not limited by the open-circuit voltage – allows for much thicker electrolyte layers, typically between 8 and 25 cm (Beljajew, Rapoport & Firsanowa Reference Beljajew, Rapoport and Firsanowa1957; Pearson & Phillips Reference Pearson and Phillips1957).

3.1 Numerical set-up and procedure

For the numerical study we used the three-dimensional multiphase DNS solver developed by Weber et al. (Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier2017b ) particularly to investigate MHD instabilities in multilayer systems. The solver is implemented in the open source computational fluid dynamics (CFD) library OpenFOAM. The physical model is sketched in appendix B, further detailed information can be found in Weber et al. (Reference Weber, Beckstein, Galindo, Herreman, Nore, Stefani and Weier2017a ,Reference Weber, Beckstein, Herreman, Horstmann, Nore, Stefani and Weier b ).

Table 3.

Used generic simulation parameters of the standard case yielding ${\mathcal{A}}=1$ .

On the basis of the theoretical results a parameter study around ${\mathcal{A}}=1$ was performed to investigate strongly coupled wave dynamics. To reach ${\mathcal{A}}=1$ , a generic standard case of equal metal-layer heights $h_{i}$ , conductivities $\unicode[STIX]{x1D70E}_{i}$ , density differences $\unicode[STIX]{x0394}\unicode[STIX]{x1D70C}$ , interfacial tensions $\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{1}}$ and $\unicode[STIX]{x1D6FE}_{\unicode[STIX]{x1D702}_{2}}$ as well as viscosities $\unicode[STIX]{x1D708}_{i}$ was defined. The parameters used are shown in table 3. Because of the small height chosen for the electrolyte, the influence of ${\mathcal{B}}$ is very weak such that it is sufficient to describe the coupling only by ${\mathcal{A}}$ . In order to capture also non-shallow water dynamics, all simulations were conducted in a aspect-ratio-one cylindrical tank of equal diameter $D$ and height $H$ of $10~\text{cm}$ . The grid was constructed of purely orthogonal hexaedra with a lateral resolution of 50 cells and an axial resolution of 60 cells. This resolution was found to be sufficient to capture long-wave instabilities by employing a grid refinement study. The results did not significantly change when using finer grids.

Various simulation series were conducted for different values of ${\mathcal{A}}$ . ${\mathcal{A}}$ was modified stepwise by adapting only the density of the upper layer $\unicode[STIX]{x1D70C}_{1}$ . As predicted by the theory, we achieved the same results by modifying other ${\mathcal{A}}$ -dependent parameters. Simulations with different cell currents $I_{0}$ were performed to find the stability onsets by successively increasing $I_{0}$ . Once the current becomes large enough ( $I_{0}>I_{c}$ ) to destabilize the system, small numerical errors will grow and form the metal pad roll after a while. After first noticeable interfacial displacements were identified for $I_{0}>I_{c}$ , we ensured to leave the system sufficient time to reach a steady wave propagation state or to develop a short-circuit. Sufficient means here that we only accepted wave propagation as converged if both, the interface amplitudes and the mean velocity in the cell, were found to be constant for at least 30 periods. By that procedure the propagation regime between the rest state and the short-circuit was identified for all investigated ${\mathcal{A}}$ in dependence on the cell current. Finally, exemplary simulations in the propagation regime were chosen to evaluate the wave frequencies and amplitude ratios. It was ensured that the amplitudes were sufficiently small to allow comparison with linear wave theory.

Figure 6.

Exemplary contour plots of the excited interfaces for coupling parameters ${\mathcal{A}}=0.25$ (a) and ${\mathcal{A}}=2.4$ (b). (c,d ) Cross-sections of the corresponding interface contours $\unicode[STIX]{x1D702}_{1}$ (blue) and $\unicode[STIX]{x1D702}_{2}$ (red) in the $r$ – $z$ planes orientated along the highest amplitude shaping. In addition, the dashed lines mark the rest position of the interfaces and the dotted lines show the radial component of the theoretical gravity–capillary modes due to (2.12a ) and (2.12b ).

3.2 Numerical results

For all investigated coupling parameters interfacial motion could be confirmed in accordance with the theoretical predictions. In general, we found two distinguishable coupling regimes in dependence on ${\mathcal{A}}$ comprising fundamentally different coupling dynamics. The first regime was found for $0.7\lesssim {\mathcal{A}}\lesssim 2$ , where the amplitudes of both interfaces are similarly strong. This regime involves complex dynamics and various possible wave modes that can evolve. We call it from now on the ‘strongly coupled regime’. Respectively, for ${\mathcal{A}}\lesssim 0.7$ and ${\mathcal{A}}\gtrsim 2$ a second regime, later on referred to as the ‘weakly coupled regime’, was found, where always one interface was dominant. In these simulations we indeed obtained only the antisymmetric coupling $\unicode[STIX]{x1D714}_{-}$ as predicted by the theory. At first, characteristic modes observed for both regimes are introduced qualitatively in the following section. For reasons of clarity, the weakly coupled regime is presented first.

3.2.2 Strongly coupled regime

Figure 7.

Exemplary contour plots of the excited interfaces for the coupling parameter ${\mathcal{A}}=0.9$ with the applied cell currents of $I=400~\text{A}$ (a) and $I=500~\text{A}$ (b). These currents lead to antisymmetric and symmetric wave coupling, respectively. (c,d ) Cross-sections of the corresponding interface contours $\unicode[STIX]{x1D702}_{1}$ (blue) and $\unicode[STIX]{x1D702}_{2}$ (red) in the $r$ – $z$ planes orientated along the highest amplitude shaping. (c) Shows in addition the fitted curves due to (4.3a ) and (4.3b ). (e,f ) Depict further the compensation current vectors in the same cross-section. A superimposed colour map highlights the radial (horizontal) component of the compensation currents.

For coupling parameters closer to one ( $0.7\lesssim {\mathcal{A}}\lesssim 2$ ) far more complex dynamics was observed. Note first that these regime thresholds are purely empirical at this stage. They may change in dependence on $h_{2}$ (affecting ${\mathcal{B}}$ ) and possibly on further parameters not occurring in ${\mathcal{A}}$ and ${\mathcal{B}}$ like the viscosities and electrical conductivities. In this regime, we found different coupling modes and excited wavenumbers as well as a non-oscillatory instability in dependence on the cell current. Further, time-dependent mode transitions were found such that one mode can merge into another one spontaneously making analysis more complicated. However, two characteristic flow states were found that are mainly present in this regime exemplarily shown in figure 7. For relatively low applied cell currents slightly above the observed stability thresholds, we always found oppositely directed radial-symmetric rotating interfaces, as exemplarily shown in figure 7(a) and figure 7(c). This state largely differs from the MPRI in the physical sense since we observed that these vortex shapes can persist for minutes or even for the complete simulation depending on the cell currents. Hence, it can be considered as a quasi-steady non-oscillatory instability that cannot be described by the wall reflection of oscillating interfacial waves, compare to § 1. Figure 7(e) shows in addition the induced radial compensation currents in cross-section. In contrast to the metal pad roll mode, two closed current loops can be identified leading to pure azimuthal Lorentz forces inducing oppositely directed circular motions in the liquid electrodes. In § 4.1 we give a detailed explanation how centripetal pressure differences can cause this instability.

For higher applied cell currents, a second characteristic state was found to develop: symmetrically coupled in-phase rotating waves with the expected wavenumber $\unicode[STIX]{x1D716}_{11}$ associated with the metal pad roll, as depicted in figure 7(b,d). Also for high cell currents close to the short-circuit onset, the radial instability (figure 7 a) arose first. But for sufficiently high cell currents this coupling state always spontaneously turned into the symmetric mode (figure 7 b) after some time, then remained preserved, or – in only very few cases – once again turned into the previously presented antisymmetric mode of the weakly coupled regime (§ 3.2.1). This instability even strongly differs from the common MPRI mechanism. It seems to be counterintuitive at first glance since the metal pad is always explained by a redistribution of the cell currents due to changes of the electrical resistivity caused by electrolyte-layer deformations. But for parallel interfaces there is no preferred direction for the currents to choose because the salt-layer height remains constant everywhere. Thus, horizontal compensation currents are not expected to arise. However, it can be seen from figure 7(f) that here horizontal currents indeed appear in the electrolyte. This is in contrast to the common assumption stating that the current remains almost vertical in the low conductivity electrolyte. A detailed explanation unravelling that issue is given in § 4.2.

For the sake of completeness, it should be noted that besides the modes presented here we sometimes observed also higher modes or mode superpositions in particular for high cell currents near the short-circuit onset. However, these waves were only temporarily stable and seemed to appear randomly. All in all, the two modes presented here were highly dominant and can be considered as the characteristic modes of the strongly coupled regime. However, the mode dynamics was found to be far more complex and cannot be precisely predicted at this stage.

Interestingly, the observed transition from antisymmetrically coupled waves to symmetrically coupled waves was also numerically found recently by Zikanov (Reference Zikanov2017) in rectangular cells. He also adjusted the density differences (corresponding to the derived parameter ${\mathcal{A}}_{gravity}$ ) and found the symmetric mode only for strongly coupled interfaces as well, where ${\mathcal{A}}_{gravity}$ is close to one. This agreement is remarkable since the used numerical approaches and cell geometries are quite different in both studies. Hence, we expect that the mode transition is indeed likely to occur in future large LMBs and may be of high practical importance.

3.3 Comparison with the potential flow theory

Both the amplitude ratios and the wave frequencies are quantitatively analysed and compared with the formulas (2.21) and (2.22). To calculate the numerical amplitude ratios, we tracked the minimum and maximum displacement of both interfaces at every time step. Since the interfaces are not necessarily evenly displaced around the initial positions, as shown in figure 6(c,d), the amplitude ratios are calculated by the differences of the maximum and minimum displacements

(3.1) $$\begin{eqnarray}{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}_{1}}{\tilde{\unicode[STIX]{x1D702}}_{2}}}={\displaystyle \frac{\max (\unicode[STIX]{x1D702}_{1})-\min (\unicode[STIX]{x1D702}_{1})}{\max (\unicode[STIX]{x1D702}_{2})-\min (\unicode[STIX]{x1D702}_{2})}}.\end{eqnarray}$$

In order to reduce uncertainties, $\tilde{\unicode[STIX]{x1D702}}_{1}/\tilde{\unicode[STIX]{x1D702}}_{2}$ was evaluated for a wide range of time steps comprising at least 50 periods and averaged. The wave periods were calculated by tracking the angle of highest position of the interfaces over time and fitting the angular velocity. Figure 8(a,b) shows the amplitude ratios determined in such a way and wave periods for all evaluated ${\mathcal{A}}/{\mathcal{B}}$ in comparison with the potential curves. For clarity, we only show the dominant symmetric mode $\unicode[STIX]{x1D714}_{+}$ in the strongly coupled regime marked by two vertical lines. Both the numerical amplitude ratios as well as the wave periods are in agreement with the theory for both modes. Consequently, though the interfacial shapes are considerably deviating from the potential solutions, it is confirmed that – at least to very good approximation – eigenmodes of pure gravity–capillary waves are excited by the induced Lorentz forces. Further it can be concluded that the wave coupling indeed is mainly determined by pressure exchanges. Contributions of induced current redistributions to the coupling as well as viscous effects seem to be of little relevance. Hence, the highly simplified (and often underestimated) irrotational potential flow theory is sufficient to approximate the coupling dynamics in three-layer LMBs driven by rotational Lorentz forces.

4.1 Bulge instability

The BI, that was found only in the strongly coupled regime, can be regarded as a special instability case since it appeared as a quasi-steady state that may persist over minutes. Hence, it can not be described by time-dependent gravity–capillary waves. Rather, it should be considered as a further fixed point existing in the system. In fact, stationary interface perturbations also have been observed before in physically similar ARCs by Bojarevics & Pericleous (Reference Bojarevics and Pericleous2006) and were later quantitatively described by Munger & Vincent (Reference Munger and Vincent2008) for cylindrical geometries. They identified differences between the centripetal pressure of both phases that have to be balanced by interface deformations, as the key mechanism. This phenomenon is widely present in everyday life and can be seen, e.g. when stirring a cup of tea. Figure 9 illustrates the driving mechanism due to Munger & Vincent (Reference Munger and Vincent2008) transferred to three-layer LMBs. At first, the explanation is the same as for the MPRI. Radial-symmetrically displaced interfaces will lead to also radial-symmetric compensation currents $\boldsymbol{J}_{c}$ since the perturbed current $\boldsymbol{J}_{p}$ chooses the path near the circular tank walls where the electrolyte is thin. In interaction with $b_{z}$ , these currents lead to pure azimuthal Lorentz forces $\boldsymbol{F}_{L}$ driving two large-scale vortices spinning clockwise in the anode layer and counterclockwise in the cathode layer. The vortices can be described by a priori unknown average radial angular velocities profiles $\unicode[STIX]{x1D6FA}_{1}(r)$ and $\unicode[STIX]{x1D6FA}_{3}(r)$ . No average rotation is expected to arise in the electrolyte layer $\unicode[STIX]{x1D6FA}_{2}\approx 0$ . Thereby, centripetal pressure jumps arise at the interfaces reading

(4.1a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x0394}P|_{\unicode[STIX]{x1D702}_{1}}={\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}_{1}\unicode[STIX]{x1D6FA}_{1}^{2}r^{2}>0,\end{eqnarray}$$

(4.1b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x0394}P|_{\unicode[STIX]{x1D702}_{2}}=-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D70C}_{3}\unicode[STIX]{x1D6FA}_{3}^{2}r^{2}<0.\end{eqnarray}$$

Figure 9.

Schematic vertical cut of the BI. The current $\boldsymbol{J}_{p}$ is redirected to the tank wall giving rise to horizontal compensation currents $\boldsymbol{J}_{c}$ which, in interaction with an external vertical magnetic field $b_{z}$ , in turn lead to azimuthal Lorentz forces $\boldsymbol{F}_{L}$ inducing counter-rotating vortices with angular velocities $\unicode[STIX]{x1D6FA}_{1}$ and $\unicode[STIX]{x1D6FA}_{3}$ . The concave and convex interface shapes then arise to equilibrate centripetal pressure differences.

These pressure drops can be equilibrated only by a concave deformation of $\unicode[STIX]{x1D702}_{1}$ locally decreasing the hydrostatic pressure $P_{g}$ and a convex deformation of $\unicode[STIX]{x1D702}_{2}$ increasing $P_{g}$ , respectively. To provide a more quantitative comprehension that the BI is caused by centripetal pressure drops, we further studied the interfacial shapes. Since it is very complicated to calculate Lorentz force-induced angular velocities profiles, we assumed that the Lorentz force is acting like a mechanical stirrer. For unbaffled stirred tanks, the angular velocity can be described by (Platzer & Noll Reference Platzer and Noll1983)

(4.2) $$\begin{eqnarray}\unicode[STIX]{x1D6FA}(r)=\unicode[STIX]{x1D6FA}_{0}\left[{\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{r}{r_{0}}}\right)^{4}\right]^{-1/2},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FA}_{0}$ and $r_{0}$ denote the rotation frequency and characteristic length of the impeller. Applied to the metal layers, this profile leads to the theoretical vortex shapes

(4.3a ) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D702}_{1}(r)=h_{1}-\tilde{\unicode[STIX]{x1D702}}_{1}-{\displaystyle \frac{(\unicode[STIX]{x1D6FA}_{0,1}r_{0,1})^{2}}{g}}\arctan \left({\displaystyle \frac{r}{r_{0,1}}}\right)^{2},\end{eqnarray}$$

(4.3b ) $$\begin{eqnarray}\displaystyle & & \displaystyle \unicode[STIX]{x1D702}_{2}(r)=h_{3}-\tilde{\unicode[STIX]{x1D702}}_{2}+{\displaystyle \frac{(\unicode[STIX]{x1D6FA}_{0,2}r_{0,2})^{2}}{g}}\arctan \left({\displaystyle \frac{r}{r_{0,2}}}\right)^{2}.\end{eqnarray}$$

$\unicode[STIX]{x1D6FA}_{0,1}$

$\unicode[STIX]{x1D6FA}_{0,2}$

$r_{0,1}$

$r_{0,2}$

Figure 10.

Schematic explanation of the STI. (a) Profile of an LMB showing schematically the paths of the unperturbed cell current $\boldsymbol{J}_{0}$ , the perturbation current $\boldsymbol{J}_{p}$ caused by the symmetrically tilted interfaces and the resulting horizontal compensation current $\boldsymbol{J}_{c}$ . $\boldsymbol{J}_{c}$ closes here in the upper layer of higher electrical conductivity $\unicode[STIX]{x1D70E}_{1}>\unicode[STIX]{x1D70E}_{3}$ . The interaction of $\boldsymbol{J}_{c}$ with an external vertical magnetic field $b_{z}$ leads to azimuthal Lorentz forces mainly in the electrolyte driving the rotating waves. (b) Electrolyte-layer section. The perturbation current $\boldsymbol{J}_{p}$ in the electrolyte can be decomposed into a component $\boldsymbol{J}_{p\Vert }$ parallel to $b_{z}$ and a component $\boldsymbol{J}_{p\bot }$ vertical to $b_{z}$ . The ratio of both components is connected by a pitch angle $\unicode[STIX]{x1D6FC}$ to the slope of the interfaces.

For ARCs Munger & Vincent (Reference Munger and Vincent2008) found that non-oscillatory instabilities evolve half as fast as the MPRI, which is why they were rarely found in the simulations. However, in three-layer LMBs strongly coupled interfaces favour this instability, because the oppositely rotating liquid metals produce stronger shearing forces in the electrolyte. That causes higher pressure drops at the interfaces in comparison to ARCs, where the aluminium and the cryolite rotate in the same direction with different angular velocities. Hence, the centripetal pressure differences are even stronger such that the growth rate of the BI can apparently exceed the growth of the MPRI, when both interfaces are simultaneously excited. Either way, since the chosen cylindrical geometry naturally favours the forming of radial-symmetric vortex shapes, we do not expect the BI being that significant in rectangular containers.

4.2 Synchronous tilting instability

As already mentioned in § 3.2.2, the discovered STI can not be explained by the Sele mechanism alone, neither (see § 1). This follows from the fact that the electrolyte-layer height remains (almost) preserved such that no preferred direction for the current is expected to develop. We observed that the STI is considerably ‘stabilizing’ in the manner that higher cell currents are necessary to provoke short-circuits in comparison to the MPRI of the weakly coupled regime. Hence, it is highly desirable to understand the physical mechanism behind the STI that can deliver practical benefits for LMB operation. Figure 10(a) clarifies schematically the driving mechanism we assume behind the STI. The key assumption is, again, that the current always takes the shortest path through the electrolyte due to its orders higher electrical resistivity. When the interfaces are displaced in parallel, the perturbed current $\boldsymbol{J}_{p}$ must pass perpendicularly through the electrolyte, as sketched in figure 10(a). This current has, in dependence on the electrolyte layer’s inclination, a horizontal component leading to Lorentz forces acting within the electrolyte and driving the rotation of both interfaces simultaneously. Further, the horizontal component can be regarded as acting like the compensation current $\boldsymbol{J}_{c}$ of the classical mechanism described in figure 1. Similarly, it has to close either in the anode or cathode layer (as long the electrical conductivities of the liquid metals are larger than those of the current collectors (Munger & Vincent Reference Munger and Vincent2006)). It will close in the better conducting metal as shown in figure 10(a) for the case $\unicode[STIX]{x1D70E}_{1}>\unicode[STIX]{x1D70E}_{3}$ . Here, the current in the anode layer is deflected to the left in order to minimize the current path through the more resistive cathode layer. Additional simulations with varying electrical conductivities, not shown here, confirmed this explanation. For equal electrode-layer conductivities, as used in the coupling study, the stronger displaced interface and thus ${\mathcal{A}}$ determined the closing layer.

Compared to the MPRI, the STI allows much higher currents without leading to a short-circuit. We will give here a possible explanation by comparing the horizontal perturbation currents of both instabilities. The horizontal currents arising in the STI can be quantitatively estimated by decomposing the perpendicular perturbation current $\boldsymbol{J}_{p}$ into a component $\boldsymbol{J}_{p\Vert }$ parallel to the magnetic field $b_{z}$ , which does not contribute to the Lorentz force, and a component $\boldsymbol{J}_{p\bot }$ vertical to $b_{z}$ . $\boldsymbol{J}_{p\bot }$ is acting like the horizontal compensation current of the Sele mechanism $\boldsymbol{J}_{p\bot }\hat{=}\boldsymbol{J}_{c}$ and drives the instability. Figure 10(b) shows these decompositions applied to a small electrolyte-layer section appearing approximately linear. By assuming that $\boldsymbol{J}_{p}$ is equally distributed in the electrolyte layer, we identify $|\boldsymbol{J}_{p}|\approx |\boldsymbol{J}_{0}|$ within that layer since the three layers can be considered as a serial connection. From figure 10(b) we find the relation

(4.4) $$\begin{eqnarray}\sin (\unicode[STIX]{x1D6FC})={\displaystyle \frac{J_{p\bot }}{J_{0}}}\end{eqnarray}$$

that can be further equalized with the interfacial slope. For the purpose of a conservative estimate the interfaces are approximated only by the radial component $\unicode[STIX]{x1D702}_{r}$ in direction of the strongest interfacial displacement. From (2.14a ) the radial slope of both interfaces is given by

(4.5) $$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D702}_{r}(r)}{\unicode[STIX]{x2202}r}}=\tilde{\unicode[STIX]{x1D702}}{\displaystyle \frac{\unicode[STIX]{x1D716}_{11}}{R}}{\displaystyle \frac{1}{2}}\left(J_{0}\left(\unicode[STIX]{x1D716}_{11}{\displaystyle \frac{r}{R}}\right)-J_{2}\left(\unicode[STIX]{x1D716}_{11}{\displaystyle \frac{r}{R}}\right)\right)\end{eqnarray}$$

leading to the $r$ -dependent horizontal perturbation current

(4.6) $$\begin{eqnarray}|\boldsymbol{J}_{p\bot }|=|\boldsymbol{J}_{0}|\tilde{\unicode[STIX]{x1D702}}{\displaystyle \frac{\unicode[STIX]{x1D716}_{11}}{R}}{\displaystyle \frac{1}{2}}\left(J_{0}\left(\unicode[STIX]{x1D716}_{11}{\displaystyle \frac{r}{R}}\right)-J_{2}\left(\unicode[STIX]{x1D716}_{11}{\displaystyle \frac{r}{R}}\right)\right).\end{eqnarray}$$

In order to estimate the global impact to the interfaces, both the maximum horizontal current $\max |\boldsymbol{J}_{p\bot }|$ at $r=0$ as well as the mean current $\langle |\boldsymbol{J}_{p\bot }|\rangle _{r}$ are evaluated

(4.7) $$\begin{eqnarray}\displaystyle & & \displaystyle \max |\boldsymbol{J}_{p\bot }|=|\boldsymbol{J}_{0}|{\displaystyle \frac{\unicode[STIX]{x1D716}_{11}}{2}}{\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}}{R}},\end{eqnarray}$$

(4.8) $$\begin{eqnarray}\displaystyle & & \displaystyle \langle |\boldsymbol{J}_{p\bot }|\rangle _{r}=|\boldsymbol{J}_{0}|J_{1}(\unicode[STIX]{x1D716}_{11}){\displaystyle \frac{\tilde{\unicode[STIX]{x1D702}}}{R}}.\end{eqnarray}$$

The maximum current $\max |\boldsymbol{J}_{p\bot }|$ can be easily exploited for comparisons with the maximum compensation currents obtained in the simulation. For the case shown in figure 7(f) we find for example an average amplitude of $\tilde{\unicode[STIX]{x1D702}}\approx 0.85~\text{cm}$ with the applied current density of $|\boldsymbol{J}_{0}|=6.37~\text{A}~\text{cm}^{-2}$ and $R=5~\text{cm}$ . These values lead to maximum horizontal currents of $\max |\boldsymbol{J}_{p\bot }|\approx 1~\text{A}~\text{cm}^{-2}$ matching well with the highest compensation currents occurring in figure 7(f). Good agreement was found as well for the other simulations within the strongly coupled regime, giving a promising indication that indeed simple vector decomposition is sufficient to approximate horizontal compensation currents caused by the symmetric wave mode.

Finally, we analyse the potential of the STI to delay the short-circuit in comparison to the MPRI present in the weakly coupled or decoupled regimes. Since it is unclear how the wave coupling of the antisymmetric mode in the weakly coupled regime exactly effects the stability, the classical Sele mechanism describing the MPRI in ARCs is used for comparisons. As known from the literature (Gerbeau et al. Reference Gerbeau, Le Bris and Lelièvre2006), in these systems the induced horizontal compensation currents $|\boldsymbol{J}_{c}|_{ARC}$ in the aluminium layer are of the order

(4.9) $$\begin{eqnarray}|\boldsymbol{J}_{c}|_{ARC}\sim {\displaystyle \frac{|\boldsymbol{J}_{0}|\tilde{\unicode[STIX]{x1D702}}_{1}2R}{h_{1}h_{2}}}.\end{eqnarray}$$

Transferred to the LMBs, this expression is a valid estimation for ${\mathcal{A}}/{\mathcal{B}}\gg 1$ (the estimation can be done just as well for ${\mathcal{A}}/{\mathcal{B}}\ll 1$ ), where only the upper interface is excited by the slow mode $|\boldsymbol{J}_{c}|_{ARC}\,\hat{=}\,|\boldsymbol{J}_{c}|_{LMB,\unicode[STIX]{x1D714}_{-}}$ . The ratio of $|\boldsymbol{J}_{c}|_{LMB,\unicode[STIX]{x1D714}_{-}}$ to the mean horizontal current produced by the STI (4.8) that is induced by the same interface displacements, scales with

(4.10) $$\begin{eqnarray}{\displaystyle \frac{|\boldsymbol{J}_{c}|_{LMB,\unicode[STIX]{x1D714}_{-}}}{|\boldsymbol{J}_{c}|_{LMB,\unicode[STIX]{x1D714}_{+}}}}\sim {\displaystyle \frac{R^{2}}{h_{1}h_{2}}}.\end{eqnarray}$$

Since the horizontal currents drive both instabilities, we can conclude that the MPRI will lead faster to short-circuits than the STI because we find $h_{1},h_{2}<R$ in practice. This effect will be even more prominent in shallow batteries, as they are desired for practical operation in order to minimize ohmic losses. Hence, this stabilizing mechanism can be exploited to potentially improve the operation safety of LMBs e.g. by modifying the density of the electrolytes to approach the strongly coupled regime. However, we provide only a first understanding of this novel instability mechanism. Quantitative three-layer stability analysis and further numerical studies, which are necessary to gain a deeper understanding, are left for future work.