Numerical simulation of rolling pad instability in cuboid liquid metal batteries

TL;DR

This paper uses numerical simulations to analyze the rolling pad instability in liquid metal batteries, revealing how density differences influence stability and wave behavior, which is crucial for large-scale energy storage applications.

Contribution

It provides a detailed numerical analysis of the electromagnetic interfacial wave dynamics in liquid metal batteries, highlighting the impact of density ratios on stability and wave coupling.

Findings

Density ratio determines stability and wave type.

Instability significantly affects battery operation.

Simulation confirms earlier theoretical conclusions.

Abstract

The rolling pad instability is caused by electromagnetic interactions in systems of horizontal layers with strongly different electric conductivities. We analyze the instability for a simplified model of a liquid metal battery (LMB), a promising device for large-scale stationary energy storage. Numerical simulations of the flow and the dynamics of electromagnetically coupled interfacial waves are performed using OpenFOAM. The work confirms the earlier conclusions that the instability is a significant factor affecting battery's operation. The critical role played by the ratio between the density differences across the two interfaces is elucidated. It is found that the ratio determines the stability characteristics and the type (symmetrically or antisymmetrically coupled) of dominant interfacial waves.

Figures

b of Ref. Weber *et al.*, [2017](#bib.bib10)).

b of Ref. Weber *et al.*, [2017](#bib.bib10)).

b of Ref. Weber *et al.*, [2017](#bib.bib10)).

b of Ref. Weber *et al.*, [2017](#bib.bib10)).

b of Ref. Weber *et al.*, [2017](#bib.bib10)).

Tables

Table 1. Table 1: Physical properties of a typical Mg-Sb LMB Kim  et al. ( 2013 ) used in the simulations. Density ρ 𝜌 \rho , kinematic viscosity ν = μ / ρ 𝜈 𝜇 𝜌 \nu=\mu/\rho and electric conductivity σ 𝜎 \sigma are shown for each layer. As discussed in the text, the electrolyte density ρ E subscript 𝜌 𝐸 \rho_{E} is different from the table value in some simulations. The other parameters used in the simulations include the dimensions of the cell L x = L y = L z = 0.1 subscript 𝐿 𝑥 subscript 𝐿 𝑦 subscript 𝐿 𝑧 0.1 L_{x}=L_{y}=L_{z}=0.1 m, the unperturbed current density J 0 subscript 𝐽 0 J_{0} =7850 A/m 2 and the interfacial tension coefficients γ A − E = 0.19 subscript 𝛾 𝐴 𝐸 0.19 \gamma_{A-E}=0.19 N/m, γ B − E = 0.095 subscript 𝛾 𝐵 𝐸 0.095 \gamma_{B-E}=0.095 N/m.

	$\rho$	$\nu$	$\sigma$
unit	kg/m3	m2/s	S/m
Metal A	1577	$6.70\times {10}^{-7}$	$3.62\times {10}^6$
Electrolyte E	1715	$6.80\times {10}^{-7}$	80
Metal B	6270	$1.96\times {10}^{-7}$	$8.66\times {10}^5$

Table 2. Table 2: Summary of completed simulations. Values of the imposed magnetic field B 0 subscript 𝐵 0 B_{0} , unperturbed electrolyte thickness H E 0 superscript subscript 𝐻 𝐸 0 H_{E}^{0} , and electrolyte density ρ E subscript 𝜌 𝐸 \rho_{E} are shown. The other physical properties remain the same in all simulations (see Table I). Also shown in the table are the value of the non-dimensional parameter ( 21 ) (computed with Δ ​ ρ A Δ subscript 𝜌 𝐴 \Delta\rho_{A} in runs 1-17, and Δ ​ ρ B Δ subscript 𝜌 𝐵 \Delta\rho_{B} in runs 18-21) and the characteristics of the flow: time period of the unstable wave, stability identification, and the type of coupling (‘A’ denotes the antisymmetric coupling, while ‘S’ denotes the symmetric, in-phase coupling, both being described in detail in the text) and the typical ratio of amplitudes of the waves at the upper and lower interfaces. Two digits are shown for the amplitude ratios, but in some cases, in particular in flows with rupture, only the first digit can be reliably determined.

Case	$B_0\left[mT\right]$	$H_E^0\left[mm\right]$	${\rho }_E\left[kg/m^3\right]$	$\beta$	$T\left[s\right]$	stability	coupling	$\Delta {\eta }^A/\Delta {\eta }^B$
1	1	5	1715	0.244		stable	A	40
2	3	5	1715	0.732		stable	A	40
3	5	5	1715	1.22	3.3	saturation	A	40
4	7.5	5	1715	1.833	2.5	saturation	A	15
5	10	5	1715	2.444	3.3	saturation	A	32
6	15	5	1715	3.666	2.9	saturation	A	32
7	20	5	1715	4.888	3.3	rupture	A	40
8	26	5	1715	6.354	3.4	rupture	A	32
9	15	3	1715	5.984	4.4	rupture	A	24
10	15	7	1715	2.675	3.0	rupture	A	31
11	15	15	1715	1.496	2.3	rupture	A	38
12	15	20	1715	1.122	2.0	rupture	A	57
13	15	5	1646	7.332	4.0	rupture	A	40
14	15	5	3452.2	0.27	0.5	stable	S	1.0
15	15	3	3452.2	0.45	0.5	stable	S	1.0
16	50	3	3452.2	1.5	0.32	stable	S	1.0
17	200	3	3452.2	6	0.24	stable	S	1.0
18	7.5	5	5994	0.917	0.5/4.7	saturation	S/A	1.3
19	15	5	5994	1.833	0.5/4.6	saturation	S/A	1.2
20	30	5	5994	3.666	0.5/4.5	rupture	S/A	1.1/0.25
21	50	5	5994	6.11	0.5/4.4	rupture	S/A	1.1/0.10

Equations

$${\mathbf{J}}_0=J_0{\mathbf{e}}_y.$$

$$\sum _{i=1}^3{\alpha }_i=1,$$

$${\alpha }_A+{\alpha }_E+{\alpha }_B=1.$$

$$\rho =\sum _{i=1}^3{\alpha }_i{\rho }_i,\mu =\sum _{i=1}^3{\alpha }_i{\mu }_i,\sigma =\sum _{i=1}^3{\alpha }_i{\sigma }_i.$$

$$\frac{\partial {\alpha }_i}{\partial t}+\nabla \cdot ({\alpha }_i\mathbf{U})=0,$$

$$\frac{\partial \left(\rho \mathbf{U}\right)}{\partial t}+\nabla \cdot (\rho \mathbf{UU})=-\nabla p_d+\nabla \cdot \left[\mu (\nabla \mathbf{U}+\nabla {\mathbf{U}}^T)\right]+{\mathbf{f}}_L+{\mathbf{f}}_{st},$$

$${\mathbf{f}}_L=\mathbf{J}\times \mathbf{B}$$

$${\mathbf{f}}_{st}=-\sum _j\sum _{i\ne j}{\gamma }_{ij}\nabla \cdot \left(\frac{{\alpha }_j\nabla {\alpha }_i-{\alpha }_i\nabla {\alpha }_j}{\left|{\alpha }_j\nabla {\alpha }_i-{\alpha }_i\nabla {\alpha }_j\right|}\right)\left({\alpha }_j\nabla {\alpha }_i-{\alpha }_i\nabla {\alpha }_j\right)$$

$$\nabla \cdot \mathbf{U}=0.$$

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

$$\mathbf{J}=\sigma (-\nabla \varphi +\mathbf{U}\times \mathbf{B})$$

$$\nabla \cdot \mathbf{J}=0,$$

$$\nabla \cdot (\sigma \nabla \varphi )=\nabla \cdot (\sigma (\mathbf{U}\times \mathbf{B})).$$

$$\mathbf{U}=0\text{ at }x=0,L_x,y=0,L_y,z=0,L_z.$$

$$\frac{\partial \varphi }{\partial n}=0\text{ at }x=0,L_x,z=0,L_z$$

$$\varphi =0\text{ at }y=0,\varphi ={\varphi }_0=const\text{ at }y=L_y.$$

$${\varphi }_0=-J_0\left({\sigma }_A{H}_A^0+{\sigma }_E{H}_E^0+{\sigma }_B{H}_B^0\right).$$

$$\sum _{f=1}^n[{\sigma }_f{\mathbf{A}}_f\cdot {(\nabla \varphi )}_f]=\sum _{f=1}^n[{\sigma }_f{\mathbf{A}}_f\cdot {(\mathbf{U}\times \mathbf{B})}_f],$$

$${\mathbf{A}}_f\cdot {\mathbf{J}}_f={\sigma }_f\left[-{\mathbf{A}}_f\cdot (\nabla \varphi {)}_f+{\mathbf{A}}_f\cdot (\mathbf{U}\times \mathbf{B}{)}_f\right]$$

$${\sigma }_f={\left(\frac{\text{raise}3.01385pt{d}_N/\text{lower}3.01385ptd}{{\sigma }_N}+\frac{\text{raise}3.01385pt{d}_P/\text{lower}3.01385ptd}{{\sigma }_P}\right)}^{-1},$$

$$\Delta {\eta }^A=0.02z,$$

$$\beta =\frac{I{B}_0}{g\Delta {\rho }_A{H}_E^0{H}_A^0},$$

$$U_0={\left[g\left({\rho }_E-{\rho }_A\right){\left({\rho }_A{\left(H_A^0\right)}^{-1}+{\rho }_E{\left(H_E^0\right)}^{-1}\right)}^{-1}\right]}^{1/2}$$

$$\Delta {\eta }^A=0.00025\cos (80\pi x)\cos (80\pi z).$$

Full text

Liquid metal battery, Interfacial waves, Instability

††preprint: AIP/123-QED

I Introduction

This paper presents the results of numerical simulations of electromagnetically modified interfacial waves in a liquid metal battery (LMB) – a promising new device for large-scale stationary energy storage. The principle of the battery’s operation and the complex electrochemical, materials, technological, and economical aspects are discussed in specialized literature, e.g., in Ref. Kim et al., 2013. The following brief description and the illustration in Fig. 1 are sufficient for our analysis. The version of the battery considered in the paper can be, in a simplified way, viewed as a cuboid vessel filled with three liquid layers: heavy metalloid B (e.g., Bi, Sb, Zn, or PbSb) at the bottom, molten salt electrolyte E in the middle, and light metal A (e.g., Na, Li, Ca, or Mg) at the top. The system is stably stratified by density, with $\rho_{A}<\rho_{E}<\rho_{B}$. The electrolyte is immiscible with the metals but conductive to positive ions of the light metal. The energy stored in the battery is the difference between the Gibbs free energies of metal A in its pure state and in the state of alloy with the metal B.

During discharging, the ions of A pass through the electrolyte and form the alloy within the bottom layer. During charging, the ions are electrochemically reduced from the alloy and pass into the top layer. The sidewalls of the vessel are electrically insulated. The top and bottom walls contain current collectors for the electrons released or consumed during the charging/discharging reactions. The resulting vertical currents flowing through the battery are quite strong, with the density $J\sim 10^{4}$ A/m2. The system is maintained at a temperature above the melting point of the three materials, which is facilitated by the strong Joule heating of the electrically poorly conducting electrolyte.

The system presents significant challenges related to hydrodynamic effects, especially when transition from small-scale laboratory prototypes to economically preferable large-scale ($\sim 1$ m3) devices is considered (see Ref. Kelley and Weier, 2018 for a review). One such challenge is the need for mixing of reactants in the bottom layer. Others are identified as hydrodynamic instabilities expected in large batteries. When considering the instabilities, it is important to realize that if the deformation of the electrolyte-metal interfaces is strong enough to bring the two metals into contact with each other, it creates short circuit between the metal layers and, thus, requires immediate stop of the battery’s operation. Any instability leading to such a deformation must be avoided.

Several instability mechanisms have been identified in recent studies: the Tayler instability,Seilmayer et al. (2012); Weber et al. (2014); Herreman et al. (2015) thermal convection caused by the Joule heating of the electrolyteShen and Zikanov (2016) or by bottom heating,Kelley and Sadoway (2014) the Marangoni instability at the electrolyte-metal interfaces,Köllner, Boeck, and Schumacher (2017) and the rolling pad (also called metal pad) instability.Zikanov (2015); Weber et al. (2017); Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018); Herreman (2019) The first three instabilities, while interesting and definitively affecting the battery’s dynamics, are not expected to result in the just described disruption of operation. EstimatesHerreman et al. (2015); Shen and Zikanov (2016) show that, at realistic battery parameters, the forces leading to the instability are balanced by the gravity force due to the density stratification between the layers so that the growing perturbations saturate at a low amplitude of interface deformation. The situation is less clear in the case of the rolling pad instability considered in this paper.

The mechanism of the rolling pad instability is known from the operation of the Hall-Héroult aluminum reduction cells, each of which can, in a schematic and drastically simplified way, be presented as the cell in Fig. 1 minus the top metal layer (see, e.g., Ref. Davidson, 2016). Many substantial features of the Hall-Héroult cell differentiate it from the liquid metal battery: two layers instead of three, a structured top current collector, the effects of gas bubbles and turbulent flow of the melts, low aspect ratio at very large (several m) horizontal dimensions, and so on. Both the systems, however, have the same key feature necessary for the rolling pad instability, namely strong electric currents passing through layers of vastly different electric conductivities. The conductivity of molten salt electrolytes is about four orders of magnitude lower than the conductivity of liquid metals. This means that even a small-amplitude deformation of the metal-electrolyte interface causing a small-amplitude variation of the local electrolyte thickness along the current path leads to significant perturbations of currents and, thus, to significant additional Lorentz forces acting on the melts.

For the Hall-Héroult cells, it has been long understood that the Lorentz forces may create melt flows enhancing the interface deformation and leading to an instability in the form of a growing interfacial wave rotating around the cell (the “rolling pad”). The key mechanism has been identified as the interaction between the perturbation currents primarily closing within the metal due to its high electric conductivity and the vertical component of the magnetic field always present in the cell.Sele (1977); Urata (2016); Bojarevics and Romerio (1994); Sneyd and Wang (1994); Davidson and Lindsay (1998); Zikanov et al. (2000); Sun, Zikanov, and Ziegler (2004); Sun, Zikanov, and Finlayson (2005)

The possibility of a similar scenario in the case of liquid metal battery has been a subject of active research recently. The first approachZikanov (2015) was based on the mechanical analogy, in which sloshing of two metal layers was modeled by oscillations of two rectangular solid pendula suspended in an electrolyte. The main result was a demonstration, albeit in the framework of a strongly simplified model, that the rolling pad instability mechanism is realizable in a three-layer system. This conclusion was confirmed using a much more realistic model in the three-dimensional simulations of a cylindrical battery cell.Weber et al. (2017) The instability was developing at the upper interface, where the density jump was much smaller than at the lower one. It was concluded that the non-dimensional instability parameter derived for the Hall-Héroult cells Sele (1977) was also a good indicator of the instability in the battery system. Further analysis was performed using two-dimensional shallow water models.Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018) While limited in their applicability to systems with low aspect ratio and low-to-moderate amplitudes of interface deformation, such models provided the advantage of easier numerical solution and stability analysis than for the three-dimensional model and allowed the authors to study various aspects of the instability. The two other recent studies remain to be mentioned: the linear stability analysis,Molokov (2018) where the critical role of the non-conducting sidewalls was elucidated, the theoretical and numerical analysis,Horstmann, Weber, and Weier (2018) where the hydrodynamic coupling between the waves at the two interfaces was analyzed, and the perturbation theory analysis of waves in a cylindrical cell with viscous, Joule dissipation, and capillary effects.Herreman (2019)

This paper presents the results of the three-dimensional simulations of the rolling pad instability in a liquid metal battery of cubic shape. As in the other studies of the instability,Zikanov (2015); Weber et al. (2017); Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018); Herreman (2019) a simplified physical model is used, in which the instability effect is separated from the effects of thermal convection, Marangoni, and Tayler instabilities. The model is presented in detail in section II. The numerical method described in section III is based on the OpenFoamOpe and, in its many features, similar to the method applied in Ref. Weber et al., 2017. The simulation procedure is described in section IV. The results presented in section V are quite different from those in Ref. Weber et al., 2017 in that they pertain to a cuboid rather than a cylindrical cell (a substantial difference for sloshing flows), extend the analysis to systems of various ratios between the density differences across the interfaces, and apply a broader array of diagnostics tools thus putting the study into the context established by the recent works.Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018)

II Physical model

The liquid metal battery is a complex system. Its operation involves various physical mechanisms, some of which are unknown. Therefore, a reasonable first step toward understanding the physics of the battery is to investigate each mechanism individually, isolated from the rest, and in the framework of a simplified model. This approach has been applied to the Tayler instability,Seilmayer et al. (2012); Weber et al. (2014); Herreman et al. (2015) thermal convection,Shen and Zikanov (2016); Kelley and Sadoway (2014) Marangoni instability,Köllner, Boeck, and Schumacher (2017) and the rolling pad instability.Zikanov (2015); Weber et al. (2017); Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018); Herreman (2019)

We apply the approach to the rolling pad instability. The simplified physical model is, in its many aspects, similar to that used in Ref. Weber et al., 2017 and based on the following main assumptions.

Simplified geometry

The geometry of the battery is simplified in our case to a cuboid with electrically perfectly insulating sidewalls and the top and bottom walls serving as equipotential current collectors (see Fig. 1a).

Simplified model of magnetic field

The typical magnetic field within an operating battery is three-dimensional and possibly unsteady. A part of the field, which is primarily horizontal, is induced by the currents flowing through the battery. The other part is truly three-dimensional and is generated by the currents flowing in the supply lines and neighboring battery cells. By the analogy with aluminum reduction cellsSele (1977); Urata (2016); Bojarevics and Romerio (1994); Sneyd and Wang (1994); Davidson and Lindsay (1998); Zikanov et al. (2000); Sun, Zikanov, and Ziegler (2004) and based on the recent results concerning the liquid metal batteries,Zikanov (2015); Weber et al. (2017); Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018); Herreman (2019) we assume that the rolling pad instability is caused by the vertical component of the magnetic field. For this reason, we neglect all components but the vertical one, which is approximated by a constant. In addition to removing the effects of the Tayler instability and the electro-vortex flow, this assumption also means that we ignore the background flow of the melts that would be induced by the interaction between the variable component of the electric current and the three-dimensional magnetic field even when the interfaces are flat. Only the electromagnetic interactions directly participating in the rolling pad instability are retained in the model.

Absence of solid phase

During the discharge stage, poor mixing at insufficiently high temperature can cause formation of intermetallic solid in the bottom layer just below the interface.Kelley and Weier (2018) This effect is ignored in our study under the assumption that proper care is taken to avoid this phenomenon generally detrimental for battery’s operation.

Immiscibility

The liquid metals forming the top and bottom layers of the liquid metal battery are assumed to be immiscible with the electrolyte forming the middle layer.

Constant temperature

The liquid metals and the electrolyte are assumed to be maintained at the same constant temperature. Accordingly, the effects of thermal and Marangoni convection are ignored. The approximation is justified by the recent resultsShen and Zikanov (2016); Köllner, Boeck, and Schumacher (2017) showing that under typical operational conditions the thermal and Marangoni convection are unlikely to perturb the interfaces strongly and to have any significant effect beyond additional mixing within each layer.

Constant fluid properties

Each fluid is assumed to be incompressible, Newtonian and having a constant electrical conductivity. The properties are, however, different among the three fluids. The coefficients of surface tension at the two interfaces are assumed to be constant.

Ignoring mass transport through interfaces

The typical time scale of the processes analyzed in this paper is tens to hundreds of seconds, which is much smaller than the typical time of charging or discharging the battery. This allows us to ignore the change of mass and physical properties of the metal layers due to the possible cross-interface diffusion and the transport of ions through the electrolyte. The only manifestation of the ion motion is the vertical electric current.

Quasi-static approximation of electromagnetic fields

Using the expected typical velocity $U$ below $0.1$ m/s, the typical size $L=0.1$ m, and the typical maximum electric conductivity $\sigma=4\times 10^{6}$ S/m, we estimate the magnetic Reynolds number as $\mbox{{Re}}_{m}\equiv UL\sigma\mu_{0}=4.5\times 10^{-2}$ (here $\mu_{0}$ is the magnetic permeability of vacuum). Even at larger $UL$, the number remains substantially smaller than one. We can apply the quasi-static approximation, according to which the magnetic field induced by fluid motion is neglected in comparison with the imposed field in the expressions for the Ohm’s law and Lorentz force (see, e.g., Ref. Davidson, 2016).

Further simplifications, such as the shallowness of the systemZikanov (2018); Tucs, Bojarevics, and Pericleous (2018) or small amplitude of perturbationsMolokov (2018); Herreman (2019) are not made in our model. Similarly to Ref. Weber et al., 2017, the problem is solved as three-dimensional, non-linear, and time-dependent.

The assumptions made above imply that in the unperturbed state, in which the interfaces between the layers are perfectly horizontal, the electric current is purely vertical and has constant density:

[TABLE]

As will be evident in the following discussion of the governing equations, the body forces and, thus, the melt velocity are identically zero in this case. A deformation of the interfaces causes a variation of the local thickness of the electrolyte, which leads to current perturbations $\bm{j}$, Lorentz forces, and melt flows (see Fig. 1b).

The governing equations are based on the one-fluid approach in its volume-of-fluid version (see, e.g., Ref. Prosperetti and Tryggvason, 2009). This is an Eulerian method, in which the kinematic and dynamic conditions at the interface are implemented on a fixed grid via direct balances of fluxes of mass, momentum, and, in the case of our system, electric charge. The presence of each of the three liquids at a given location is determined by its phase fraction $0\leq\alpha_{i}(\bm{x},t)\leq 1$, which is defined as the volume fraction of the liquid $i$ in a computational cell. Only two phase fractions are independent, since at every location we have

[TABLE]

or, in terms of our notation,

[TABLE]

The physical properties of the mixture are defined as weighted averages. For the density, dynamic viscosity and electrical conductivity, we have, respectively,

[TABLE]

The mass conservation in the absence of diffusion and mass transfer between the phases requires that each phase fraction satisfies

[TABLE]

where $\bm{U}$ is the velocity of the mixture assumed to be shared by all fluids.

The momentum equations is

[TABLE]

where

[TABLE]

is the Lorentz force discussed in detail below, and

[TABLE]

is the interface tension force approximated as a volumetric force active in the interface area,Brackbill, Kothe, and Zemach (1992) with $\gamma_{ij}$ being the constant coefficient of interfacial tension.

The velocity is subject to the incompressibility constraint:

[TABLE]

The constraint is satisfied in the standard manner via solution of the Poisson equation for the modified pressure $p_{d}=p+\rho gy$.

As we have already mentioned, only the vertical component of the magnetic field is included:

[TABLE]

The electric current density is found using the Ohm’s law

[TABLE]

and the charge conservation condition

[TABLE]

which results in the elliptic equation for the electric potential $\phi$:

[TABLE]

The electric current can be viewed as the sum of the base current (1) and the current perturbations caused by the deformation of the interfaces and the flow of conducting fluids in a constant magnetic field.

The no-slip boundary conditions are imposed at the walls:

[TABLE]

The conditions on the electric potential include those of perfect insulation of the sidewalls:

[TABLE]

and those of equipotential top and bottom surfaces:

[TABLE]

The potential drop is calculated for the unperturbed state of the battery with constant thicknesses $H_{A}^{0}$, $H_{E}^{0}$, $H_{B}^{0}$ of the layers as

[TABLE]

The condition (15) is an approximation based on the assumption that the perturbations of the electric current form loops closing entirely within the battery cell. The same condition was used in Ref. Weber et al., 2017. It is slightly different from the condition applied, e.g., in Ref. Zikanov, 2018, where the potential drop was adjusted continuously to account for the change of the total electric resistance of the cell. The difference only appears at strong deformations of the interfaces and has no critical impact on the instability even in that case.

III Numerical method

Time-dependent three-dimensional solutions of the problem are computed using the finite-volume CFD tool OpenFOAM.Ope The method is similar to that applied to the rolling pad instability in a cylindrical battery cell in Ref. Weber et al., 2017. It is based on the solver for simulations of incompressible multiphase flows multiphaseInterFoamOpe ; Ubbink (1997); Deshpande, Anumolu, and Trujillo (2012) modified to include electromagnetic fields and account for sharp variations of the electrical conductivity.

The finite-volume technique is described in detail in, e.g., Ref. Deshpande, Anumolu, and Trujillo, 2012, so only a brief description is given here. The spatial discretization is of the second order. The first-order simple implicit scheme is applied for time discretization. The PISO algorithmIssa (1986) is applied to iteratively solve the momentum and pressure equations. The advection equation for the phase fraction (4) is solved using the implicit Euler time discretization and the spatial discretization by the MULES schememarquez (2013) based on the upwind interpolation outside the interface areas and a delimiter-based combination of the linear and upwind interpolations inside it. The use of the compressive fluxes allows the solver to restrict the numerically diffused interface thickness to a few (typically two or three in our simulations) cell widths.

The discretization of the electromagnetic equations is based on the conservative scheme.Ni et al. (2007a, b) The potential equation (12) integrated over the finite volume cell $\Omega$ with faces $\bm{A}_{f}$, $f=1,\ldots,n$ is written as

[TABLE]

where the subscript $f$ indicates the values at the midpoint of the cell faces $f$, $(\nabla\phi)_{f}$ is the face-normal gradient evaluated by central differences, and ${({\bm{U}}\times{\bm{B}})}_{f}$ is computed using linear interpolation. The face-normal components of the electric current are calculated at the cell faces using the Ohm’s law as

[TABLE]

and then interpolated linearly to the cell’s centers to be used in the evaluation of the Lorentz force (6).

The extreme (four orders of magnitude) difference between the electrical conductivities of the electrolyte and the metals presents a challenge to accurate calculation of the electric potential and currents, in particular to satisfying the charge conservation condition (11) in the interface area. We have determined that the difficulties are resolved by two adjustments. One is the use of harmonic interpolation for the evaluation of the electric conductivity of the multiphase mixture $\sigma$ (see (3)) at the cell faces in (17) and (18):

[TABLE]

where $P$ and $N$ are the two neighboring cell centers, $d=|PN|$ is the distance between them, and $d_{N}$ and $d_{P}$ are the distances from cell centers to the midpoint $f$ of the cell face separating them from each other ($d_{N}=|fN|$, $d_{P}=|fP|$). As discussed in the context of heat conduction in Ref. Patankar, 1980, the harmonic interpolation avoids the unacceptably large error in the evaluation of fluxes at cell faces in the situation of strong variation of the transport coefficient within the cell.

We have also found that satisfying the charge conservation and avoiding unacceptably large spurious electric currents near the interfaces requires convergence of the iterative solution of (17) to very low tolerance, at the level of $10^{-20}$. In the simulations presented in section V of this paper, this typically required about 500 iterations.

The numerical model was verified by repeating several results presented in the earlier three-dimensional simulations of the rolling pad instability in a cylindrical battery.Weber et al. (2017) Consistently good agreement was found. In particular, we were able to accurately reproduce the period of interfacial waves in an unstable system and the value of the threshold magnetic field above which the system experiences short circuit (see Fig. 10b of Ref. Weber et al., 2017).

IV Simulation procedure

IV.1 Problem setting and diagnostics

The interior of the battery is modeled as a cube of side length $L_{x}=L_{y}=L_{z}=0.1$ m (see Fig. 1). This is approximately the size of the currently tested commercial LMB prototypes, but substantially smaller than the size ($\sim 1$ m) of the envisioned large-scale batteries. Various values of the unperturbed thickness of the electrolyte layer $H_{E}^{0}$ are used in the simulations. The layer is always located in the middle of the cell, so the unperturbed thicknesses of the metal layers are $H_{A}^{0}=H_{B}^{0}=\left(L_{y}-H_{E}^{0}\right)/2$.

The square horizontal cross-section of the battery cell differentiates our work from the recent studies, where cells of circularWeber et al. (2017) or rectangularZikanov (2018); Tucs, Bojarevics, and Pericleous (2018) cross-sections were considered. The analogy with the Hall-Héroult reduction cells leads us to expect the instability characterized by (i) the threshold much lower than for rectangular cells, and (ii) possible spatial complexity and multitude of instability modes.Davidson and Lindsay (1998) We will return to these expectations in the following analysis.

The physical properties of the materials corresponding to the Mg-Sb battery with NaCl-KCl-MgCl2 electrolyteKim et al. (2013) (see table 1) are used in some of the simulations. This relatively well studiedZikanov (2015); Weber et al. (2017); Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018) combination of material properties is an example of the system with small density difference across the upper interface $\Delta\rho_{A}\equiv\rho_{E}-\rho_{A}\ll\Delta\rho_{B}\equiv\rho_{B}-\rho_{E}$. In order to see how the instability changes in batteries with other ratios between $\Delta\rho_{A}$ and $\Delta\rho_{B}$, the electrolyte density $\rho_{E}$ is varied in the other simulations, while the rest of the properties are kept the same.

The imposed vertical current $J_{0}=7850$ A/m2 is used in all the simulations. The value $\phi_{0}$ of the potential at the top wall is computed according to (16). The positive value of $J_{0}$ implies that the process of battery charging is simulated. Based on the physical mechanism established in earlier studiesZikanov (2015); Weber et al. (2017); Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018) and confirmed in our work, the instability considered in the framework of our simplified model is expected to have the same behavior in the course of charging, except for the opposite direction of the precession of interfacial waves.Zikanov (2018)

Each simulation starts with the initial state, in which $\bm{U}=0$, $\bm{J}=J_{0}\bm{e}_{y}$, and $\phi$ is given by the piecewise-linear function of $y$ derived as a solution of (12) with the boundary conditions (14)–(16). A weak perturbation of the electrolyte thickness is added to initiate the flow. It has the shape

[TABLE]

which implies the maximum amplitude of 2 mm. Each simulation is continued forward in time until the behavior of the flow becomes clear.

To investigate the effect of the initial conditions, we have repeated several simulations of unstable systems with different shapes of the initial perturbations as well as with the initial perturbations applied to the bottom rather than top interface. As discussed and illustrated in section V.3, these changes affect the first stage of the flow evolution, but not the qualitative and quantitative characteristics of the fully developed waves.

The solutions are visualized using the snapshots of the interfaces, distributions of velocity and Lorentz force fields. The profiles of the upper and lower interfaces $y=\eta^{A}(x,z,t)$ and $y=\eta^{B}(x,z,t)$ are determined during the simulations as the locations where, respectively, $\alpha_{B}=0.5$ and $\alpha_{A}=0.5$. We also use the distributions of the vertical electric current density $J_{y}(x,z,t)$ in the mid-plane of the electrolyte layer $y=0.05$ m, which is nearly reciprocal to the local thickness of the electrolyte.

The quantitative characteristics of the interfacial waves (period, growth or decay of amplitude) are determined from the time signals of $\eta^{A}$, $\eta^{B}$ and of the velocity and electric current components at several points in the metal and electrolyte layers. The points are chosen along the vertical line at $x$=0.05 m, $z=0.016$ m, although any other line would give similar characteristics.

As in the previous simulations of the instability,Weber et al. (2017); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018) the problem is solved in dimensional units. This is justified by the multitude of relevant physical properties and the general desire to keep the analysis focused on the realistic values of these properties. At the same time, it is useful and interesting to consider the possible forms of the non-dimensional parameters controlling the instability. From the three-dimensional analysis of the cylindrical cellsWeber et al. (2017) and the shallow water analysis of rectangular cellsZikanov (2018) we conclude that in the case of small density difference across the upper interface the key parameters are likely to be the shape of the horizontal cross-section of the cell and the parameterSele (1977)

[TABLE]

where $I$ is the total current through the battery ($=J_{0}L_{x}L_{z}$ in our system). This non-dimensional parameter was originally introduced by SeleSele (1977) to characterize the metal pad rolling in aluminum reduction cells. It estimates the typical ratio between the destabilizing Lorentz force and the stabilizing effect of density stratification. A similar parameter with $\Delta\rho_{B}$ in place of $\Delta\rho_{A}$ can be formally introduced in the case of small density difference across the lower interface. There is no reason to expect the parameter to accurately describe the instability in the cases where the density differences across the two interfaces are comparable to each other.

The role of the other non-dimensional parameters, such as the Reynolds, Hartmann, and Weber numbers remains largely unknown. We will leave this question to future studies, but note that the expected large typical length scale of growing perturbations suggests that these parameters somewhat affect the growth rate and the shape of the unstable modes, but do not change them in a profound way.

The results reported in this paper can be recalculated into non-dimensional form using, for example, the velocity of a large-scale gravitational wave on the upper interface

[TABLE]

as the typical velocity scale, and $L_{x}$, $L_{x}/U_{0}$, $\rho_{E}U_{0}^{2}$, $J_{0}$, and $J_{0}L_{x}\sigma_{A}^{-1}$ as the typical scales of length, time, pressure, electric current density, and electric potential, respectively.Zikanov (2018)

IV.2 Computational grid and grid sensitivity tests

Constant time step and a structured Cartesian grid are used in all the simulations. The design of the grid is based on the following considerations. Firstly, we expect particularly strong vertical gradients of velocity, electric potential, and other variables in the areas around the interfaces, especially when the interfaces are significantly deformed. Secondly and more importantly, accuracy of determining the position of the interfaces is critical for the overall accuracy of the model. An error in the evaluation of the local relative change of $H_{E}$ leads to an approximately $\sigma_{A}/\sigma_{E}\sim\sigma_{B}/\sigma_{E}\gg 1$ times bigger error in the evaluation of the perturbation currents and, thus, of the Lorentz forces that drive the instability. An accurate volume-of-fluid simulation of the interface movement requires fine resolution in the direction perpendicular to the interface. Since long waves are expected to dominate in the flows created by the rolling pad instability, this means a requirement of fine resolution in the vertical direction.

Constant grid steps $\Delta x=\Delta z$ are used in the horizontal plane. In the $y$-direction, a constant small grid step $\Delta y_{min}$ is used within the horizontal layer that extends on both sides sufficiently far beyond the area of expected locations of the two interfaces. The position of this layer is not necessarily symmetric with respect to the middle of the cell, since in some of the simulations, only one interface is expected to move significantly. In the metal domains below and above this layer, $\Delta y$ increases in geometric progression with the expansion rate chosen so that it reaches $\Delta y_{max}=1$ mm at the bottom and top boundaries of the cell.

Grid sensitivity tests were performed to determine the optimal values of $\Delta x=\Delta z$, $\Delta y_{min}$, and the time step $\Delta t$. The grid was deemed sufficient if further refinement did not result in significant changes of the time period and peak values of the signals of velocity and electric current density. A detailed description of the results can be found in Ref. Xiang, 2018. In summary, the spatial discretization with $\Delta y_{min}=0.2$ mm and $\Delta x=\Delta z=2$ mm is sufficient for accurate simulations (see an example in Fig. 2). The situation with the effect of the time step $\Delta t$ is more complex. We have found that, at the spatial grid described above, $\Delta t=4\times 10^{-3}$ s is sufficient to avoid numerical instabilities and accurately predict the velocity field. At the same time, the signals of the vertical electric current $J_{y}$ inside and near the electrolyte, which is very sensitive to the accuracy of reproduction of interface movement, only become truly $\Delta t$-independent at $\Delta t$ below $2\times 10^{-3}$ s. Furthermore, $\Delta t$ must be decreased in strongly unstable regimes, when the amplitudes of the interface deformation, electric current perturbations, velocities and velocity gradients become large, and the electrolyte layer is about to be ruptured. The strategy in the simulations was to use $\Delta t=4\times 10^{-3}$ s, while periodically conducting verification runs at $\Delta t=2\times 10^{-3}$ s and reducing $\Delta t$ in the cases of rapidly growing perturbations.

V Results

V.1 Summary of completed simulations

The completed simulations are summarized in Table II. Each case is identified by the values of the parameters varied in the study: the imposed vertical magnetic field $B_{0}$, the unperturbed thickness of electrolyte $H^{0}_{E}$, and the electrolyte density $\rho_{E}$. We also show the values of the non-dimensional parameter $\beta$ (see (21)) and the main period of the wave if it can be determined from the time signals of the flow variables. The last three columns show the characteristics of the flow regimes introduced later in the text: the stability identification, the type of wave coupling, and the typical ratio between the wave amplitudes at the upper and lower interfaces.

V.2 Examples of flow behavior

In this section, we show three examples typifying the three principal scenarios possible in the system. The electrolyte density is $\rho_{E}=1715$ kg/m3 (a Mg-Sb battery), which corresponds to $\Delta\rho_{A}\ll\Delta\rho_{B}$, and allows us to consider the solutions in the context of earlier studies.Weber et al. (2017); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018)

The first example is the flow at $H_{E}^{0}=5$ mm, $B_{0}=1$ mT, and $\rho_{E}=1715$ kg/m3 (case 1 in table II). We see in Fig. 3 that the amplitudes of the interfacial waves created by the initial perturbation decay with time. Similar decay is observed for all the other flow variables, such as the kinetic energy, perturbations of electric currents, etc. This allows us to identify the system as stable.

The second example illustrated in Figs. 4 and 5 is the flow of the case 5, in which the system parameters are the same as in the case 1, except that the strength of the vertical magnetic field $B_{0}$ is increased to 10 mT. Fig. 4 clearly shows that there are two stages of the flow evolution. The first stage, lasting for about 60 s, is characterized by growth of the fluctuation amplitudes of all variables. During the second stage, which continues for the rest of the simulated flow evolution, the amplitudes remain approximately constant. The situation is identifiable as that of an instability leading to nonlinear saturation at some level of finite-amplitude perturbations.

The quantitative characteristics of the interfacial wave in the saturated flow regime can be determined from the plots in Fig. 4 and other simulation data not shown here. The dominant period of the wave is about 3.3 s. The amplitude of the fluctuation is slightly less than 2.5 mm for the upper interface, and much smaller, about 0.07 mm for the lower interface. As we see in, respectively, Figs. 4b and 4c, the resulting variations of the local thickness of electrolyte result in strong (up to 5000 A/m2 in amplitude) fluctuations of the vertical current and strong (up to 2000 A/m2) horizontal currents in the metal layer. The melt flows associated with the wave are not strong. The typical magnitude of the horizontal velocity components is $\sim 1$ cm/s in the electrolyte and smaller than that in the metal layers (see Figs. 4d-f). It is interesting that the velocity in the electrolyte fluctuates around zero, while the signals of velocity in the metal layers show additional strong variations on the typical time scale of many tens of s. Our analysis of the evolution of the velocity structure indicates that this behavior is a reflection of a slow transformation of the rotating pattern of horizontal velocity associated with the wave.

We need to stress that the just mentioned fluctuation amplitudes refer to the signals recorded at the specific locations used in Fig. 4. Signals recorded at other points would produce different values. The order of magnitude and the general characteristics of the wave, such as the lengths of growth and saturation stages, time period, and the qualitative shape of the signal curves would, however, remain the same.

The structure of the flow with a fully developed wave (the saturation stage) is illustrated in Fig. 5. The interfacial waves presented as three-dimensional (in Fig. 5a) and two-dimensional (in Figs. 5c and d) plots of $\eta^{A}(x,z)$ and $\eta^{B}(x,z)$ are qualitatively similar to the typical waves observed in the two layer systems (aluminum reduction cells)Zikanov et al. (2000); Sun, Zikanov, and Ziegler (2004); Lukyanov, El, and Molokov (2001) and three-layer systems with $\Delta\rho^{A}\ll\Delta\rho^{B}$.Weber et al. (2017); Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018); Herreman (2019) The waves are characterized by the wavelength comparable with the horizontal size of the cell and the zone of large-amplitude interface deformation located near the sidewall. The pattern rotates clockwise if viewed from above (see animationsXiang (2019)). The rotation is responsible for the nearly periodic oscillations of $\eta^{A}$, $\eta^{B}$ and other variables shown in Fig. 4. Two interesting and useful for our further discussion observations can be made. One is that the wave at the upper interface has much larger amplitude than the lower interface wave. Another observation is that the two waves are antisymmetrically (with the phase shift of 180 degree) coupled.

The distribution of the vertical electric current $J_{y}$ in the electrolyte (see Fig. 5b) follows closely the variation of the local electrolyte thickness, thus the high-amplitude interfacial wave at the upper interface. The distributions of other variables, such as flow velocity (see Fig. 5e) or Lorentz force (see Fig. 5f) cannot be explained that simply but are evidently parts of the flow associated with the electromagnetically coupled interfacial waves.

In addition to causing local current perturbations, the variation of the electrolyte thickness caused by interfacial waves also changes the overall impedance $R_{LMB}$ of the battery. This affects the energy efficiency. In our model simulations, in which the voltage drop across the battery remains constant, the impedance changes are manifested as the reciprocal variations of the total (integrated over a horizontal cross-section) vertical current $J_{LMB}$. Furthermore, the mass conservation of liquid phases implies that the integral of the perturbation of electrolyte thickness and, thus, of the perturbation of local electric resistance are zero. One should expect decreased overall impedance and increased total current in this case.

An example shown in Fig. 6 illustrates the typical behavior found in the simulated flows with instability. The time-averaged overall current is about 89.6 A, which is noticeably larger than 78.5 A for the unperturbed battery. The interfacial wave significantly reduces the overall impedance of the cell. We also see fluctuations correlating with the evolution of the wave. This is to be expected for a wave of finite amplitude, which changes its shape in the course of its progression around the battery.

We note that in the case 5 as in all the other cases identified in table II as those with saturation, the interface deformation never grows as strong as to cause a rupture of the electrolyte layer and short circuit between the metal layer. The system is hydrodynamically unstable, but this does not result in operational failure of the battery.

Our third example is the case 7 in table II. The system’s parameters are the same as in the cases 1 and 5, but the vertical magnetic field $B_{0}$ is now increased to 20 mT. We see in Fig. 7 that the system is unstable and that the outcome of the instability is different from that observed in the previous example. The perturbations do not saturate but continue to grow until the electrolyte layer ruptures at some point. This implies zero overall impedance of the electrolyte layer and a short circuit between the metal layers, which makes further operation of the battery impossible. The situation can be identified as that of not just hydrodynamic, but also operational instability.

The structure of the flow about 0.5 s before the rupture is illustrated in Fig. 8. We see current perturbations, velocity and velocity gradients of the amplitudes comparable to those in case 5. The situation changes shortly before the rupture. At illustrated in Fig. 9, the interfaces nearly touch each other in a small area in the corner of the cell (the corner location was observed in all our simulations with rupture). The vertical electric current grows to about $10^{6}$ A/m2 in this corner, which results in very strong localized horizontal current perturbations and Lorentz forces.

At the moment of rupture, the vertical current increases further. The numerical model fails to converge, so the simulations must be stopped. It should be stressed that during the short time interval before the rupture, the small-scale, strong-gradient features appearing in the flow imply that the accuracy of the model’s prediction is lower than estimated in our grid sensitivity studies. Additional simulations performed on a refined grid and with substantially reduced time step produced a qualitatively and quantitatively similar flow. The rupture event was delayed by a just fraction of a second. This indicates that the electrolyte rupture observed in the simulations is a real physical effect and not an artifact of the numerical model.

V.3 Effects of system parameters

In this section we analyze how the system’s behavior is affected by the key parameters of the rolling pad instability: the strength of the vertical magnetic field $B_{0}$, the electrolyte thickness $H_{E}^{0}$, and the electrolyte density $\rho_{E}$. The discussion is based on the data summarized in Table II and the illustrations in Figs. 3-14.

The following discussion will benefit from a reference to the recent analysis of interfacial waves in a stably-stratified three-layer system.Horstmann, Weber, and Weier (2018) Flow in a cylindrical vessel without electromagnetic effects was considered. Two possible solution modes were identified. One was the fast mode with high oscillation frequency and symmetric (in-phase) coupling between the waves on the upper and lower interfaces. Another was the slow mode with lower frequency and antisymmetric (with 180∘ degree phase shift) coupling. The solutions of the two types also differ from each other by the typical ratios between the amplitudes of the upper and lower interface waves.

The system considered in this paper is evidently very different from that of Ref. Horstmann, Weber, and Weier, 2018 due to the different geometry and the effects of viscosity and surface tension. The electromagnetic effects provide an additional coupling and forcing mechanisms and may influence the modes themselves and the selection between them in a particular solution. Nevertheless, as we will see below, the classification remains valid and useful in the analysis.

We start with the cases 1-13, in which $\Delta\rho_{A}\ll\Delta\rho_{B}$. The three flows used as examples in section V.2 (see Figs. 3-9) belong to this group. Further information is provided in Table II and Fig. 10. We clearly see that the solutions have the form of a slow mode in the classification of Ref. Horstmann, Weber, and Weier, 2018. The typical time period $T$ is large in comparison to the other cases. The waves at the two interfaces are coupled antisymmetrically (see Figs. 3, 4, 7, 10). Finally, the wave amplitude on the upper interface is much (20 to 60 times) larger than the amplitude of the lower interface wave. This is in agreement with the results of calculations for the material properties of a Mg-Sb cell presented in Fig. 4 of Ref. Horstmann, Weber, and Weier, 2018. For the typical horizontal wavenumber $\sim 20m^{-1}$ (observed in our simulations as illustrated in Figs. 5 and 8), the slow modes are expected to have $\Delta\eta^{A}$ one to two orders of magnitude larger than $\Delta\eta^{B}$. The growth of $\Delta\eta^{A}/\Delta\eta^{B}$ with increasing $H_{E}^{0}$ (see Table II) is also consistent with the behavior of the slow mode solutions illustrated in Fig. 4 of Ref. Horstmann, Weber, and Weier, 2018.

The combination of the material properties in the cases 1-12 corresponds to the Mg-Sb battery, which is relatively well studied in literature.Weber et al. (2017); Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018) It has been found (see, e.g. Refs. Weber et al., 2017; Zikanov, 2018) that the value of the parameter $\beta$ defined by (21) is an important, albeit not unique, factor of the instability in a cell of a given shape of the horizontal cross-section. The data in Table II are partially consistent with this conclusion. As we increase $B_{0}$, while keeping $H_{E}^{0}=5$ mm and $\rho_{E}=1715$ kg/m3, the state of the system changes from stable at approximately $\beta$ less than one (cases 1 and 2) to saturated at $\beta$ up to about 4 (cases 3-6) and then to ruptured (cases 7 and 8). The threshold values of $\beta$ are reasonably close to those of Ref. Weber et al., 2017 with the difference attributable to different cell geometries (a cylinder on Ref. Weber et al., 2017 and a cube in our case).

The situation with the effect of the electrolyte thickness $H_{E}^{0}$ is less clear. Five simulations (cases 6, 9-12) are performed at $B_{0}=15$ mT, $\rho_{E}=1715$ kg/m3, and $H_{E}^{0}$ increasing from 3 to 20 mm. This corresponds to decrease of $\beta$, and, if $\beta$ is the parameter controlling the instability, should lead to less unstable behavior. This is observed in the transition from case 9 ($H_{E}^{0}=3$ mm, $\beta=5.984$, rupture) to case 6 ($H_{E}^{0}=5$ mm, $\beta=3.666$, saturation). At larger $H_{E}^{0}$, however, rupture of the electrolyte layer is observed (cases 10-12), even though increase of $H_{E}^{0}$ within this group leads to slower growth of the unstable wave (compare the signals for cases 10 and 12 in Fig. 10). We note that special care has been given to assure adequate grid resolution of the flow in these cases. For example, in the case 12, the fine mesh with $\Delta y=\Delta y_{min}=2\times 10^{-4}$ m was used at $0.02\leq y\leq 0.08$. This can be compared with the highest (just before the rupture) position of the upper interface at about $y=0.0668$. We conclude that the most plausible explanation of the behavior in the cases 6, 9-12 is that $\beta$ is an imperfect instability parameter for such systems.

The analysisHorstmann, Weber, and Weier (2018) also shows that, for solutions with $\Delta\eta^{A}\gg\Delta\eta^{B}$, the behavior of the wave on the upper interface is close to that of the wave developing in a two-layer system. The presence of the lower interface has only small impact. The fact that such solutions are found in our and otherWeber et al. (2017); Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018); Herreman (2019) studies of the Mg-Sb cells is, therefore, consistent with the detected at least partial role of $\beta$ derived for a two-layer aluminum reduction cell as the instability parameter.

In the cases 14-17, the electrolyte density $\rho_{E}$ is such that the density differences across the lower and upper interfaces are the same. The system is stable in all these cases even though $\beta$ can be as large as 6 (in case 17). No consistent growth of perturbations is found. Instead, the perturbations decay initially and then the flow shows nearly periodic oscillations of approximately constant small amplitude. An illustration in the form of the signals of interface locations is provided in Fig. 11. We see that the waves are coupled symmetrically (in-phase with each other) and have much smaller time period (between 0.24 and 0.5 s) than the waves observed in the cases 1-13. This allows us to identify them as the fast modes according to the classification of Ref. Horstmann, Weber, and Weier, 2018. We also see in Fig. 11 that the amplitudes of the upper and lower interfacial waves are nearly equal to each other.

The spatial structure of the flow is illustrated in Fig. 12 as spatial distributions of $H_{E}(x,z,t)$ on upper ($\eta^{A}$) and lower interface ($\eta^{B}$). The combination of symmetric coupling and nearly equal oscillation amplitudes implies that the local thickness of the electrolyte layer is nearly constant. Some variations of $H_{E}(x,z,t)$ and $J_{y}(x,z,t)$ on the small length scale are found, but they are of small and not growing amplitude. We conclude that the fast mode solutions realized in our system at $\Delta\rho_{A}=\Delta\rho_{B}$ do not generate significant Lorentz forces and, thus, cannot lead to the rolling pad instability.

Another feature of the interfacial waves in flows with $\Delta\rho_{A}=\Delta\rho_{B}$ is that, unlike in the unstable flows at $\Delta\rho_{A}\ll\Delta\rho_{B}$, the shape of the wave remains strongly correlated with the shape imposed as an initial perturbation. We see in Fig. 12 that, apart from effect of rotation and development of a wave at the lower interface, the shape resembles the initial constant slope wave 60 s after the start.

The just reported findings can be compared with the results of the two-dimensional model of shallow battery cells of rectangular cross-section,Zikanov (2018) to our knowledge the only other study where the case $\Delta\rho_{A}=\Delta\rho_{B}$ is considered. Despite the evident differences in the geometry and modeling approach, the results of the two studies are in clear quantitative agreement with each other. Stable solutions in the form of symmetrically coupled waves of approximately equal amplitudes are found in both studies.

The situation when $\Delta\rho_{A}\gg\Delta\rho_{B}$ is represented by the cases 18-21. As illustrated in Fig. 13, the flow can be unstable. Growth of perturbations occurs and ends in saturation in the flow with $\beta$ (based on $\Delta\rho_{B}$) equal to 0.917 (case 18) and 1.833 (case 19). Growth leading to rupture of electrolyte occurs at $\beta=3.666$ (case 20) and 6.11 (case 21). The analysis of the distributions of flow variables and time signals, such as those in Fig. 13, leads us to the conclusion that the solutions are superpositions of fast and slow modes. The fast modes cause the high-frequency (period of about 0.5 s) symmetric oscillations with approximately equal amplitudes on both interfaces. The slow modes are responsible for low-frequency (period between 4.4 and 4.7 s) antisymmetric oscillations. Their amplitudes cannot be accurately determined from the data in Fig. 13. For this reason only the amplitude ratios of the fast mode oscillations are listed in Table II. Nevertheless we clearly see that the amplitude of the slow mode is significantly higher on the lower than on the upper interface. Interestingly, the slow modes are only noticeably strong when the perturbations are growing, i.e. during the initial growth stage in the saturation cases 18 and 19 and during the entire flow evolution in the rupture cases 20 and 21. During the saturation phase in cases 18 and 19, the amplitude of the slow mode oscillations gradually decreases until they practically disappear.

The typical spatial structure of the flows in the cases 18-21 is illustrated in Fig. 14. We see that the upper ($\eta^{A}$) and lower interface ($\eta^{B}$) are coupled and oscillated symmetrically while the amplitude of oscillation differs from upper and lower interface. This implies variable local electrolyte thickness and, thus, Lorentz forces and the possibility of unstable behavior.

We conclude that in the simulations with $\Delta\rho_{A}\gg\Delta\rho_{B}$ (cases 18-21), the rolling pad instability is caused by the slow modes developing at the lower interface. Unlike the flows with $\Delta\rho_{A}\ll\Delta\rho_{B}$, however, the unstable waves evolve in the presence of strong fast mode oscillations.

We will now discuss the effect of the form of the initial perturbations. All the simulations presented so far were initialized by a slight incline of the upper interface (20). The simulations show that shape of the wave is retained by stable flows, in particular by the flows with $\Delta\rho_{A}=\Delta\rho_{B}$ (cases 14-17). The situation is quite different for unstable flows, in which the initial conditions only affect the early evolution. The final fully developed state is determined by the instability mechanism.

Our first illustration is the test simulations in which the initial perturbation of the same shape are applied to the lower rather than upper interface. This does not change the initial distribution of the Lorentz forces, but affects the initial potential energy of the perturbations, specifically its distribution between the waves on the upper and lower interfaces and its total magnitude, which is proportional to the density difference across the perturbed interface. One may justifiably ask whether this would affect the selection between fast and slow modes or affect the solution in some other way.

The results of the numerical experiments carried out for cases 5 and 19 are shown in Fig. 15. We see that the change of the initial conditions does not influence the solution in any major way.

In the second test, the perturbations are imposed at the upper interface but have a different shape, smaller horizontal length scale, and smaller (with the maximum of 0.25 mm) amplitude:

[TABLE]

The results obtained for the case 5 are presented in Fig. 16. We see that the instability takes much longer to develop at the new initial conditions. At the fully developed stage, however, the parameters of the wave, such as its period, amplitude, and the spatial shape (compare Figs. 16b and 5c) remain practically the same.

VI Concluding remarks

The simulation results reported in this paper confirm the conclusion reached earlier in the studies of other geometriesWeber et al. (2017); Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018); Herreman (2019) that the rolling pad instability is a significant factor affecting performance of liquid metal batteries. In particular, batteries with much smaller density difference across the upper interface ($\Delta\rho_{A}\ll\Delta\rho_{B}$), for example the better studied Mg-Sb battery, become unstable and face the danger of a ruptured electrolyte layer if the vertical component of the magnetic field is too large. The instability is caused by growing antisymmetrically coupled interfacial waves, which can be identified as the slow modes in the classification of Ref. Horstmann, Weber, and Weier, 2018. The wave amplitude is much larger on the upper than on the lower interface in such solutions, and the behavior is qualitatively similar to the behavior of unstable waves in two-layer aluminum reductions cells.

Much less expected and, to our knowledge, entirely new results are found for the much less explored configurations in which $\Delta\rho_{A}=\Delta\rho_{B}$ or $\Delta\rho_{A}\gg\Delta\rho_{B}$. We see that the analogy with the aluminum reduction cells would be misleading for such systems. In addition to the slow modes, the symmetrically coupled fast modesHorstmann, Weber, and Weier (2018) appear in the solution. At $\Delta\rho_{A}=\Delta\rho_{B}$, only the fast modes are present. Having nearly equal amplitudes, the interfacial waves do not cause significant variations of electrolyte thickness and, therefore, cannot cause the rolling pad instability. The interfacial waves that appear in batteries with $\Delta\rho_{A}\gg\Delta\rho_{B}$ are superpositions of fast and slow modes, with the latter causing the instability.

The results are consistent with the findings of earlier works,Weber et al. (2017); Horstmann, Weber, and Weier (2018); Zikanov (2018); Tucs, Bojarevics, and Pericleous (2018); Molokov (2018); Herreman (2019) in which batteries of other geometries (cylinders or shallow parallelepipeds) were considered. Specifically, qualitative agreement is found with Refs. Weber et al., 2017; Horstmann, Weber, and Weier, 2018; Zikanov, 2018; Tucs, Bojarevics, and Pericleous, 2018; Molokov, 2018; Herreman, 2019 for the case $\Delta\rho_{A}\ll\Delta\rho_{B}$ and with Ref. Zikanov, 2018 for the case $\Delta\rho_{A}=\Delta\rho_{B}$.

Our parametric study is incomplete, since only a small sample of possible combinations of material properties could be covered. Furthermore, comparison of our results with the results of earlier studies confirms the rather evident conclusion that the interfacial wave solutions are critically affected by the shape of the horizontal cross-section of the battery cell. An exhaustive parametric study of the rolling pad instability would be unjustifiably difficult and costly to complete. Analyzing specific shapes and material combinations of the battery versions selected for commercialization appears to be a more promising path of future research.

A practically relevant numerical model of a working battery cell would require accurate description of a large number of diverse physical effects. In addition to the effects addressed in this paper, the list includes the thermal and concentration convection, diffusive and convective mass transport, surface tension, Tayler instability, electrovortex effect, etc. Full three-dimensional magnetic field has to be included. No such comprehensive simulations have been done so far. The question whether such a simulation is possible at the current level of computer modeling techniques remains open. The OpenFOAM seems a promising tool due to the abundance of already developed and verified solver modules and the relative easiness of implementation of physical models. Our work and the earlier studiesWeber et al. (2017); Horstmann, Weber, and Weier (2018) suggest that one computationally particularly challenging task, namely the description of the interface motion and the associated electric currents, can be accomplished with the help of the OpenFOAM. Acceptable, albeit not perfect accuracy is achieved on a properly refined grid if appropriate modifications of the discretization scheme (see section III) are made.

Acknowledgements.

The authors are thankful to Norbert Weber for help with the numerical method and to Norbert Weber, Tom Weier, Valdis Bojarevics, and Gerrit Horstmann for interesting and stimulating discussions. Financial support was provided by the US NSF (Grants CBET 1435269 and CBET 1803730) and the University of Michigan - Dearborn.