Magneto-hydrodynamic stability of a liquid metal battery in discharge

With the advancement of renewable intermittent energy sources the idea of Liquid Metal Battery (LMB) as a device for stationary energy storage is gaining attention [1]. The development of small-scale (about few cm) prototypes demonstrates technical feasibility of the battery and its advantages to existing solid batteries [1,2]. The liquid state of the main components necessitates consideration of the fluid dynamics in the LMBs. A number of recent publications are devoted to the problem: a possibility of the Tayler instability [3,4], the thermal convection [5], observation of a vortical flow in the LMB model [6], as well as the electromagnetically driven wave instability in the stratified three-liquid-layer system [7–12]. The main motivation of these investigations is to gain understanding of the liquid mixing, the charge and mass tansfer, and a possibility of electro-magnetically driven destabilizing interface motion.

LMBs are thought to be easily scalable due to their simple construction using the natural density stratification of the liquid layers. The magneto-hydrodynamic stability of the liquid filled cells could be problematic with the upscaling. As was shown in [7], large cells of several cubic meters total volume can be designed to operate safely at very high power value (${\sim}10^{5}$ W). The same high power and storage capacity could be achieved using a large cluster of small LMBs [1,8]. However in both cases high current densities in the cells will be coupled to the total magnetic field (created by the currents in the cell, the supply bars and the neighbour cells) leading to significant electromagnetic forces. Such forces in stratified liquid layers with free surface areas may cause the long-wave interfacial instability as is well known in the case of Hall-Héroult Cells (HHC) [13–15]. The manifestation of this instability in the small-scale LMB at various stages of discharge, when the metal levels change, is the subject of this paper. The previously developed theoretical and numerical models [12] will be applied to analyze the interfacial waves in the small-scale LMB filled with three liquid layers (fig. 1) composed of Na∣ NaF-NaCl-NaI∣ Bi [1,2]. The liquid sodium is considered as one of attractive and relatively safe anodic materials to be used in the upscaling physical experiments.

Zoom In Zoom Out Reset image size

The numerical models will be applied to answer the questions about the two interface coupling, the importance of the initial perturbation conditions and the destabilizing effect of the metal level change during the charging/discharging process.

The model of battery consisting of three density stratified electrically conductive liquid layers contained in a rectangular box is schematically represented in fig. 1. The full mathematical model derivation is presented in [12], therefore we will emphasize here only the main points leading to the coupled wave and electric current equations.

Coupled interfacial dynamics

The incompressible laminar fluid flow in the presence of electro-magnetic fields is described by the equations

$$\begin{eqnarray} \rho \partial_{t}u_{i}+\rho u_{j}\partial_{j}u_{i}&= -\partial_{i}(p+\rho gz)+\rho \nu \partial_{jj}u_{i}+f_{i}, \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} \partial_{i}u_{i}&= 0, \end{eqnarray} \tag{ 2 }$$

where the indices $i,j=1,2,3$ correspond to the coordinates (x, y, z), the velocity components are given as (u₁, u₂, u₃), the summation over repeated indices is implied, ρ represents the density, ν the kinematic viscosity, g the gravitational acceleration, p the pressure and $\bm{f}=\bm{j}\times \bm{B}$ stands for the electromagnetic force vector, whereas j is the electric current density and B the magnetic field. The horizontal dimensions of the cell are assumed to be much larger compared to the vertical depth of layers, so that the description can be based on a systematically derived shallow-water approximation. The velocity components can be represented as an expansion in a small aspect ratio parameter $\delta =\max (h)/\min(L)$, where h is a typical depth (e.g., of the unperturbed layer) and L is the characteristic horizontal dimension (width of the cell):

$$\begin{equation} u_{i}=u_{0i}( x,y,t)+{\delta u}_{1i}( x,y,z,t)+O ( \delta^{2}). \end{equation} \tag{ 3 }$$

The leading horizontal components of (1) for $i=1,2$ are

$$\begin{equation} \rho \partial_{t}u_{0i}+\rho u_{0j}\partial_{j}u_{0i}=-\partial_{i}p+{\rho k_{f}U}_{i}+f_{i}, \end{equation} \tag{ 4 }$$

for the depth averaged velocity components when using the shallow-water linear friction law (the value of k_f is given below in (26)). The depth averaged horizontal velocity components are defined as

$$\begin{equation} U_{ik}=h_{k}^{-1}\int_{H_{k-1}}^{H_{k}}u_{ik}(x,y,z)\mathrm{d}z, \end{equation} \tag{ 5 }$$

where $k=1,2,3$ is the layer number, $h_{k}( x,y,t)=H_{k}-H_{k-1}$ is the local variable depth, see the fig. 1.

The vertical component i = 3 of the eq. 1 gives the hydrostatic pressure distribution in the liquid layers accounting for the variable interface shapes H₁(x, y, t) and H₂(x, y, t) [12]. The depth averaged continuity eq. (2) is accurate to all orders of approximation:

$$\begin{equation} \partial_{t}h_{k}+\partial_{i}( U_{ik}h_{k} )=0. \end{equation} \tag{ 6 }$$

The equations can be linearized if an additional approximation of a small wave amplitude is introduced: $h_{k}( x,y,t )=h_{0k}+\varepsilon h_{k}^{\prime}( x,y,t )$ for the layer thickness, or equivalently $H_{k}( x,y,t )=H_{0k}+\varepsilon H_{k}^{\prime}( x,y,t )$ for the interface position, $\varepsilon =\max {(A)/h}$, where A is a typical wave amplitude and h_0k, H_0k are the unperturbed values. For typical geometries considered in this paper $\delta =h /L=(0.01\div 0.05)\ \text{m}/0.2\ \text{m}=0.05\div 0.25<1$, and $\varepsilon =A/h= 0.001\ \text{m}/0.01\ \text{m}=0.1\ll 1$. The small depth approximation $\delta \ll 1$ is just marginally satisfied in some cases considered below, however the 3-layer wave frequencies are mostly affected by the presence of the shallow intermediate layer ($h_{2}=0.01\ \text{m}$), therefore permitting a meaningful approximation in the 3-layer problem.

The wave equations for the coupled interfaces H₁ and H₂ can be derived following the procedure described in [12]:

$$\begin{eqnarray} & \alpha_{1}\partial_{tt}H_{1}+k_{fe1}\partial_{t}H_{1}-\rho_{2}/h_{02}\partial_{tt}H_{2}\nonumber\\ & {}-{\rho_{2}k_{f2}}/ h_{02}\partial_{t}H_{1}=R_{1}\partial_{ii}H_{1}+\partial_{i}( F_{i2}-F_{i1})\nonumber\\ & {}+[ {\rho_{1}\partial }_{i}( U_{j1}\partial_{j}U_{i1} )-{\rho_{2}\partial }_{i}( U_{j2}\partial_{j}U_{i2} )], \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} & \alpha_{2}\partial_{tt}H_{2}+k_{fe2}\partial_{t}H_{2}-\rho_{2}/ h_{02}\partial_{tt}H_{1}\nonumber\\ & {}-{\rho_{2}k_{f2}}/ h_{02}\partial_{t}H_{1}=R_{2}\partial_{ii}H_{2}+\partial_{i}( F_{i3}-F_{i2} )\nonumber\\ & {}+ [ {\rho_{2}\partial }_{i}( U_{j2}\partial_{j}U_{i2} )-{\rho_{3}\partial }_{i}( U_{j3}\partial_{j}U_{i3} )], \end{eqnarray} \tag{ 8 }$$

where the depth averaged force density F_i is defined similarly to (5). The electromagnetic terms are of the same order of magnitude as the leading terms, while the nonlinear wave motion terms are of lower order (${\sim}\varepsilon$). The boundary condition for the normal velocity $u_{n}=0$ at the side walls leads to the condition [12]:

$$\begin{equation} \partial_{n}H_{m}=( F_{n,m}-F_{n,m+1} )/ R_{m}, \end{equation} \tag{ 9 }$$

for $m=1,2$. The constants introduced in the above equations are defined as

$$\begin{eqnarray} \alpha_{m}&= \rho_{m}/h_{0m}+\rho_{m+1}/ h_{0m+1},\nonumber\\ R_{m}&= ( \rho_{m}-\rho_{m+1} )g,\\ k_{fe,m}&= {\rho_{m}k_{fm}}/h_{0m}+ {\rho_{m+1}k_{fm+1}}/{h_{0m+1}}.\nonumber \end{eqnarray} \tag{ 10 }$$

Note that in (7) and (8) the summation over repeated indices is limited to the two horizontal dimensions $i,j=1,2$. Equations (7), (8) include additional terms to the similar wave equations obtained in [11], particularly the viscous friction, generalised electromagnetic force and the horizontal velocity effects. The MHD wave model for the aluminium electrolysis cell [15] with the single interface H₁(x, y, t) can be recovered from (7) if the top interface is constant in time: $H_{2}=H_{2}(x,y)$.

If neglecting the electro-magnetic and the dissipation terms, eqs. (7)–(9) can be solved to obtain the expressions for the 2 coupled interface gravity wave frequencies [9,11,12]:

$$\begin{align} &\hspace{-7pt}\omega_{12,\bm{k}}^{2}=\nonumber\\ &\hspace{-7pt}\frac{-(\omega_{1,\bm{k}}^{2}+\omega_{2,\bm{k}}^{2})\pm[(\omega_{1,\bm{k}}^{2}\!-\!\omega_{2,\bm{k}}^{2})^{2}+4R_{c,1}R_{c,2}\omega_{1,\bm{k}}^{2}\omega_{2,\bm{k}}^{2}]^{1/2}}{2( 1-R_{c,1}R_{c,2} )}, \end{align} \tag{ 11 }$$

$$\begin{equation*} R_{c,m}=\rho_{2}/ ( {\alpha_{m}h}_{02} ),\qquad\omega_{m,\bm{k}}^{2}={R_{m}\bm{k}^{2}} / \alpha_{m}, \end{equation*}$$

showing the correlation between the uncoupled shallow-layer gravity wave frequencies $\omega_{m,\bm{k}}$ and the coupled two interface frequencies $\omega_{12,\bm{k}}$. The wave vector k is defined as

$$\begin{equation} \bm{k}=(k_{x},k_{y}),\quad k_{x}={m\pi } / L_{x},\quad k_{y}={n\pi } / L_{y};\quad m, n\in N. \end{equation} \tag{ 12 }$$

Equation (11) describes the two sets of self-sustained wave frequencies on both interfaces with the "−" sign corresponding to the so-called "slow" mode and the "+" sign to the "fast" mode [9,12]. This expression will be useful to compare the pure gravity wave frequencies to the electromagnetically modified waves and the frequency at which instability starts

Electric current flow

The LMB operates in a dynamic regime of charging or discharging, resulting in the current flowing upwards or downwards. The current flow at the discharging stage in the layered structure is illustrated in fig. 1. We assume that the characteristic time-scale for the wave motion is much larger than the diffusion time of the magnetic field to satisfy the low magnetic Reynolds number approximation and the flow effect on the instantaneous current distribution can be neglected, so that the electric current can be described by a set of coupled Laplace equations for the electric potential $\varphi_{k}( x,y,z,t )$:

$$\begin{equation} \partial_{ii}\varphi_{k}=0, \end{equation} \tag{ 13 }$$

where $k=1,2,3$ corresponds to the layer number. The continuity conditions for the electric potential and the normal current component $\bm{j}\cdot \bm{n}$ at the interfaces $z=H_{m}\ (m=1,2)$ are

$$\begin{eqnarray} &&\varphi_{m}=\varphi_{m+1},\qquad \sigma_{m+1}\partial_{n}\varphi_{m+1}=\sigma_{m}\partial_{n}\varphi_{m}. \end{eqnarray} \tag{ 14 }$$

The normal derivatives at the deformed interfaces are defined as

$$\begin{equation} \partial_{n}\varphi_{k}=\frac{\partial_{z}\varphi_{k}-\partial_{i}H_{k}\partial_{i}\varphi_{k}}{( 1+\partial_{i}H_{k}\partial_{i}H_{k} )^{1/2}}, \end{equation} \tag{ 15 }$$

assuming the summation over the repeated index i only. The side walls of the domain are considered to be electrically insulating. The applied current distributions at the top and the bottom of the cell could be given from the external circuit solution [13], however here for simplicity we assume the supply current to be uniform and equal:

$$\begin{equation} ( j )_{z=H_{0}}=( j )_{z=H_{3}}=j_{0}. \end{equation} \tag{ 16 }$$

The shallow-layer approximation for the 3d electric current requires the expansion of the electric potential in terms of the parameters δ and ε [12,15]:

$$\begin{eqnarray} \varphi_{k}( x,y,z,t )&= \varphi_{0k}+\delta^{2}\varphi_{1k}+o( \delta^{2} ). \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} \varphi_{0k}&= ( a_{k}+\varepsilon A_{k} )z+( b_{k}+\varepsilon B_{k} ), \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} \varphi_{1k}&= ( c_{k}+\varepsilon C_{k} )z+z^{2}\partial_{ii}( b_{k}+\varepsilon B_{k} ), \end{eqnarray} \tag{ 19 }$$

where the 2d functions a_k, A_k, b_k, B_k, c_k, C_k(x, y) are determined from the set of interfacial conditions (14). After a lengthy derivation the final equations for the perturbation potentials defined as ${\Phi }_{1}=\varepsilon B_{1}$ and ${\Phi }_{3}=\varepsilon B_{3}$ can be obtained [12]:

$$\begin{eqnarray} h_{01}h_{02}\partial_{ii}{\Phi }_{1}+\frac{\sigma_{2}}{\sigma_{1}}( {\Phi }_{3}-{\Phi }_{1} )&= \frac{j_{0}}{\sigma_{1}}( h_{2}-h_{02} ), \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} h_{02}h_{03}\partial_{ii}{\Phi }_{3}+\frac{\sigma_{2}}{\sigma_{3}}( {\Phi }_{3}-{\Phi }_{1} )&= -\frac{j_{0}}{\sigma_{3}}( h_{2}-h_{02} ). \end{eqnarray} \tag{ 21 }$$

According to (20) and (21), the perturbed current flow is determined by the electrolyte thickness variations and the ratio of the electrical conductivities. The dimensional current components in the three layers can be expressed as

$$\begin{eqnarray} \bm{j}_{1}&= -\sigma_{1}\biggl( \partial_{x}{\Phi }_{1},\partial_{y}{\Phi }_{1},\frac{j_{0}}{\sigma_{1}}+( H_{0}-z )\partial_{ii}{\Phi }_{1} \biggr), \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} \bm{j}_{2}&= -\sigma_{2}\biggl( 0,0,\biggl( 1-\frac{h_{2}-h_{02}}{h_{02}} \biggr)\frac{j_{0}}{\sigma_{2}}+\frac{{\Phi }_{3}-{\Phi }_{1}}{h_{02}} \biggr), \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} \bm{j}_{3}&= -\sigma_{3}\biggl( \partial_{x}{\Phi }_{3},\partial_{y}{\Phi }_{3},\frac{j_{0}}{\sigma_{3}}+( H_{3}-z )\partial_{ii}{\Phi }_{3} \biggr). \end{eqnarray} \tag{ 24 }$$

To the given perturbation accuracy the horizontal currents in liquid metal layers are depth independent and the vertical current is linearly varying from the imposed supply value to the wave perturbed distribution at the electrolyte. The current distribution in the electrolyte is almost purely vertical due to the relatively low conductivity $\sigma_{2}\ll \sigma_{1}\sim \sigma_{3}$. The terms in the left side of (20), (21) are of equal order proportional to $h^{2}\sim \sigma_{2}/\sigma_{1}\sim \sigma_{2}/\sigma_{3}$. Similar equations to (20), (21), however without the terms containing the perturbed potential differences, have been derived in [11].

Numerical solution construction

The solution of the problem can be constructed using the spectral representation in Fourier space [12,15]. The following set of functions is introduced:

$$\begin{align} &\{ H_{1},H_{2},{\Phi }_{1},{\Phi }_{3} \}( x,y,t )=\nonumber\\ &{} \sum_{\bm{k}} {\{ \hat{H}_{1,\bm{k}},\hat{H}_{2,\bm{k}},\hat{\Phi }_{1,\bm{k}},\hat{\Phi }_{3,\bm{k}}\}( t ) \epsilon_{\bm{k}}\cos ( k_{x}x )\cos {( k_{y}y ).}} \end{align} \tag{ 25 }$$

$\hat{H}_{\bm{k}}( t )$ and $\hat{\Phi }_{\bm{k}}( t )$ are the spectral wave amplitudes and the perturbed potentials in Fourier space, k is defined in (12) and $\epsilon_{\bm{k}}$ is the normalization coefficient. Equations (7), (8), (20), (21) are represented in the spectral space using the expansions (25). The second-order implicit central finite difference scheme is used to approximate the time derivatives when solving the wave evolution problem. The physical variables are reconstructed using (25).

The following numerical study is intended to study a specific liquid metal battery proposed as a candidate for physical experiments and possible commercial implementation. The battery contains sodium metal as the anodic top layer, bismuth-sodium alloy as the bottom cathode and the liquid salt composed of NaF-NaCl-NaI (the material properties [1] are given in table 1). The cell is square in the horizontal section with the dimensions: $L_{x}\times L_{y}=0.2\times 0.2\ \text{m}^{2}$, similar as proposed in [8] for a circular cell of radius 0.2 m. The applied total electric current is fixed at the value I = 130 A ($j_{0}= 3250\ \text{A/m}^{2}$). The magnetic field typical of this cell can be estimated on the dimensional grounds: $| \bm{B}| \sim {\mu_{0}I}/ L\approx 1 \ \text{mT}$. In accordance to the previous studies [8–15] the interfacial stability is crucially affected by the vertical component of the magnetic field, therefore we choose the field as $\bm{B}=( 0,0,B_{z} )$, with $B_{z}=1\ \text{mT}$. It is assumed that the flow is purely laminar, which permits to use the shallow-layer friction coefficient derived from the gravity wave damping rate for long waves according to [16] including the bottom and the side wall friction:

$$\begin{equation} k_{fk}={( 2h_{0k}+L )\sqrt {\nu_{k}\omega_{12,(1,0)}} }/ ( 2\sqrt 2 h_{0k}L ). \end{equation} \tag{ 26 }$$

Table 1:.  Material properties of a Na∥Bi battery at the operating temperature T = 560 °C [1].

Liquid	$\rho_{k}$ (kg/m3)	$\nu_{k}$ (m2/s)	$\sigma_{k}$ (S/m)
Bi	9720	$1.1\cdot10^{-6}$	$0.69\cdot10^{6}$
NaF-NaCl-NaI	2540	$0.67\cdot10^{-6}$	200
Na	831	$0.26\cdot10^{-6}$	$3.5\cdot10^{6}$

During the charge and discharge the liquid metal levels are changing, which can affect the interfacial stability at various stages. The effect was analysed considering three different metal layer thickness combinations: the fully charged case, intermediate situation, and the final discharged case. The total liquid height is assumed to be constant: $h_{01}+h_{02}+h_{03}=8\ \text{cm}$, and the electrolyte thickness is constant: $h_{02}=1\ \text{cm}$ in all cases. A small initial wave perturbation corresponding to the long gravity wave of the mode (m = 1, n = 0) of amplitude $A=0.001\ \text{m}$ was applied at the upper interface: $H_{2}( x,y,0 )=H_{02}+A\cos( k_{x}x )\cos ( k_{y}y )$, while the lower interface was initially unperturbed, but involved in the dynamic simulation following the top perturbation. The results of the direct wave evolution problem solution are summarized in fig. 2. On the left side figs. 2(a), (c), (e) show interfacial oscillations at a fixed position in the left corner of the cell (x = 0, y = 0 in fig. 1). On the right side figs. 2(b), (d), (f) show the respective computed Fourier power spectra appearing as vertical narrow peaks. The triangles located at the horizontal axis represent the theoretical three-layer gravity frequencies defined by (11) in the absence of the electromagnetic interaction.

Zoom In Zoom Out Reset image size

When the cell is at the fully charged state the bottom metal layer is thinner than the top one: $h_{01}=2\ \text{cm}$, $h_{03}=5\ \text{cm}$. The resulting motion of the waves generated at both interfaces is stable in this case (figs. 2(a), (b)). There is a single spectral peak which corresponds to the upper metal interface gravity wave mode ((11) with "+" sign) at $f=0.542\ \text{Hz}$. The bottom metal interface oscillates at the same frequency, while the corresponding eigenmode frequency ((11) with "−" sign) is not excited. The laminar damping in this case is sufficient to suppress the instability which would otherwise appear in a purely frictionless case for the battery of square horizontal section $L_{x}=L_{y}$ in the presence of the electric current [12].

In an intermediate stage of discharge the top layer decreases while the bottom layer increases in thickness. When the metal layers are of the same thicknesses: $h_{01}=h_{03}=3.5\ \text{cm}$, the wave oscillation development is changing, see figs. 2(c), (d). In this case the unstable interfacial motion slowly sets in the upper metal. The Fourier spectrum indicates the instability onset frequency $f=0.549\ \text{Hz}$, which is the respective three-layer gravity frequency for the upper metal interface defined by (11) for the present depth values.

When the cell is approaching fully discharged state, the lower electrode is thicker than the upper one: $h_{01}=5$, $h_{03}=2\ \text{cm}$. The results shown in figs. 2(e), (f) demonstrate the unstable interfacial motion of the upper metal, which eventually leads to the destabilization of the lower metal interface as well. In this case a complete rupture of the electrolyte, giving the current short-circuit is reached after 1062 seconds. The Fourier spectrum, similarly to the previous example, contains the single low-frequency peak at the corresponding three-layer gravity wave frequency $f=0.542\ \text{Hz}$ (11). The oscillation of the lower metal is exactly at the same frequency.

From these results it can be concluded that the interfacial motion becomes unstable when the cell approaches the discharged state. These results demonstrate also that the instability is developing on the upper metal interface. This agrees with findings in [9] where the condition for the upper interface dominant role was quantified as $\omega_{1,\bm{k}}^{2}/ \omega _{2,\bm{k}}^{2}\ge 2$. In our case at $h_{01}=5$, $h_{03}=2\ \text{cm}$: $\omega_{1,(1,0)}^{2} /\omega_{2,(1,0)}^{2}=2.2\ge 2$.

The square-section cell is particularly sensitive to the magnetic interaction, as noted in [10], where the instability onset was characterised by the parameter: ${\Pi }={j_{0}B_{z}L_{x}L_{y}}/ ([ \rho_{2}-\rho_{3}]h_{02}h_{03}g )$. For the square cell the inviscid instability threshold was found at $\Pi \to 0$. The results in fig. 2 show that the viscous friction is the crucial factor between the stable and unstable states even for the critical square-section cell.

The role of the hydrodynamic coupling between the top and bottom interfaces can be further analysed if assuming that the bottom interface is fixed as $H_{1}=\textrm{const}$ at all time steps while the electric current is computed for all three layers according to (20), (21). The results for the nearly fully discharged battery ($h_{01}=5$, $h_{03}=2\ \text{cm}$) are presented in fig. 3. In this case, when the hydrodynamic coupling is neglected, the instability growth rate is even higher than in the coupled case (compare to fig. 2(e)). The respective Fourier spectra (fig. 2(f) and fig. 3(b)), give different frequency values for the oscillation at which the instability sets in. The three-layer frequency ($f=0.542\ \text{Hz}$) is lower than the corresponding two-layer frequency ($f=0.596\ \text{Hz}$).

Zoom In Zoom Out Reset image size

The role of the different types of initial perturbations can be investigated by introducing a more complex set of the initial wave shapes defined as ${H_{2}( x,y,0 )=H_{02}+A^{\prime}\sum_{m,n=0}^{4,4} {\cos ( k_{x}x )\cos ( k_{y}y )}}$ at the upper metal interface. The numerical results are shown in fig. 4, which can be compared to figs. 2(e), (f) for the single-mode perturbation. When the complex perturbation is applied, the growth of the oscillation amplitude is considerably slower compared to the case with the single-mode perturbation at the critical frequency. However the unstable wave frequencies in both cases are exactly the same ($f=0.542\ \text{Hz}$). A similar outcome was observed when the lower metal interface perturbation is additionally included. These findings indicate that regardless of the initial condition combination the instability sets in at the critical frequency determined by the upper metal interface.

Zoom In Zoom Out Reset image size

The unstable interfacial motion is illustrated in fig. 5 at consecutive time increments of a quarter of the period. It can be seen that the instability is in the form of a rotating wave along the cell perimeter. The top and bottom interfaces move in opposite phases. This observation is similar to results in [11]. Figure 6 shows the electric current distribution at the time moment when the interfaces are in the position shown in fig. 5(c). It can be seen that the horizontal current flow in both layers is in opposite directions (fig. 6(a), (b)), and the vertical current j_z distribution in the electrolyte layer correlates mostly with the upper interface due to its more pronounced deformation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The liquid metal battery of 20 cm size with sodium metal anode, which is a candidate for experimental and commercial implementation, is stable to magneto-hydrodynamic waves for certain combination of the liquid layer depths and becomes unstable at the discharge state when the top metal layer depth is reduced below a critical level. The instability does not depend on the initial interface perturbation type. It was found that when the system is close to the critical discharge state when the top metal depth is reduced to approximately 3.5 cm (equal to the bottom metal depth) it will experience the rotating wave instability unless the metal depth reduction is strictly controlled.