# Layer coupling between solutal and thermal convection in liquid metal
batteries

Paolo Personnettaz111The two first authors contributed equally to the
publication. Tanja Sophia Klopper222The two first authors contributed equally
to the publication. Norbert Weber  Tom Weier  Helmholtz-Zentrum Dresden –
Rossendorf, Bautzner Landstr. 400, 01328 Dresden, Germany  Carl von Ossietzky
University Oldenburg, Ammerländer Heerstr. 114-118, 26129 Oldenburg, Germany

###### Abstract

For longer than one decade, liquid metal batteries (LMBs) are developed with
the primary aim to provide economic stationary energy storage. Featuring two
liquid metal electrodes separated by a molten salt electrolyte, LMBs operate
at elevated temperature as simple concentration cells. Therefore, efficient
mass transfer is a basic prerequisite for their economic operation.
Understanding these mechanisms cannot be limited at the single layer level.
With this motivation, the effects of solutal- and thermally-driven flow are
studied, as well as the flow coupling between the three liquid layers of the
cell. It is shown that solutal convection appears first and thermal convection
much later. While the presence of solutal flow depends on the mode of
operation (charge or discharge), the occurrence of thermal convection is
dictated by the geometry (thickness of layers). The coupling of the flow
phenomena between the layers is intriguing: while thermal convection is
confined to its area of origin, i.e. the electrolyte, solutal convection is
able to drive flow in the positive electrode and in the electrolyte.

###### keywords:

liquid metal battery , mixing , solutal convection , Rayleigh-Bénard
convection , internally heated convection , layer coupling

††journal: International Journal of Heat and Mass Transfer

thermal convection appears first in the electrolyte – at charge and discharge

solutal convection emerges much faster and stronger than thermal convection

solutal convection can mix all layers by viscous coupling, thermal convection
not

solutal convection affects internally heated convection in the electrolyte

coupled top interface allows for better mixing in the positive electrode than
no-slip

##  1 Introduction

Since the turn of the millennium, liquid metal batteries (LMBs) are discussed
as low-cost stationary electrical energy storage for fluctuating renewable
sources. Figure 1 illustrates the typical setup of such a device. The negative
electrode on top is separated from the positive electrode by a thin molten
salt electrolyte layer [1]. In order to prevent a short-circuit between the
electrodes, violent fluid flow needs to be avoided [2].

Figure 1: Setup and operation of a Li$||$Bi LMB at charge (a) and discharge
(c) as well as the typical vertical temperature profile during operation (b).
Sub-figure (d) illustrates a charge-discharge cycle with typical effects due
to mass-transport overvoltage [3].

LMBs are simple concentration cells. At discharge, the electro-active species
migrates from the negative electrode through the electrolyte before being
alloyed into the positive electrode. It has already been noticed during the
first days of LMB research that non-ideal mixing during this alloying process
leads to mass transport overpotentials. A highly concentrated solute region
located near the interfaces will develop especially at high current density
[4]. In extreme cases, even solid intermetallic phases might form at the
electrolyte-positive electrode interface [5, 6]. Similarly early, thermal
convection [6] and mechanical stirring [5, 7] have been proposed to mitigate
concentration overpotentials.

In recent years, many flow effects have been investigated aiming to enhance
mass transfer in liquid metal electrodes. Electro-vortex flow (EVF), which is
driven by converging or diverging current lines, has been explored by Ashour
et al. [8]. Shortly later, it has been shown that the placement of the feeding
lines as well as stray magnetic fields have a critical influence on the flow
structure and magnitude of EVF [9]. Using an energy balance and numerical
modelling, Herreman et al. concluded that non-swirling electro-vortex flow is
too weak to break the stable density stratification during discharge inside
the positive electrode [10].

Recently, Weber et al. combined convection with an electrochemical model
showing that the influence of EVF on the cell voltage is indeed limited. On
the other hand, a vertical magnetic background field, which drives strong
swirling flow, is able to increase the cell efficiency substantially [11]. The
latter has been confirmed theoretically and with DNS simulations by Herreman
et al. as well [12].

Thermal convection in LMBs has been subject of various investigations.
Motivated by mixing the positive electrode alloy, Kelley & Sadoway heated a
PbBi layer from below measuring the flow velocity [13]; the experiment was
later reproduced numerically by Beltrán et al. [14]. Shen & Zikanov were the
first who modelled thermal convection in a complete LMB. They show that it is
caused mainly by the heat released in the electrolyte, as illustrated in
Figure 1b. Consequently, flow is expected to appear in the negative electrode
and electrolyte, while the positive electrode is stably stratified. Flow
structures and magnitude depend very much on the thickness of the single
layers [15]. Köllner et al. [16] confirmed these early investigations and
pointed out the presence of three additional thermally driven flows: Marangoni
flows at both interfaces and anti-convection in the lower layer.

Finally, Personnettaz et al. discussed the presence of electrochemical heat,
and highlighted that in realistic LMBs with thin electrolyte layers (5 mm)
thermal convection will be dominant in the thick negative electrode, but
negligible in the electrolyte itself [17].

While thermal convection and EVF appear independently of the direction of the
current, solutal convection can be observed only at charge. As illustrated in
Figure 1a for a Li$||$Bi cell, Li is transferred then from the positive
electrode alloy towards the negative electrode. The heavy (Li-poor) Bi-Li
alloy will sink down in plumes. At discharge, this process is reversed: Li is
alloyed on top of the Bi-rich layer forming a stable density stratification
(Figure 1c). This “asymmetry” of the effects leads to an “asymmetry” within
the corresponding charge-discharge curve, shown in Figure 1d. At discharge a
considerable mass-transport overpotential is observed, which is not present
while charging. When switching off the current, the voltage returns
immediately to the equilibrium potential after a charging phase (solutal
flow), but needs a certain diffusive relaxation time after a discharge phase
(stable density stratification). The crucial relevance of solutal convection
and the stable density stratification has to our best knowledge first been
discussed by Kelley & Weier [2] and has for the first time been evidenced by
Personnettaz et al. [3]. The findings have later been confirmed by Herreman et
al. [10] performing two- and three-dimensional numerical simulations of mass
transport and EVF. Interestingly, the two mass-transfer polarization effects,
highlighted in Figure 1d can also be observed in older measurements. For
example, the charge-discharge curve of a Na$||$Bi cell from 2015 shows an
additional overpotential at discharge [18], and Fig. 8b in Kim et al. [19] the
relaxation of the cell potential after a discharge cycle. Tab. 1 gives a few
more examples of related effects in the old literature.

Table 1: Experimental results showing mass transport overvoltage in the literature of liquid metal electrodes. Chronological order. observation | source
---|---
polarization effects at high discharge current in Li-Cd and Li-Zn cells | [20, p. 202]
clear concentration polarization, when switching off current after discharge | [21, p. 75]
strong asymmetry of interface concentration between charge and discharge | [22, p. 955]
clear concentration polarization effects after switching off the current | [23, p. 159]
charge-discharge curve shows additional (concentration) losses at discharge | [18]

The aim of this article is to elucidate the interaction between solutal and
thermal convection in particular during the transient start-up of the battery.
In this context, we pay special attention to the coupling of the flow between
the three distinct layers. In the first part of the article, we state the
relevant equations that describe thermal and solutal convection in three-layer
LMBs. From them, we deduce the relevant dimensionless number and we provide
estimations of the expected orders of magnitude for different chemistries.
Previous studies are then used to compute the critical current density at
which thermal convection appears in the two top layers. Then, we move to the
numerical modelling of 2D Li$||$Bi cells during charge and discharge. There,
we discuss the effect of sign and magnitude of the current on heat, mass and
momentum transfer.

##  2 Theory

The distinct regions $\Omega^{(i)}$ (N-negative electrode, E-electrolyte,
P-positive electrode) of an LMB are coupled at their interfaces
$\partial\Omega^{(i|j)}$, see Figure 3. They permit the exchange of momentum,
heat and – for certain ions – mass through electrochemical reactions. The
transport processes inside the regions depend on the boundary conditions and
material properties. These vary widely from layer to layer and so do the
transport phenomena. In this section we will discuss the equations describing
these mechanisms, we will rearrange the relevant parameters in the form of
dimensionless numbers and provide indications on the presence of convection in
the different layers.

###  2.1 Mass transport

Mass transport is vital for the operation of electrochemical cells. In the
negative electrode of LMBs, concentration gradients are typically absent
because it consists in many cases of a pure substance (but see Ca-Bi$||$Ca-Sb
[24], Ca-Mg$||$Bi [25] and Li-Mg$||$Sb-Pb [26] for counterexamples). In this
work, we neglect the presence of advective and diffusive transport of ions in
the electrolyte. Assuming a perfectly blended electrolyte seems justified by
experimental evidence [27]. Therefore, concentration distributions are
considered only in the positive electrode, assuming a binary mixture A-B (in
the present study Li-Bi). The mass concentration of component A,
$\rho_{\mathrm{A}}$ ($\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$), is described by
the advection-diffusion equation [28, 3, 10]

$\frac{\partial\rho_{\mathrm{A}}}{\partial
t}+\mathbf{u}\cdot\nabla\rho_{\mathrm{A}}=\mathcal{D}_{\mathrm{AB}}\nabla^{2}\rho_{\mathrm{A}}\
\mathrm{in}\ \Omega^{\mathrm{(P)}},$ (1)

in which $\mathbf{u}$ is the velocity field and $\mathcal{D}_{\mathrm{AB}}$ is
the mass transport diffusion coefficient. We neglect the volume change of the
positive electrode due to the density variation and concentration dependent
mass transport coefficient [29]. As the solid boundaries are impermeable, we
impose an homogeneous Neumann condition for the concentration. The interface
between the electrolyte and the positive electrode allows for the transfer of
mass, because ions are reduced or oxidised there as
$\mathrm{A}^{+}+\mathrm{e}^{-}\leftrightharpoons\mathrm{A(inB)}$ through
electrochemical reactions. The mass flux of the component A under amperostatic
condition, $\dot{m}^{\prime\prime}_{\mathrm{A}}$, can be computed with the
Faraday law of electrolysis. The resulting boundary condition can be
formulated then as

$-\mathcal{D}_{\mathrm{AB}}\nabla\rho_{\mathrm{A}}\cdot\textbf{n}=\dot{m}^{\prime\prime}_{\mathrm{A}}=\frac{j\mathcal{M}_{\mathrm{A}}}{n_{\mathrm{el}}\mathrm{F}}\
\mathrm{at}\ \partial\Omega^{\mathrm{(\mathrm{P|E})}}\ ,$ (2)

where $j$, $\mathcal{M}_{\mathrm{A}}$, n, $n_{\mathrm{el}}$, $\mathrm{F}$ are
current density, molar mass of the component A (Li in our application), the
surface normal vector, the number of charges exchanged by the electrochemical
reaction and the Faraday constant, respectively. We consider only the primary
current distribution, which is uniform across all the domain, assuming that
the top and bottom plate are ideal current collectors and that the lateral
walls are insulated.

###  2.2 Heat transport

In LMBs, heat is mainly generated by the passage of the current through the
electrolyte, and then transported into the full cell. For an overview of
thermal phenomena in LMB and the relevant assumptions, please refer to
Personnettaz et al. [17]. The transport equation for the temperature $T$ in
the $i$th layer reads

$\frac{\partial T}{\partial t}+\mathbf{u}\cdot\nabla
T=\mathcal{D}^{(i)}_{T}\nabla^{2}T+\frac{\mathcal{D}_{\mathrm{T}}^{(i)}}{k^{(i)}}(\dot{q}^{\prime\prime\prime})^{(i)}\
\mathrm{in}\ \Omega^{(i)}\ \mathrm{for}\ i=\mathrm{N,E,P},$ (3)

in which $\mathcal{D}_{\mathrm{T}}$ and $k$ are thermal diffusivity and
conductivity. The ohmic internal heat generation is described by
$\dot{q}^{\prime\prime\prime}=\rho_{\mathrm{el}}j^{2}$, in which
$\rho_{\mathrm{el}}$ is the electrical resistivity. Electrochemical heat
generation is completely neglected. The top and bottom boundary are held at
uniform temperature $T_{\mathrm{ref}}$. Lateral sides are assumed adiabatic.
At the interfaces between the layers, $\partial\Omega^{(i|j)}$, continuity of
temperature and heat flux must be enforced. Most of these assumptions are
shared by previous works [17, 16, 15].

###  2.3 Fluid Flow

In LMBs we have three immiscible liquid phases; however, the flow in each
layer is characterised by an independent topology. We work in the assumption
of fixed interfaces and of the Boussinesq approximation, similar as other
authors [16, 15]. The momentum balance is defined as

$\frac{\partial\mathbf{u}}{\partial
t}+\mathbf{u}\cdot\nabla\mathbf{u}=-\frac{1}{\rho_{\mathrm{ref}}^{(i)}}\nabla
p_{\mathrm{d}}+\nu^{(i)}\nabla^{2}\mathbf{u}\\\
+(\beta_{\mathrm{T}}^{(i)}(T-T_{\mathrm{ref}})+\beta_{\rho_{\mathrm{A}}}^{(k)}(\rho_{\mathrm{A}}-\rho_{\mathrm{A,ref}}))\mathbf{g}\
\mathrm{in}\ \Omega^{(i)}\ \mathrm{for}\ i=\mathrm{N,E,P}\ \mathrm{and}\
k=\mathrm{P},$ (4)

where $\rho_{\mathrm{ref}}$ is the density in the reference thermodynamic
condition, $\nu$ is the kinematic viscosity, $p_{\mathrm{d}}$ is the modified
pressure, $\mathbf{g}$ is the gravity acceleration vector and
$\beta_{\mathrm{T}}$ and $\beta_{\rho_{\mathrm{A}}}$ are thermal and solutal
expansion coefficient [30]. We neglect the presence of other body forces, like
the one generated by magneto-hydrodynamic effects. The velocity field
satisfies the continuity equation in the incompressible form
$\nabla\cdot\mathbf{u}=0$. At the interfaces, $\partial\Omega^{(i|j)}$, the
normal velocity is enforced to be zero. The tangential components are assumed
continuous across the interfaces. Furthermore, the stress balance

$(\rho_{\mathrm{ref}}\nu)^{\mathrm{(i)}}\frac{\partial u_{k}}{\partial
z}=(\rho_{\mathrm{ref}}\nu)^{\mathrm{(E)}}\frac{\partial u_{k}}{\partial z}\
\mathrm{at}\ \partial\Omega^{(i|\mathrm{E})}\ \mathrm{for}\ i=\mathrm{N,P}\
\mathrm{and}\ k=x,y$ (5)

must be satisfied. In this context, solutal and thermal Marangoni effects are
neglected. The walls are rigid and impermeable and described by a no-slip
boundary condition.

###  2.4 Dimensionless numbers

The system of equations that describes convective heat and mass transport in a
3-layer LMB contains 14 different material properties, the reference state
($T_{\mathrm{ref}}$, $\rho_{\mathrm{A,ref}}$), and one operating parameter:
the current density $j$. Furthermore, we have to consider the geometrical
dimensions: three layer heights $\Delta h^{\mathrm{(i)}}$ and the lateral
extension of the cell $L_{x}$. This large number of parameters motivates
performing a dimensionless analysis. Here, we consider every layer per se and
we use as example the positive electrode in which both heat and mass transfer
take place. At least five distinct time scales can be defined: three based on
diffusion ($\mathcal{T}_{\mathcal{D}i}=\mathcal{L}^{2}\mathcal{D}_{i}^{-1}$)
and two based on the free fall velocity (e.g.
$\mathcal{T}_{\beta\Theta}=\mathcal{L}(g\beta_{\mathrm{T}}\Delta\Theta\mathcal{L})^{-0.5}$),
in which $\mathcal{D}_{i}$, $\mathcal{L}$ and $\Delta\Theta$ are the
diffusivity of the $i$-th component, the reference height and temperature
difference, respectively. If we substitute them in the system of equations,
the following five dimensionless numbers appear [28]:

$\displaystyle\text{Prandtl number:\qquad}Pr$
$\displaystyle=\frac{\nu}{\mathcal{D}_{\mathrm{T}}}$
$\displaystyle=\frac{\mathcal{T}_{\mathcal{D}\Theta}}{\mathcal{T}_{\mathcal{D}u}},$
(6) $\displaystyle\text{Schmidt number:\qquad}Sc$
$\displaystyle=\frac{\nu}{\mathcal{D}_{\mathrm{AB}}}$
$\displaystyle=\frac{\mathcal{T}_{\mathcal{D}\mathcal{C}}}{\mathcal{T}_{\mathcal{D}u}},$
(7) $\displaystyle\text{Lewis number:\qquad}Le$
$\displaystyle=\frac{\mathcal{D}_{\mathrm{T}}}{\mathcal{D}_{\mathrm{AB}}}$
$\displaystyle=\frac{\mathcal{T}_{\mathcal{D}\mathcal{C}}}{\mathcal{T}_{\mathcal{D}\Theta}},$
(8) $\displaystyle\text{Thermal Rayleigh number:\qquad}Ra_{\Theta}$
$\displaystyle=\frac{g\beta_{\mathrm{T}}\Delta\Theta\mathcal{L}^{3}}{\mathcal{D}_{\mathrm{T}}\nu}$
$\displaystyle=\frac{\mathcal{T}_{\mathcal{D}u}\mathcal{T}_{\mathcal{D}\Theta}}{\mathcal{T}_{\beta\Theta}^{2}},$
(9) $\displaystyle\text{Solutal Rayleigh number:\qquad}Ra_{\mathcal{C}}$
$\displaystyle=\frac{g\beta_{\rho_{\mathrm{A}}}\Delta\mathcal{C}\mathcal{L}^{3}}{\mathcal{D}_{\mathrm{AB}}\nu}$
$\displaystyle=\frac{\mathcal{T}_{\mathcal{D}u}\mathcal{T}_{\mathcal{D}\mathcal{C}}}{\mathcal{T}_{\beta\mathcal{C}}^{2}},$
(10) $\displaystyle\text{Density ratio:\qquad}R_{\rho}$
$\displaystyle=\frac{\beta_{\rho_{\mathrm{A}}}\Delta\mathcal{C}}{\beta_{\mathrm{T}}\Delta\Theta}$
$\displaystyle=\frac{\mathcal{T}_{\beta\Theta}^{2}}{\mathcal{T}_{\beta\mathcal{C}}^{2}}.$
(11)

Here, $\Delta\mathcal{C}$ denotes the reference concentration difference
across the layer. The first three dimensionless numbers are ratios of couples
of diffusivities and related time scales. The two Rayleigh numbers quantify
the strength of the buoyancy with respect to dissipative effects of viscosity
and component diffusion. Lewis number and density ratio are not essential
parameters of the problem, but they are helpful to understand the relative
importance of heat and mass transport.

The concentration difference across the positive electrode can be estimated a
priori using the boundary conditions of Eq. 2:

$\Delta\mathcal{C}_{\mathrm{apriori}}=\frac{\dot{m}^{\prime\prime}_{\mathrm{A}}\Delta
h^{\mathrm{(P)}}}{\mathcal{D}_{\mathrm{AB}}}=\frac{j\mathcal{M}_{\mathrm{A}}\Delta
h^{\mathrm{(P)}}}{n_{\mathrm{el}}\mathrm{F}\mathcal{D}_{\mathrm{AB}}}\ .$ (12)

This value largely overestimates the real difference, as the concentration
gradient is present only in a thin region near the interface. A better
estimation can be done a posteriori, applying a 1D diffusion model, using the
semi-infinite approximation:

$\Delta\mathcal{C}_{\mathcal{D}}(t)=\frac{j\mathcal{M}_{\mathrm{A}}}{n_{\mathrm{el}}\mathrm{F}\mathcal{D}_{\mathrm{AB}}}\sqrt{\frac{4\mathcal{D}_{\mathrm{AB}}t}{\pi}}.$
(13)

This formulation is time dependent; for the analytical solution of the
diffusion problem, please refer to [3] and [10]. If solutal convection takes
place, the variation across the layer is even lower and we can compute the
latter a posteriori by numerical simulations.

The temperature differences within each layer can be estimated applying a 1D
pure conduction model, like the ones presented by Köllner et al. [16] and
Personnettaz et al. [17]. Here we use a simplified model based on the
assumption that the thermal resistances ($\frac{\Delta
h^{\mathrm{(i)}}}{k^{\mathrm{(i)}}}$) of the two liquid metal electrodes are
comparable. In this condition, the maximum of the temperature is situated at
the mid height of the electrolyte, where we expect no vertical heat flux.
Therefore, we can split the domain in two in order to decouple the electrodes.
The temperature jump across the full height of the $i$-th electrode can be
expressed as

$\Delta\Theta^{(i)}=\frac{\rho_{\mathrm{el}}^{\mathrm{(E)}}j^{2}\Delta
h^{\mathrm{(E)}}}{2}\frac{\Delta h^{(i)}}{k^{(i)}}\ \mathrm{for}\
i=\mathrm{N,P}.$ (14)

If additional heat fluxes are present at the interfaces, due to
electrochemical reactions, it is sufficient to add them to the first term. In
the electrolyte the maximum temperature variation spans across half of the
electrolyte and can be estimated as

$\Delta\Theta^{\mathrm{(E)}}=\frac{\rho_{\mathrm{el}}^{\mathrm{(E)}}j^{2}(\Delta
h^{\mathrm{(E)}})^{2}}{8k^{\mathrm{(E)}}}.$ (15)

This quantity was used as the temperature unit by Köllner et al. [16]. From
these three simple formulas, we deduce that the temperature difference will
increase quadratically with the current density, but the concentration
difference only linearly. We can use Eq. 14 and Eq. 13 to estimate thermal and
solutal Rayleigh numbers in the liquid metal electrodes at a fixed current.
Tab. 2 gives an overview on the estimated dimensionless numbers for different
binary mixtures of the positive electrode. In calculating these numbers, we
assumed a layer height $\Delta h$ of $1\text{\,}\mathrm{cm}$ for the
electrodes and the electrolyte, a current density of
$1\text{\,}\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$ and electrical resistivity
of the salt of
$5\text{\cdot}{10}^{-3}\text{\,}\mathrm{\SIUnitSymbolOhm}\text{\,}{\mathrm{m}}^{-1}$.
Most of the material has a Prandtl number substantially lower than 1. The
diffusion of heat in the electrodes is very efficient. It becomes clear that
the solutal Rayleigh number is at least two orders of magnitude larger than
the thermal Rayleigh number. This suggests that solutal convection will be
able to mix the thermally stably stratified positive electrode layer easily
during charge. On the other hand, flow induced in the positive electrode by
viscous coupling with the electrolyte is expected to be strongly dampened by
the compositionally stably stratified layer during discharge. Two exceptions
are electrodes made of Selenium or Tellurium. As both have a considerably
lower thermal diffusivity compared to other liquid metals, we expect in such
electrodes an interesting competition between heat and mass transfer.

Table 2:  Heat and mass transport dimensionless numbers for different combinations of active materials in the positive electrode. For the sources of the material properties, refer to [31, 32]. The properties of the alloys (A-B) are computed at molar concentration of $10\text{\,}\mathrm{\char 37\relax}$ of component A. Due to the scarcity of mixture properties we use blending rules, as described in §4. The mass transport diffusivity is estimated with the model proposed by Roy and Chhabra [33]. $Ra_{\mathrm{T,1}}$, $Ra_{{\rho_{A}},1}$ are computed considering a $\Delta\Theta=$1\text{\,}\mathrm{K}$$, $\Delta\mathcal{C}=$1\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$$, respectively. All values of Ra are given as absolute values. During discharge, $Ra_{{\rho_{A}}}$ is negative, Eq. (13) is evaluated at $t=$$10\text{\,}\mathrm{s}$. Due to the stable temperature distribution, $Ra_{\mathrm{T}}$ is negative as well. name | T | Pr | Sc | Le | $Ra_{\mathrm{T,1}}$ | $Ra_{\mathrm{T}}$ | $Ra_{{\rho_{A}},1}$ | $Ra_{{\rho_{A}}}$ | $R_{\rho}$
---|---|---|---|---|---|---|---|---|---
Ca-Bi | $550\text{\,}$ | $1.3\text{\cdot}{10}^{-02}\text{\,}$ | $1.6\text{\cdot}{10}^{01}\text{\,}$ | $1.2\text{\cdot}{10}^{03}\text{\,}$ | $7\text{\cdot}{10}^{02}\text{\,}$ | $1\text{\cdot}{10}^{03}\text{\,}$ | $4\text{\cdot}{10}^{06}\text{\,}$ | $3\text{\cdot}{10}^{08}\text{\,}$ | $9\text{\cdot}{10}^{00}\text{\,}$
Ca-Sb | $700\text{\,}$ | $1.7\text{\cdot}{10}^{-02}\text{\,}$ | $4.0\text{\cdot}{10}^{01}\text{\,}$ | $2.3\text{\cdot}{10}^{03}\text{\,}$ | $4\text{\cdot}{10}^{02}\text{\,}$ | $5\text{\cdot}{10}^{02}\text{\,}$ | $4\text{\cdot}{10}^{06}\text{\,}$ | $4\text{\cdot}{10}^{08}\text{\,}$ | $20\text{\,}$
K-Hg | $250\text{\,}$ | $3.6\text{\cdot}{10}^{-03}\text{\,}$ | $8.6\text{\cdot}{10}^{00}\text{\,}$ | $2.4\text{\cdot}{10}^{03}\text{\,}$ | $1\text{\cdot}{10}^{04}\text{\,}$ | $2\text{\cdot}{10}^{04}\text{\,}$ | $1\text{\cdot}{10}^{08}\text{\,}$ | $4\text{\cdot}{10}^{10}\text{\,}$ | $20\text{\,}$
K-Tl | $250\text{\,}$ | $1.8\text{\cdot}{10}^{-02}\text{\,}$ | $1.6\text{\cdot}{10}^{02}\text{\,}$ | $9.2\text{\cdot}{10}^{03}\text{\,}$ | $3\text{\cdot}{10}^{02}\text{\,}$ | $3\text{\cdot}{10}^{02}\text{\,}$ | $2\text{\cdot}{10}^{07}\text{\,}$ | $8\text{\cdot}{10}^{09}\text{\,}$ | $40\text{\,}$
Li-Bi | $450\text{\,}$ | $1.3\text{\cdot}{10}^{-02}\text{\,}$ | $1.9\text{\cdot}{10}^{01}\text{\,}$ | $1.5\text{\cdot}{10}^{03}\text{\,}$ | $1\text{\cdot}{10}^{03}\text{\,}$ | $2\text{\cdot}{10}^{03}\text{\,}$ | $2\text{\cdot}{10}^{07}\text{\,}$ | $6\text{\cdot}{10}^{08}\text{\,}$ | $1\text{\cdot}{10}^{00}\text{\,}$
Li-Cd | $493\text{\,}$ | $9.2\text{\cdot}{10}^{-03}\text{\,}$ | $6.8\text{\cdot}{10}^{01}\text{\,}$ | $7.4\text{\cdot}{10}^{03}\text{\,}$ | $3\text{\cdot}{10}^{02}\text{\,}$ | $1\text{\cdot}{10}^{02}\text{\,}$ | $2\text{\cdot}{10}^{07}\text{\,}$ | $1\text{\cdot}{10}^{09}\text{\,}$ | $7\text{\cdot}{10}^{00}\text{\,}$
Li-Pb | $500\text{\,}$ | $1.5\text{\cdot}{10}^{-02}\text{\,}$ | $2.1\text{\cdot}{10}^{01}\text{\,}$ | $1.5\text{\cdot}{10}^{03}\text{\,}$ | $7\text{\cdot}{10}^{02}\text{\,}$ | $9\text{\cdot}{10}^{02}\text{\,}$ | $1\text{\cdot}{10}^{07}\text{\,}$ | $4\text{\cdot}{10}^{08}\text{\,}$ | $1\text{\cdot}{10}^{00}\text{\,}$
Li-Sb | $450\text{\,}$ | $2.4\text{\cdot}{10}^{-02}\text{\,}$ | $6.9\text{\cdot}{10}^{01}\text{\,}$ | $2.8\text{\cdot}{10}^{03}\text{\,}$ | $3\text{\cdot}{10}^{02}\text{\,}$ | $3\text{\cdot}{10}^{02}\text{\,}$ | $1\text{\cdot}{10}^{07}\text{\,}$ | $5\text{\cdot}{10}^{08}\text{\,}$ | $2\text{\cdot}{10}^{00}\text{\,}$
Li-Se | $375\text{\,}$ | $1.0\text{\cdot}{10}^{02}\text{\,}$ | $6.3\text{\cdot}{10}^{03}\text{\,}$ | $6.2\text{\cdot}{10}^{01}\text{\,}$ | $9\text{\cdot}{10}^{02}\text{\,}$ | $8\text{\cdot}{10}^{04}\text{\,}$ | $3\text{\cdot}{10}^{05}\text{\,}$ | $2\text{\cdot}{10}^{07}\text{\,}$ | $0.1\text{\,}$
Li-Sn | $400\text{\,}$ | $9.8\text{\cdot}{10}^{-03}\text{\,}$ | $8.4\text{\cdot}{10}^{01}\text{\,}$ | $8.5\text{\cdot}{10}^{03}\text{\,}$ | $3\text{\cdot}{10}^{02}\text{\,}$ | $2\text{\cdot}{10}^{02}\text{\,}$ | $4\text{\cdot}{10}^{07}\text{\,}$ | $2\text{\cdot}{10}^{09}\text{\,}$ | $4\text{\cdot}{10}^{00}\text{\,}$
Li-Te | $475\text{\,}$ | $1.7\text{\cdot}{10}^{-01}\text{\,}$ | $8.4\text{\cdot}{10}^{01}\text{\,}$ | $4.8\text{\cdot}{10}^{02}\text{\,}$ | $3\text{\cdot}{10}^{03}\text{\,}$ | $3\text{\cdot}{10}^{04}\text{\,}$ | $2\text{\cdot}{10}^{07}\text{\,}$ | $9\text{\cdot}{10}^{08}\text{\,}$ | $0.4\text{\,}$
Li-Zn | $486\text{\,}$ | $2.8\text{\cdot}{10}^{-02}\text{\,}$ | $6.3\text{\cdot}{10}^{02}\text{\,}$ | $2.2\text{\cdot}{10}^{04}\text{\,}$ | $2\text{\cdot}{10}^{02}\text{\,}$ | $8\text{\cdot}{10}^{01}\text{\,}$ | $5\text{\cdot}{10}^{07}\text{\,}$ | $5\text{\cdot}{10}^{09}\text{\,}$ | $10\text{\,}$
Mg-Sb | $700\text{\,}$ | $1.7\text{\cdot}{10}^{-02}\text{\,}$ | $3.0\text{\cdot}{10}^{01}\text{\,}$ | $1.8\text{\cdot}{10}^{03}\text{\,}$ | $5\text{\cdot}{10}^{02}\text{\,}$ | $6\text{\cdot}{10}^{02}\text{\,}$ | $3\text{\cdot}{10}^{06}\text{\,}$ | $4\text{\cdot}{10}^{08}\text{\,}$ | $20\text{\,}$
Na-Bi | $550\text{\,}$ | $1.0\text{\cdot}{10}^{-02}\text{\,}$ | $6.1\text{\cdot}{10}^{01}\text{\,}$ | $6.0\text{\cdot}{10}^{03}\text{\,}$ | $1\text{\cdot}{10}^{03}\text{\,}$ | $2\text{\cdot}{10}^{03}\text{\,}$ | $5\text{\cdot}{10}^{07}\text{\,}$ | $9\text{\cdot}{10}^{09}\text{\,}$ | $10\text{\,}$
Na-Hg | $275\text{\,}$ | $3.7\text{\cdot}{10}^{-03}\text{\,}$ | $6.4\text{\cdot}{10}^{00}\text{\,}$ | $1.7\text{\cdot}{10}^{03}\text{\,}$ | $1\text{\cdot}{10}^{04}\text{\,}$ | $2\text{\cdot}{10}^{04}\text{\,}$ | $1\text{\cdot}{10}^{08}\text{\,}$ | $1\text{\cdot}{10}^{10}\text{\,}$ | $10\text{\,}$
Na-Pb | $575\text{\,}$ | $1.2\text{\cdot}{10}^{-02}\text{\,}$ | $2.8\text{\cdot}{10}^{01}\text{\,}$ | $2.4\text{\cdot}{10}^{03}\text{\,}$ | $7\text{\cdot}{10}^{02}\text{\,}$ | $8\text{\cdot}{10}^{02}\text{\,}$ | $1\text{\cdot}{10}^{07}\text{\,}$ | $2\text{\cdot}{10}^{09}\text{\,}$ | $10\text{\,}$
Na-Sn | $625\text{\,}$ | $6.2\text{\cdot}{10}^{-03}\text{\,}$ | $2.1\text{\cdot}{10}^{01}\text{\,}$ | $3.4\text{\cdot}{10}^{03}\text{\,}$ | $3\text{\cdot}{10}^{02}\text{\,}$ | $2\text{\cdot}{10}^{02}\text{\,}$ | $1\text{\cdot}{10}^{07}\text{\,}$ | $9\text{\cdot}{10}^{08}\text{\,}$ | $10\text{\,}$

###  2.5 Onset of convection in the electrolyte and negative electrode

The onset of thermal convection in the negative electrode and electrolyte
depends on current density, material properties and layer thicknesses; it can
roughly be estimated using critical Rayleigh numbers from stability analysis
available in the literature. The Rayleigh number in the negative electrode can
be estimated in first approximation combining Eq. (9) and Eq. (14):

$Ra^{\mathrm{(\mathrm{N})}}\approx\frac{g\beta_{\mathrm{T}}^{\mathrm{(N)}}\rho_{\mathrm{el}}^{\mathrm{(E)}}j^{2}(\Delta
h^{\mathrm{(N)}})^{4}(\Delta
h^{\mathrm{(E)}})}{2k^{\mathrm{(N)}}\mathcal{D}_{\mathrm{T}}^{\mathrm{(N)}}\nu^{\mathrm{(N)}}}.$
(16)

In negative electrode, the temperature profile is linear for pure conduction,
see Figure 1. The top surface is rigid (wall) and at constant temperature, the
bottom surface can be well approximated by a uniform flux from the
electrolyte. This configuration with rigid bottom was studied by Sparrow et
al. [34]. They found a critical Rayleigh number of $Ra_{\mathrm{cr}}=1296$.
Köllner et al. [16] performed a linear stability analysis of the three-layer
system, and found for the negative electrode a critical Ra of 1290, when the
electrolyte is one tenth of the height of the electrodes.

The electrolyte has a parabolic temperature profile due to the Joule heating,
as shown in Figure 1. There, the Rayleigh number is well approximated assuming
an unstable length of half the layer thickness and the temperature difference
proposed in Eq. (15):

$Ra^{\mathrm{(\mathrm{E})}}\approx\frac{g\beta_{\mathrm{T}}^{\mathrm{(E)}}\rho_{\mathrm{el}}^{\mathrm{(E)}}j^{2}(\Delta
h^{\mathrm{(E)}})^{5}}{64k^{\mathrm{(E)}}\mathcal{D}_{\mathrm{T}}^{\mathrm{(E)}}\nu^{\mathrm{(E)}}}=Ra_{\mathrm{H}}.$
(17)

This formulation is free of parameters of the other layers and consistent with
previous work on internally heated convection [35]. Kulacki and Goldstein [36]
provided a complete description of the stability of internally heated layers
subject to different combinations of thermal and mechanical boundary
conditions. In our application, the electrolyte is mechanical and thermally
coupled with good heat conductors (liquid metal electrodes). The
conductivities and thicknesses of the electrodes, as well as any localised
cooling or heating play a role for the onset of thermal convection within the
molten salt. Due to the low and comparable thermal resistance of the two
layers, we assume an idealised condition in which the two interfaces are
assumed isothermal at the same temperature. The critical Rayleigh numbers in
this condition of rigid-rigid and free-free interfaces are respectively 583
and 266 [36]. In a similar arrangement, Sparrow et al. [34] show that the
critical Rayleigh number for rigid interfaces is bounded between 560 and 583
for sufficiently conductive boundaries.

In order to get an impression of the relevance of thermal convection in
realistic LMBs, we apply the critical Rayleigh numbers to the Li$||$Bi cell
studied here. The threshold current density, which corresponds to the critical
Rayleigh number can be computed as function of the layer thickness and
material properties using for the temperature difference an analytical model
of pure conduction (e.g. [17]), or by Eq. (14) and Eq. (15). The results, as
illustrated in Figure 2, show clearly that thermal convection in the negative
electrode is not guaranteed. It will appear only in thick layers. Using
typical moderate current densities in the order of
$0.25\text{\,}\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$ the electrolyte layer
needs to be fairly thick as well for convection to set in, see Figure 2b. An
efficient operation of the cell in terms of a low ohmic overpotential
($\eta_{\Omega}=\rho_{\mathrm{el}}j\Delta h^{\mathrm{(E)}}$) [1] requires thin
electrolytes, which does not guarantee flow will appear in the latter. The two
simple formulae (Eq. (14), Eq. (15)) provide a good approximation compared to
the more complex analytical model of pure conduction and can be used to have a
first prediction of temperature differences in the cells.

Figure 2: Critical current density and layer heights for onset of convection
in the negative electrode (a) and electrolyte (b). If not otherwise specified,
the layers are 10 mm thick; for the material properties, see Tab. 3.
Continuous lines are computed with the pure diffusion model of [17]. The dark
grey rhombi correspond to the parameters used for the numerical studies LARC.
The black hexagon in the left figure correspond to the HARC simulation.

In order to proceed in our investigation we simulate heat and mass transport
in LMB model.

##  3 Numerical model

The equations described in §2.1, §2.2 and §2.3 are discretized using the
finite volume method, in the framework provided by the C++ library OpenFOAM
[37, 38]. The parent-child mesh strategy is applied to model the three-layer
system with flat interfaces [39, 40]. This means that besides of the global
mesh taking care of the full domain, child meshes for each layer are
generated. Material properties can be seen as a field and they are
consistently assigned to each mesh point. A technique called mapping is used
to exchange field data across the different domains and meshes. For every time
step, the system of equations is solved consecutively, with the following
procedure: (a) The advection-diffusion equation for the temperature, Eq. (3),
is first solved on the parent mesh (full domain). Continuity of temperature
and heat flux at the interfaces are automatically enforced by this equation.
(b) In the second step, the mass transport equation, Eq. (1), is solved only
in the positive electrode (child mesh). The interface between the positive
electrode and the electrolyte is treated as a boundary for this equation. (c)
Following these two steps, the new density for the buoyancy force is computed
and mapped to each domain. (d) Finally, the Navier-Stokes equations are solved
sequentially in the negative electrode, the electrolyte and the positive
electrode. In order to ensure the mechanical coupling at the common
interfaces, appropriate velocity interface conditions are required, as
described in §3.1. Knowing that that the mechanical coupling is performed in
an explicit way (as boundary conditions), the Navier-Stokes equations need to
be solved several times in each phase until reaching convergence.

###  3.1 Discretisation of the mechanical layer coupling

In the OpenFOAM framework, the mechanical layer coupling must be ensured by ad
hoc interface conditions. The latter we base here on a simplified version of
Eq. 29 by Tukovic and Jasak [41], neglecting the effects of a moving mesh and
of surface tension. The exact derivation is the following: the interface
conditions for the velocity of a fluid 1 and 2 at a coupled interface read

$\displaystyle\boldsymbol{u}_{\mathrm{f1}}$
$\displaystyle=\boldsymbol{u}_{\mathrm{f2}},$ (18)
$\displaystyle\boldsymbol{u}\cdot\boldsymbol{n}$ $\displaystyle=0,$ (19)
$\displaystyle\mu_{1}\nabla\boldsymbol{u}_{\mathrm{t1}}\cdot\boldsymbol{n}_{1}$
$\displaystyle=\mu_{2}\nabla\boldsymbol{u}_{\mathrm{t2}}\cdot\boldsymbol{n}_{2},$
(20)

with $\boldsymbol{u}_{\mathrm{f}}$ denoting the velocity at the face between
fluid 1 and 2, $\boldsymbol{u}_{\mathrm{t}}$ the tangential velocity in the
cell centre, $\boldsymbol{u}_{\mathrm{tf}}$ the tangential velocity at the
face, $\boldsymbol{n}$ the face normal vector and $\mu$ the dynamic viscosity.
By combining the last two boundary conditions, we find

$\mu_{1}\frac{\boldsymbol{u}_{\mathrm{t1}}-\boldsymbol{u}_{\mathrm{tf}}}{\delta_{1}}=\mu_{2}\frac{\boldsymbol{u}_{\mathrm{tf}}-\boldsymbol{u}_{\mathrm{t2}}}{\delta_{2}},$
(21)

with $\delta$ denoting the distance between cell centre and face. This leads
to

$\boldsymbol{u}_{\mathrm{tf}}=w\cdot\boldsymbol{u}_{\mathrm{t1}}+(1-w)\cdot\boldsymbol{u}_{\mathrm{t2}}$
(22)

with the weighting factor

$w=\frac{\delta_{2}\cdot\mu_{1}}{\delta_{1}\mu_{2}+\delta_{2}\mu_{1}}.$ (23)

As a tangential velocity can be defined as

$\boldsymbol{u}_{\mathrm{t}}=\boldsymbol{u}\cdot(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n}),$
(24)

the boundary condition for the velocity at the interface might be written as

$\boldsymbol{u}_{\mathrm{tf}}=w\cdot\boldsymbol{u}_{1}\cdot(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n})+(1-w)\cdot\boldsymbol{u}_{2}\cdot(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n})=(\boldsymbol{I}-\boldsymbol{n}\boldsymbol{n})\cdot(w\cdot\boldsymbol{u}_{1}+(1-w)\cdot\boldsymbol{u}_{2}),$
(25)

in which $\boldsymbol{I}$ denotes the identity matrix.

##  4 Material properties and geometry

We focus in the following on the well-investigated Li$||$Bi LMB at the
reference temperature $T_{\mathrm{ref}}$ of
$450\text{\,}\mathrm{\SIUnitSymbolCelsius}$ in order to facilitate comparison
of our results with previous literature [17]. Considering the excellent
availability of material properties, we use a eutectic LiCl-KCl electrolyte,
although KCl is known to be unstable in contact with Li at higher temperature
[42]. Tab. 3 gives an overview on the material properties of the three
liquids. Most of the properties of the Li-Bi alloy are not exactly known.
Therefore, the alloy density and solutal expansion coefficient have been
estimated using Vegard’s law as suggested by Fazio et al. [31]. We assume a
reference molar fraction of $x_{\mathrm{Li}}=$ $20\text{\,}\mathrm{\char
37\relax}$, which corresponds to a Li-mass concentration of
$70.62\text{\,}\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$. While the specific heat
capacity is mass-weighted between Li and Bi according to the Neummann-Kopp’s
law [32], the kinematic viscosity is molar weighted [17]. Thermal and
electrical conductivity of the alloy are assumed to be equal to the one of
pure Bi at the reference temperature.

Table 3: Thermodynamic and transport properties for a Li$||$Bi LMB at $T_{\mathrm{ref}}=$450\text{\,}\mathrm{\SIUnitSymbolCelsius}$$. The alloy properties have been calculated for a Li-molar fraction of $20\text{\,}\mathrm{\char 37\relax}$. $c_{p}$ is the specific heat capacity, which is used to compute the thermal diffusivity as $\mathcal{D}_{\mathrm{T}}=k(\rho_{\mathrm{ref}}c_{p})^{-1}$. quantity | unit | N - Li | E - LiCl-KCl | P - Li-Bi | source
---|---|---|---|---|---
$\rho_{\mathrm{ref}}$ | $\mathrm{kg}\text{\,}{\mathrm{m}}^{-3}$ | $4.911\text{\cdot}{10}^{2}\text{\,}$ | $1.648\text{\cdot}{10}^{3}\text{\,}$ | $8.576\text{\cdot}{10}^{3}\text{\,}$ | [43, 44, 31]
$c_{p}$ | $\mathrm{J}\text{\,}{\mathrm{kg}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ | $4.237\text{\cdot}{10}^{3}\text{\,}$ | $1.330\text{\cdot}{10}^{3}\text{\,}$ | $1.7\text{\cdot}{10}^{2}\text{\,}$ | [43, 44, 31]
$k$ | $\mathrm{W}\text{\,}{\mathrm{m}}^{-1}\text{\,}{\mathrm{K}}^{-1}$ | $5.218\text{\cdot}{10}^{1}\text{\,}$ | $6.904\text{\cdot}{10}^{-1}\text{\,}$ | $1.421\text{\cdot}{10}^{1}\text{\,}$ | [43, 45, 31]
$\rho_{\mathrm{el}}$ | $\mathrm{\SIUnitSymbolOhm}\text{\,}{\mathrm{m}}^{-1}$ | $3.27\text{\cdot}{10}^{-7}\text{\,}$ | $6.358\text{\cdot}{10}^{-03}\text{\,}$ | $1.39\text{\cdot}{10}^{-6}\text{\,}$ | [43, 44, 31]
$\mathcal{D}_{\mathrm{T}}$ | ${\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$ | $2.51\text{\cdot}{10}^{-5}\text{\,}$ | $3.15\text{\cdot}{10}^{-7}\text{\,}$ | $9.75\text{\cdot}{10}^{-6}\text{\,}$ |
$Pr$ |  | $2.84\text{\cdot}{10}^{-2}\text{\,}$ | $6.29\text{\,}$ | $1.34\text{\cdot}{10}^{-2}\text{\,}$ |
$\nu$ | ${\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$ | $7.13\text{\cdot}{10}^{-7}\text{\,}$ | $1.983\text{\cdot}{10}^{-06}\text{\,}$ | $1.304\text{\cdot}{10}^{-7}\text{\,}$ | [43, 44, 31]
$\beta_{\mathrm{T}}$ | ${\mathrm{K}}^{-1}$ | $1.923\text{\cdot}{10}^{-4}\text{\,}$ | $3.2\text{\cdot}{10}^{-4}\text{\,}$ | $1.362\text{\cdot}{10}^{-4}\text{\,}$ | [43, 46, 31]
$\beta_{\rho_{\mathrm{A}}}$ | ${\mathrm{m}}^{3}\text{\,}{\mathrm{kg}}^{-1}$ |  |  | $2.122\text{\cdot}{10}^{-3}\text{\,}$ | [43, 31]
$\mathcal{D}_{\mathrm{AB}}$ | ${\mathrm{m}}^{2}\text{\,}{\mathrm{s}}^{-1}$ |  |  | $7\text{\cdot}{10}^{-9}\text{\,}$ | [3]
$\mathcal{M}_{\mathrm{A}}$ | $\mathrm{kg}\text{\,}{\mathrm{mol}}^{-1}$ | $6.94\text{\cdot}{10}^{-3}\text{\,}$ |  |  | [47]
$Sc$ |  |  |  | $18.6\text{\,}$ |

In order to resolve the flimsy plumes of solutal convection sufficiently – and
taking into account the performance limitations of the finite volume solver –
we confine our study to two-dimensional simulations. The geometry resembles
the one shown in Figure 3.
The thickness of the electrolyte is bounded by safety and efficiency reasons
(see [17] details); here we use 5 and $10\text{\,}\mathrm{mm}$. The necessary
negative and positive electrode layer heights are determined by the desired
composition of the positive electrode in the discharged state. Assuming that
all Bi of the positive electrode can be converted into the intermetallic phase
Li3Bi, a negative electrode height of about 20 mm would be sufficient. An
alternative approach is to discharge only in the liquid phase region of the
Li-Bi phase diagram ($x_{\text{Li}}\lesssim$$40\text{\,}\mathrm{\char
37\relax}$ at $450\text{\,}\mathrm{\SIUnitSymbolCelsius}$). This would result
in a minimum negative electrode height of approximately 4 mm required to fully
discharge a cell containing a 10 mm thick positive electrode. Here, we
investigate two configurations both with a $10\text{\,}\mathrm{mm}$ positive
electrode, referred to as HARC (high aspect ratio cell) and LARC (low aspect
ratio cell), respectively. The HARC has a $40\text{\,}\mathrm{mm}$ thick
negative electrode, see Figure 3. In this configuration, thermal convection is
also expected in the negative electrode at a current density of
$0.5\text{\,}\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$ (see the black hexagon in
Figure 2a). The low aspect ratio cell LARC with $10\text{\,}\mathrm{mm}$ thick
negative electrode is used to study the effect of solutal convection on
internally heated convection. Geometrical dimensions and applied current
densities of the different configurations studied are collected in Tab. 4.

Figure 3: Sketch of the 2D liquid metal battery model with temperature and concentration boundary conditions. The temperature is solved on the parent mesh (full domain), and the other equations on the child meshes (single regions). Table 4: Simulations list, geometry and current density. HT: heat transfer, MT: mass transfer. name | $L_{x}$ | $\Delta h^{\mathrm{(P)}}$ | $\Delta h^{\mathrm{(E)}}$ | $\Delta h^{\mathrm{(N)}}$ | $j$ | properties
---|---|---|---|---|---|---
| $\mathrm{mm}$ | $\mathrm{mm}$ | $\mathrm{mm}$ | $\mathrm{mm}$ | $\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$ |
HT_PS | 135 | 45 | 10 | 45 | 0.5 | [16]
MT_SEM | 40 | 10 | 10 | 10 | 0.5 | Tab. 3
HARC | 80 | 10 | 10 | 40 | 0.5, -0.5 | Tab. 3
LARC_H10 | 80 | 10 | 5 | 10 | 0.5, 0.75, 1.0, 1.5, 2.0 | Tab. 3
LARC_H5 | 80 | 10 | 10 | 10 | 0.5, 0.75, 1.0, 1.5, 2.0 | Tab. 3

##  5 Comparison with spectral solvers

In the absence of a suitable comparative case for the full problem, validation
of our model is performed separately for thermal and solutal convection, each
by comparing the OpenFOAM simulations with results of a pseudo-spectral and a
Fourier-spectral-element code respectively.

###  5.1 Thermal convection

The simulation of thermal convection in a three-layer LMB done with a pseudo-
spectral code by Köllner et al. was employed here as a reference. This test
case has already been used by Personnettaz at. [17]. The solver relies on a
spatial discretization in Fourier modes on the horizontal plane and Chebychev
polynomials in the vertical direction [16, 48]. The horizontal discretization
requires periodic boundary conditions at the lateral walls. We enforce the
same boundary conditions also in our OpenFOAM simulation to allow for
comparison.

The same cell as the one studied by Köllner et al. [16] is used here: a
Li$||$Pb-Bi LMB operating at $500\text{\,}\mathrm{\SIUnitSymbolCelsius}$; for
the material properties, see [16]. The geometry and current density are
collected in Tab. 4, with the simulation being denoted as HT_PS. The velocity
and temperature distributions, illustrated in Figure 4a-b, shows a clear
resemblance between both numerical codes. By averaging the velocity spatially
over the x-coordinate, and thereafter over a time of 50 s, we obtain the
vertical velocity profile as shown in Figure 4c. Very similar as in a recent
publication [17], we find that a mesh resolution of 300 control volumes over
the lateral size yields a converged velocity profile. Taking the average
velocity in steady state condition (over space, and 50 s again) we basically
obtain the same result – see Figure 4d.

Figure 4: Validation of thermal convection: (a) velocity and temperature
distribution of OpenFOAM, (b) the same for the pseudo-spectral code at t =
$998\text{\,}\mathrm{s}$. (c) Volume averaged velocity and (d) space and time-
averaged vertical velocity profile as a function of the lateral mesh
resolution. Simulation HT_PS.

###  5.2 Solutal convection

As a test case for solutal convection, we model the positive electrode of a
Li$||$Bi LMB, as illustrated in Figure 5. The dimensions are given in Tab. 4
(test case MT_SEM); the applied current density is
$0.5\text{\,}\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$. The comparison is
performed against a 2D simulation in Cartesian coordinates done with the
Fourier-spectral element code Semtex [49]. In this framework, a grid
independent solution was reached using 32 elements in vertical direction with
12x12 polynomials per element. In order to perform the comparison, the
boundary condition at the positive electrode/electrolyte interface is set to
no-slip in both solvers.

Figure 5: Comparison of solutal convection: (a) concentration distribution
obtained with Semtex and (b-d) OpenFOAM at different grid resolutions at
$t=$6.5\text{\,}\mathrm{s}$$, (e) time evolution of the mean mass
concentration of Li at the interface $\partial\Omega^{\mathrm{(P|E)}}$ and (f)
error of the Li interface concentration compared to Semtex, if averaged over
the last $15\text{\,}\mathrm{s}$.

The required grid resolution for solutal convection is substantially higher
than for thermal convection due to the lower diffusivity. Especially the plume
number and their width depend crucially on the mesh resolution, as can be
observed in Figure 5a-d. Fortunately, integral quantities, as the average mass
concentration at the electrolyte interface, are much less sensitive to the
grid resolution. This means, even if the plumes structure is not perfectly
resolved, the positive electrode is well mixed in any way. Therefore, the Li
interface-concentration results to be similar even for coarser meshes – as
illustrated in Figure 5e-f.

Considering that within each time step, the equation system needs to be solved
several times due to the explicit interface coupling, it is extremely
challenging to obtain the high accuracy of the spectral codes with OpenFOAM in
reasonable computational time. Furthermore, solutal convection in the regime
studied presents a strong dependence on initial condition and level of
numerical noise. In conclusion, we select a mesh able to capture the main flow
structures and as well as integral quantities. The cell is meshed with at
least 600 control volumes in the horizontal direction and 200 over the
positive electrode’s vertical extension. The boundaries, the positive
electrode area as well as the interfaces between the fluids are strongly
refined.

##  6 Results & discussion

We first present results of pure solutal convection in the three-layer system.
Then, we move to the interaction of heat and mass transfer.

###  6.1 Solutal convection

We first present solutal convection alone in the HARC cell – for the
dimensions, see Tab. 4. As reported before, the formation of the initial
plumes involves three distinct phases: first, we observe a quasi-stable state,
where a concentration boundary layer builds up. This layer is well described
by the pure diffusion solution, see Figure 6. Shortly after, the boundary
layer becomes unstable, and finally small plumes break rapidly out of the
layer as shown in Figure 6b [50]. This complete process takes only 4 s in our
conditions. As often observed before, the spatial wavelength increases then
with time due to merging of several plumes into larger convection cells [51,
52] – this happens in our case within only 3 seconds (Figure 6e). After less
than 10 s, the complete positive electrode is already well mixed. The
secondary vortices, which develop due to viscous coupling in the electrolyte
and negative electrode layer turn as expected in opposite direction (Figure
6a) [53]. The mean averaged flow velocity in the positive electrode reaches 5
$\mathrm{mm}\text{\,}{\mathrm{s}}^{-1}$, in the electrolyte 1
$\mathrm{mm}\text{\,}{\mathrm{s}}^{-1}$ and in the negative electrode only 0.1
$\mathrm{mm}\text{\,}{\mathrm{s}}^{-1}$ (Figure 6c). Figures 6f-h illustrate
the vertical component of the velocity over time along a horizontal line in
the middle of each layer. Obviously, the number of convection cells is very
similar within the three phases. Flow starts first in the positive electrode,
slightly later in the electrolyte and finally in the negative electrode.
Momentum is effectively transferred across the layers.

The most important conclusions of our example are that solutal convection
needs only few seconds to drive strong flow in the positive electrode and that
it mixes the electrolyte layer efficiently as well. There, the velocities
reach 20% of the mean averaged velocity of the positive electrode.

Figure 6: Flow structure and magnitude at 30 s (a), mass concentration
distribution at 5.5 s (b), maximum and volume averaged velocity (c), relative
velocity compared to the one in the positive electrode (d), mass concentration
at the positive electrode/electrolyte interface (e) and vertical velocity over
a horizontal line in the negative electrode (f), electrolyte (g) and positive
electrode (h). Simulation HARC.

###  6.2 Thermal convection and mass transport

The effect of thermal convection on mixing and cell efficiency has been
discussed controversially in the past. Kelley & Sadoway [13] as well as
Beltrán [54] suggested that bottom heating might lead to mm-scale velocities
in the positive electrode. Shen & Zikanov [15] found that even the flow in the
electrolyte might induce velocities in the order of 0.3
$\mathrm{mm}\text{\,}{\mathrm{s}}^{-1}$ in the positive electrode by viscous
coupling in small LMBs; they further predicted that convection will be much
faster in large cells. Moreover, it has been shown that electrochemical
heating might even lead to an unstable temperature profile in the positive
electrode [17] and that anti-convection might appear, as well [16]. However,
all these studies neglected the interaction with mass transport that is always
present in the positive electrode during operation.

We focus our first investigations on the HARC setup, with the dimensions given
in Tab. 4, starting with discharge and then studying charge. As already
discussed, during discharge, a stable density stratification forms in the
positive electrode due to the concentration, and to smaller extent due to the
temperature distribution. We study this phase using a current density of
$-0.5\text{\,}\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$. The coupling with the
electrolyte is the only source of motion in the positive electrode. As the
latter is too weak, the flow is strongly damped by the presence of the stable
stratification. In Figure 7a we observe that the velocity in the positive
electrode at discharge (red dots) is more than one order of magnitude less
than in the case of pure thermal convection (red line). The concentration
distribution in the positive electrode is unperturbed and follows exactly the
one predicted by an analytical solution, see Figure 7b [17]. The mean
velocities in the other layers are not substantially affected, as shown in
Figure 7a. Finally, also the temperature distribution in the cell is similar
to the one in the presence of pure thermal convection, Figure 7c-d, as well.

Figure 7: Volume averaged velocity (a), volume averaged temperature (c),
maximum temperature (d) in the three layers in presence of pure thermal
convection (HT), at discharge (HT + MT (d)) and charge (HT + MT (c)). Average
concentration difference between bottom wall and positive
electrode/electrolyte interface (b) during discharge, charge and charge with a
rigid interface. HARC setup.

When switching to charge, solutal flow occurs after only few seconds and
establishes vigorous motion in the positive electrode as shown in Figure 7a.
In previous investigations, the electrolyte/positive electrode interface had
been replaced by a no-slip boundary condition, when modelling the positive
electrode alone [3, 12]. Our more realistic multi-layer model shows that the
concentration difference within the positive electrode is now smaller compared
to the single region simulation, see Figure 7b. Obviously, replacing the
electrolyte/positive electrode interface by a no-slip boundary condition
attenuates flow and mixing too much, and overestimates therefore the mass
transport overpotential.

It is also immediately mixing the electrolyte layer, flow there appears ten
times earlier than in the scenario with pure thermal convection. The mean
temperature in the electrolyte is reduced by the presence of solutal
convection, as shown in Figure 7c.

Overall, solutal convection plays two roles: first it is efficiently
transporting and dissipating heat through the positive electrode, such as a
forced convection source at the positive electrode/electrolyte interface. And,
second, it is blending the electrolyte through the mechanical coupling. In the
HARC setup, the height of the negative electrode of $40\text{\,}\mathrm{mm}$
leads to a thermal Rayleigh number that exceeds the critical value for onset
of Rayleigh-Bénard convection even for the relatively low Joule heat produced
by a current density of $0.5\text{\,}\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$
(Figure 2, left). The electrolyte layer then essentially acts as a heater at
the lower interface of the negative electrode. Such a configuration has been
studied already Goluskin (case RB3, [55]). The typical regular convection
rolls of low-Ra Rayleigh-Bénard convection develop, as shown in Figure 8b. The
growth phase needs considerable time to develop, because first a sufficient
temperature gradient has to be established. However, solutal convection is
able to introduce disturbances in the top layer, which can trigger an earlier
onset.

Viscous coupling at the interfaces is important for the HARC geometry, as
well. Especially for the fully developed flow, the negative electrode takes
then the driving role and influences the flow in the electrolyte to a large
extent. This “dragging mode” [56] contributes to the formation of relatively
broad vortex cells in the electrolyte.

Computing the temporal mean value of the volume-averaged kinetic energy as
[15]

$\frac{E_{\text{kin}}}{\rho}=\frac{1}{V}\int_{V}|{\boldsymbol{u}}|^{2}dV$ (26)

we find that the lowest energy can be observed in the electrolyte for the HARC
case, see Figure 8. Since Joule heating is proportional to the square of the
current density, but mass transfer depends only linearly on $j$, we expect
that the maximum velocities in the negative electrode exceed those in the
positive electrode for the HARC geometry if the current density exceeds $0.5$
A/cm2.

###  6.3 Interaction of thermal and solutal convection in low aspect ratio
cells

We limit the discussion in this section to the charging phase in vertically
symmetric cells. The geometry is denoted as “low aspect ratio cell” (LARC) in
the text. Tab. 4 provides the dimensions and current density, which are also
illustrated by grey rhombi in the right diagram of Figure 2. The setup used
here is much closer to real cells, as already discussed in §3. Looking at
Figure 8a it can be seen that the temporal mean value of the volume-averaged
kinetic energy decreases monotonically from the positive electrode over the
electrolyte to the negative electrode.

Figure 8: Kinetic energy in the three layers of the cell for different current
densities and aspect ratios (a). Snapshots of two flow fields for the cell
with the high aspect ratio negative electrode (HARC, b) and for a LARC cell
(c) under charge with $0.5\text{\,}\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$.

Withing all simulations, the solutal convection’s flow velocity in the
positive electrode exceeds that in the other layers by at least an order of
magnitude in terms of kinetic energy.

In all cases, the layer coupling is dominated by viscous forces as already
found by other authors [15, 16] for thermal convection in liquid metal
batteries. The horizontal velocity components near the interfaces go into the
same direction on both sides of the interface, and interface normal velocity
components mirror each other.

The total heat released in the electrolyte grows with the electrolyte
thickness and the current density. While the negative electrode’s Rayleigh
number is smaller than the critical value for all current densities
considered, the electrolyte’s Rayleigh number exceeds the critical value only
for $j=2$ Acm-2 and a 5 mm thick electrolyte and for all current densities, if
the electrolyte thickness amounts to 10 mm.

Thermally driven convection in the cell means classical Rayleigh-Bénard
convection in the negative electrode, but internally heated convection in the
electrolyte. The former is well investigated, meticulously described in almost
all details and still subject of intense research; the latter received much
less attention but is nevertheless important for topics as diverse as mantle
convection, nuclear reactor engineering, and astrophysics. For a relatively
recent comprehensive review of both, see [55].

Solutal convection provides a momentum source to the electrolyte via viscous
coupling at the positive electrode/electrolyte interface.

Figure 9: Convective cooling parameter $T^{\star}/\overline{T}_{\text{max}}$
as a function of Ra${}_{\text{H}}$ for two different electrolyte heights. The
two insets show snapshots of the temperature and velocity distribution of the
fully developed flow in the electrolyte at the two parameter sets encircled in
grey. In horizontal direction only the central part of the electrolyte is
shown. LARC setup.

The influence of this momentum source on the flow in the electrolyte is
visible in the flow structure and the integral properties of the thermal
convection as can be seen in Figure 9. It shows the convective cooling
parameter $T^{\star}/\overline{T}_{\text{max}}$ as defined by Peckover and
Hutchinson [57] vs. the Rayleigh number based on the Joule heat release and
the electrolyte layer’s half height [35] described by Eq. (17). The convective
cooling parameter is the ratio of the maximum temperature $T^{\star}$ that
would occur if heat transport were purely conductive, to the maximum value of
the (height dependent) temperature profile $\overline{T}_{\text{max}}$ in the
fully developed flow averaged over time and in horizontal direction. The
location of $T^{\star}$ and $\overline{T}_{\text{max}}$ is marked in the
temperature profiles of Figure 10. Even for the smallest Ra${}_{\text{H}}$
($\text{Ra}_{\text{H}}=56$), the cooling parameter
($T^{\star}/\overline{T}_{\text{max}}=1.26$) exceeds one. This signals that
convection is present in the electrolyte intensifying heat transfer. The lower
inset in Figure 9
($\text{Ra}_{\text{H}}=226,T^{\star}/\overline{T}_{\text{max}}=1.37$) displays
snapshots of the instantaneous flow and temperature fields in the midsection
($-20\,\text{mm}<x<20\,\text{mm}$) of the electrolyte. Despite the Rayleigh
number being considerably lower than its critical value, four convection cells
are visible in the flow field. They are driven by viscous coupling at the
positive electrode/electrolyte interface ($z=$10\text{\,}\mathrm{mm}$$) where
the highest velocities can be observed. The temperature field is modulated by
the influence of convection and displays cooler regions where the flow
transports cold fluid from the boundaries into the bulk. Thus, the four
convection cells leave their footprints in the temperature field and alter it
considerably with respect to the purely conductive case. While the velocities
in the upper part of the electrolyte are smaller than that in the lower part,
the cells as well as their thermal footprint occupy the entire electrolyte
region.

Figure 10: Horizontally averaged vertical temperature profiles in the
electrolyte normalised with the maximum conductive temperature T⋆ for the 5 mm
and the 10 mm thick electrolyte and different current densities. LARC setup.

The upper inset in Figure 9
($\text{Ra}_{\text{H}}=28979,T^{\star}/\overline{T}_{\text{max}}=2.2$) shows
flow and temperature fields of well developed internally heated convection.
The flow field looks now quite different and so does the temperature field.
While the influence of viscous coupling is still visible at the positive
electrode/electrolyte interface, the highest velocities now occur in the
plumes descending from the electrolyte/negative electrode interface because of
the unstable temperature and thereby density distribution. As characteristic
for internally heated convection with heat transfer over both horizontal
boundaries, convection is very intense in the upper half of the layer, but
damped in the lower one. This is a consequence of the – on average – stable
stratification in the lower and the unstable stratification in the upper part
of the layer. This asymmetry results in penetrative convection (see, e.g.
[58]) where fluid undergoes convective motion in the upper part and penetrates
into the lower stably stratified layer. Such a situation leads to irregular
motion even for Rayleigh numbers close to the critical one as reported by a
number of authors, e.g., [59, 60, 61, 62] and as observed here as well.

Figure 10 displays the horizontally averaged vertical temperature profiles
normalised with the maximum conductive temperature for the two electrolyte
heights investigated. As already discussed in connection with Figure 9, we
observe convection for the smallest Rayleigh number investigated
(Ra${}_{\text{H}}=56$, corresponding to $j=0.5$
$\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$ in the 5 mm thick electrolyte). This
results in a moderate reduction of the maximum temperatures
$\overline{T}_{\text{max}}$ with respect to $T^{\star}$. However, the shape of
the temperature profiles is still relatively symmetric for the 5 mm thick
electrolyte with somewhat larger deviations for
$j=$$1.5\text{\,}\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$ and 2
$\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$. This changes for the 10 mm thick
electrolyte with its much higher Rayleigh numbers due to the stronger heat
release in the thicker layer. The reduction of $\overline{T}_{\text{max}}$
with respect to $T^{\star}$ is more pronounced as is the asymmetry of the
temperature profiles. The location of $\overline{T}_{\text{max}}$ shifts to
almost $(z-\Delta h^{\mathrm{(P)}})/\Delta h^{\mathrm{(E)}}=0.8$ for
$j=$$2.0\text{\,}\mathrm{A}\text{\,}{\mathrm{cm}}^{-2}$. The increase in
asymmetry with the heating rate in the electrolyte is again a feature that is
to be expected and known from the literature, e.g., [35, 57, 61, 62].

Figure 11: Nusselt numbers for the top and bottom interface of the electrolyte
with the electrodes (left) and Nusselt number ratio (right) vs. Rayleigh
number for different electrolyte thicknesses. LARC setup.

Corresponding to these changes in the mean temperature profiles, the heat
transfer conditions at the interfaces change. This is shown in Figure 11.
Figure 11a displays the Nusselt numbers at the electrolyte/negative electrode
interface (Nu${}_{\text{top}}$) and at the positive electrode/electrolyte
interface (Nu${}_{\text{bottom}}$). The Nusselt number definition is that used
by Kulacki and Goldstein [35]

$\text{Nu}=\frac{\Delta
h^{\mathrm{(E)}}|dT/dz|_{\text{interface}}}{T_{\text{max}}-T_{\text{interface}}},$
(27)

with the electrolyte height $\Delta h^{\mathrm{(E)}}$, the temperature
gradient at the interface $|dT/dz|_{\text{interface}}$, the maximum
temperature in the layer $T_{\text{max}}$, and the temperature of the
interface $T_{\text{interface}}$. Internally heated convection with equal
temperatures at the horizontal boundaries features typically higher Nusselt
numbers on the top of the layer and lower ones on the bottom, see, e.g. [35].
This is in contrast to what we observe for the 5 mm thick electrolyte layer.
Here, Nu${}_{\text{bottom}}$ exceeds Nu${}_{\text{top}}$ for all investigated
Rayleigh numbers. Since viscous coupling at the lower interface is the only
source of convective motion for Ra${}_{\text{H}}<\text{Ra}_{\text{H, crit}}$,
stronger convective heat transfer at the lower interface is an expected result
for these cases. But even for the higher Ra in the 5 mm thick electrolyte
convective heat transport due to viscous coupling outweighs that due to
thermal convection.

The conditions are reversed for the 10 mm thick electrolyte. As in purely
internally heated convection, Nu${}_{\text{top}}$ is larger than
Nu${}_{\text{bottom}}$ for all Ra${}_{\text{H}}$ investigated. Thermal
convection clearly dominates momentum transfer by viscous coupling. This
transition can be observed as well by looking at the flux ratio
Nu${}_{\text{top}}$/(Nu${}_{\text{top}}$+Nu${}_{\text{bottom}}$) in Figure
11b. The flux ratio would be 0.5 for the conductive case but lies below this
value for the 5 mm thick electrolyte, i.e., more heat is transferred through
the positive electrode/electrolyte interface than through the
electrolyte/negative electrode interface. For the 10 mm thick electrolyte, the
flux ratio is larger than 0.5 because more heat leaves the electrolyte through
the top interface than through the bottom one.

##  7 Summary

As liquid metal batteries (LMBs) operate as concentration cells, mixing is
extremely important. After discussing that concentration overpotentials might
appear in the positive electrode, we briefly gave an overview on the various
flow phenomena being present in LMBs. We showed that the effects of solutal
convection on the cell voltage have already been observed 60 years ago.
However, the source of these effects was explained only very recently by
Personnettaz et al. [3].

With this background and motivation, we studied the interaction of solutal and
thermal convection in a Li$||$Bi LMB. A brief theoretical analysis suggests
that solutal convection will always be able to mix the stable thermal
stratification. On the other hand, the thermal flow is not expected to be
strong enough to mix the compositionally stably stratified positive electrode.
Using critical Rayleigh numbers from literature, we predicted further that
thermal convection will appear only in very thick electrolyte layers and
seldomly in the negative electrode.

In order to study the interplay of solutal and thermal convection in more
detail, we developed a numerical model with fixed interfaces, solving for the
temperature and flow field within the complete cell, and the concentration
distribution in the positive electrode. The validation, which was performed by
comparison with two (pseudo-)spectral codes, suggests that the mesh resolution
for solutal convection needs to be considerably greater than for thermal flow.

In a first step, we simulated solutal convection, which is present only in the
positive electrode during the charge phase. We found that it almost instantly
produces a very strong mixing effect in the cell. The velocities, which are
induced in the electrolyte via viscous coupling, reach still 20% of the mean
average velocity compared to the positive electrode. Therefore, we deduced
that compositional flow mixes both the positive electrode and electrolyte. We
further found that it is important to model the mechanically coupled
electrolyte/positive electrode interface. Replacing the interface with a
simplified no-slip boundary condition attenuates flow and mixing and will
therefore overestimate the mass transport overvoltage [10, 3].

In a second step, we simulated thermal convection in discharge. As it appears
first in the electrolyte, it is able to mix this region well. However, the
viscous coupling is not strong enough to induce any relevant flow in the
positive electrode.

Finally, we studied thermal and solutal convection together. We showed that
thermal convection needs roughly ten times longer to develop, compared to the
solutal flow. The layer coupling is always viscous. While thermal convection
does not develop in thin electrolyte layers, except for the highest
investigated current density, thick electrolyte layers are always dominated by
thermally driven flow. Accordingly, thin electrolyte layers feature relatively
regular convection cells during charge. These extend over the full height and
are driven by viscous coupling from the positive electrode that is well mixed
by solutal convection. In contrast, the Rayleigh numbers of thick electrolyte
layers are well above the critical ones for all cases considered. Averaged
temperature profiles show a pronounced asymmetry typical for penetrative
convection. Nusselt numbers at the upper boundary of the electrolyte layer
exceed that at the lower boundary for thick electrolytes. The opposite holds
true for thin electrolytes, where viscous coupling drives the most intense
flow of the lower boundary and the cooling parameter is smaller than 0.5.

In summary, we would like to stress that solutal convection is a very powerful
effect. Although it appears only in the positive electrode at charge, it is
able to mix then the electrolyte via viscous coupling at the interface. In
contrast, thermal convection can appear only in comparably thick electrolyte
layers, but will then allow mixing them during charge and discharge.

## Acknowledgements

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under grant agreement No 963599, from the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) by award
number 338560565 and in frame of the Helmholtz - RSF Joint Research Group
“Magnetohydrodynamic instabilities: Crucial relevance for large scale liquid
metal batteries and the sun-climate connection”, contract No. HRSF-0044 and
RSF-18-41-06201. The authors are grateful to T. Köllner for the pseudo-
spectral simulation of thermal convection in a three-layer LMB, used here for
validation. Fruitful discussions with S. Bénard, W. Herreman and C. Nore are
gratefully acknowledged.

## References

  * [1] H. Kim, D. A. Boysen, T. Ouchi, D. R. Sadoway, Calcium - bismuth electrodes for large - scale energy storage (liquid metal batteries), J. Power Sources 241 (2013) 239–248.
  * [2] D. H. Kelley, T. Weier, Fluid mechanics of liquid metal batteries, Appl. Mech. Rev. 70 (2) (2018) 020801.  doi:10.1115/1.4038699.
  * [3] P. Personnettaz, S. Landgraf, M. Nimtz, N. Weber, T. Weier, Mass transport induced asymmetry in charge/discharge behavior of liquid metal batteries, Electrochem. Commun. 105 (2019) 106496.  doi:10.1016/j.elecom.2019.106496.
  * [4] B. Agruss, H. R. Karas, V. L. Decker, Design and development of a liquid metal fuel cell, Tech. Rep. ASD-TDR-62-1045, General Motors Corporation (1962).
  * [5] E. J. Cairns, C. E. Crouthamel, A. K. Fischer, M. S. Foster, J. C. Hesson, C. E. Johnson, H. Shimotake, A. D. Tevebaugh, Galvanic Cells with Fused-Salt Electrolytes, ANL-7316, Argonne National Laboratory, 1967\.
  * [6] R. C. Vogel, M. Levenson, E. R. Proud, J. Royal, Chemical engineering division research highlights, Tech. Rep. ANL-7350, Argonne National Laboratory (1967).
  * [7] M. S. Foster, Laboratory Studies of Intermetallic Cells, in: C. E. Crouthamel, H. L. Recht (Eds.), Regenerative EMF Cells, American Chemical Society, 1967, pp. 136–148.
  * [8] R. Ashour, D. H. Kelley, A. Salas, M. Starace, N. Weber, T. Weier, Competing forces in liquid metal electrodes and batteries, J. Power Sources 378 (2018) 301–310.  doi:10.1016/j.jpowsour.2017.12.042.
  * [9] N. Weber, M. Nimtz, P. Personnettaz, A. Salas, T. Weier, Electromagnetically driven convection suitable for mass transfer enhancement in liquid metal batteries, Appl. Therm. Eng. 143 (2018) 293–301.  doi:10.1016/j.applthermaleng.2018.07.067.
  * [10] W. Herreman, S. Bénard, C. Nore, P. Personnettaz, L. Cappanera, J.-L. Guermond, Solutal buoyancy and electrovortex flow in liquid metal batteries, Phys. Rev. Fluids 5 (7) (2020) 074501.  doi:10.1103/PhysRevFluids.5.074501.
  * [11] N. Weber, M. Nimtz, P. Personnettaz, T. Weier, D. Sadoway, Numerical simulation of mass transfer enhancement in liquid metal batteries by means of electro-vortex flow, J. Power Sources Adv. 1 (2020) 100004.  doi:10.1016/j.powera.2020.100004.
  * [12] W. Herreman, C. Nore, L. Cappanera, J.-L. Guermond, Efficient mixing by swirling electrovortex flows in liquid metal batteries, J. Fluid Mech. 915 (2021) A17.  doi:10.1017/jfm.2021.79.
  * [13] D. H. Kelley, D. R. Sadoway, Mixing in a liquid metal electrode, Phys. Fluids 26 (5) (2014) 057102.  doi:10.1063/1.4875815.
  * [14] A. Beltrán, MHD natural convection flow in a liquid metal electrode, Appl. Therm. Eng. 114 (2016) 1203–1212.  doi:10.1016/j.applthermaleng.2016.09.006.
  * [15] Y. Shen, O. Zikanov, Thermal convection in a liquid metal battery, Theor. Comput. Fluid Dyn. 30 (4) (2016) 275–294.  doi:10.1007/s00162-015-0378-1.
  * [16] T. Köllner, T. Boeck, J. Schumacher, Thermal Rayleigh-Marangoni convection in a three-layer liquid-metal-battery model, Phys. Rev. E 95 (2017) 053114.  doi:10.1103/PhysRevE.95.053114.
  * [17] P. Personnettaz, P. Beckstein, S. Landgraf, T. Köllner, M. Nimtz, N. Weber, T. Weier, Thermally driven convection in Li$——$Bi liquid metal batteries, J. Power Sources 401 (2018) 362–374.  doi:10.1016/j.jpowsour.2018.08.069.
  * [18] F. Stefani, V. Galindo, C. Kasprzyk, S. Landgraf, M. Seilmayer, M. Starace, N. Weber, T. Weier, Magnetohydrodynamic effects in liquid metal batteries, IOP Conf. Ser. Mater. Sci. Eng. 143 (2016) 012024.  doi:10.1088/1757-899X/143/1/012024.
  * [19] H. Kim, D. A. Boysen, J. M. Newhouse, B. L. Spatocco, B. Chung, P. J. Burke, D. J. Bradwell, K. Jiang, A. A. Tomaszowska, K. Wang, W. Wei, L. A. Ortiz, S. A. Barriga, S. M. Poizeau, D. R. Sadoway, Liquid Metal Batteries: Past, Present, and Future, Chem. Rev. 113 (3) (2013) 2075–2099.  doi:10.1021/cr300205k.
  * [20] S. Lawroski, R. C. Vogel, V. H. Munnecke, Chemical engineering division summary report, Tech. Rep. ANL-6543, Argonne National Laboratory (1962).
  * [21] C. E. Crouthamel, H. L. Recht (Eds.), Regenerative EMF Cells, Vol. 64, American Chemical Society, 1967.
  * [22] H. Shimotake, E. J. Cairns, Bimetallic galvanic cells with fused-salt electrolytes, in: Advances in Energy Conversion Engineering, ASME, Florida, 1967, pp. 951–962.
  * [23] D. J. Bradwell, Liquid Metal Batteries: Ambipolar Electrolysis and Alkaline Earth Electroalloying Cells, Ph.D. thesis, Massachusetts Institute of Technology (2011).
  * [24] T. Ouchi, H. Kim, X. Ning, D. R. Sadoway, Calcium-Antimony Alloys as Electrodes for Liquid Metal Batteries, J. Electrochem. Soc. 161 (12) (2014) A1898–A1904.
  * [25] T. Ouchi, H. Kim, B. L. Spatocco, D. R. Sadoway, Calcium-based multi-element chemistry for grid-scale electrochemical energy storage, Nat. Commun. 7 (2016) 10999.  doi:10.1038/ncomms10999.
  * [26] A. Blanchard, Enabling multi-cation electrolyte usage in LMBs for lower cost and operating temperature, Ph.D. thesis, Massachusetts Institute of Technology (2013).
  * [27] B. Agruss, The Thermally Regenerative Liquid-Metal Cell, J. Electrochem. Soc. 110 (11) (1963) 1097–1103.
  * [28] R. B. Bird, W. E. Stewart, E. N. Lightfoot, Transport Phenomena, Wiley, 2002.  doi:10.1115/1.1424298.
  * [29] N. Weber, C. Duczek, G. M. Horstmann, S. Landgraf, M. Nimtz, P. Personnettaz, T. Weier, D. R. Sadoway, Cell voltage model for Li-Bi liquid metal batteries, Appl. Energy, under revision (2021).
  * [30] L. G. Leal, Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes, Cambridge Series in Chemical Engineering, Cambridge University Press, 2007.  doi:10.1017/CBO9780511800245.
  * [31] C. Fazio, Handbook on Lead-bismuth Eutectic Alloy and Lead Properties, Materials Compatibility, Thermal-hydraulics and Technologies, Tech. Rep. 7268, Nuclear Energy Agency (2015).
  * [32] T. Iida, R. I. L. Guthrie, The Thermophysical Properties of Metallic Liquids, Oxford University Press, United Kingdom, 2015.
  * [33] A. K. Roy, R. P. Chhabra, Prediction of Solute Diffusion Coefficients in Liquid Metals, Metall. Trans. A 19 (1988) 273–279.
  * [34] E. M. Sparrow, R. J. Goldstein, V. K. Jonsson, Thermal instability in a horizontal fluid layer: Effect of boundary conditions and non-linear temperature profile, J. Fluid Mech. 18 (4) (1964) 513–528.  doi:10.1017/S0022112064000386.
  * [35] F. A. Kulacki, R. J. Goldstein, Thermal convection in a horizontal fluid layer with uniform volumetric energy sources, J. Fluid Mech. 55 (02) (1972) 271–287.
  * [36] F. A. Kulacki, R. J. Goldstein, Hydrodynamic instability in fluid layers with uniform volumetric energy sources, Appl. Sci. Res. 31 (2) (1975) 81–109.  doi:10.1007/bf01795829.
  * [37] H. Jasak, T. Uroić, Practical Computational Fluid Dynamics with the Finite Volume Method, Springer International Publishing, 2020, pp. 103–161.  doi:10.1007/978-3-030-37518-8_4.
  * [38] H. G. Weller, G. Tabor, H. Jasak, C. Fureby, A tensorial approach to computational continuum mechanics using object-oriented techniques, Comput. Phys. 12 (6) (1998) 620–631.
  * [39] S. B. Beale, H.-W. Choi, J. G. Pharoah, H. K. Roth, H. Jasak, D. H. Jeon, Open-source computational model of a solid oxide fuel cell, Comput. Phys. Commun. 200 (2016) 15–26.  doi:10.1016/j.cpc.2015.10.007.
  * [40] N. Weber, P. Beckstein, V. Galindo, M. Starace, T. Weier, Electro-vortex flow simulation using coupled meshes, Comput. Fluids 168 (2018) 101–109.  doi:10.1016/j.compfluid.2018.03.047.
  * [41] Ž. Tuković, H. Jasak, A moving mesh finite volume interface tracking method for surface tension dominated interfacial fluid flow, Comput. Fluids 55 (2012) 70–84.  doi:10.1016/j.compfluid.2011.11.003.
  * [42] M. S. Foster, S. E. Wood, C. E. Crouthamel, Thermodynamics of Binary Alloys. I. The Lithium-Bismuth System, Inorg. Chem. 3 (10) (1964) 1428–1431.
  * [43] R. W. Ohse (Ed.), Handbook of Thermodynamic and Transport Properties of Alkali Metals, no. 30 in IUPAC Chemical Data Series, Blackwell Scientific, Oxford, 1985.
  * [44] C. J. Raseman, H. Susskind, G. Farber, W. McNulty, F. Salzano, Engineering experience at Brookhaven National Laboratory in handling fused chloride salts, Tech. Rep. BNL 627, Brookhaven National Laboratory (1960).
  * [45] G. J. Janz, R. P. T. Tomkins, C. B. Allen, Molten Salts: Volume 4, Part 4 Mixed Halide Melts Electrical Conductance, Density, Viscosity, and Surface Tension Data, J. Phys. Chem. Ref. Data 8 (1) (1979) 125–302.  doi:10.1063/1.555590.
  * [46] G. J. Janz, Thermodynamic and transport properties for molton salts: Correlation equations for critically evaluated density, surface tension, electrical conductance, and viscosity data, J. Phys. Chem. Ref. Data 17 (1988).
  * [47] P. J. Linstrom, W. G. Mallard, NIST Chemistry WebBook, NIST Standard Reference Database Number 69, National Institute of Standards and Technology, Gaithersburg MD, 20899, 2018.  doi:10.18434/T4D30.
  * [48] T. Köllner, M. Rossi, F. Broer, T. Boeck, Chemical convection in the methylene-blue–glucose system: Optimal perturbations and three-dimensional simulations, Phys. Rev. E 90 (5) (2014) 053004.  doi:10.1103/PhysRevE.90.053004.
  * [49] H. Blackburn, D. Lee, T. Albrecht, J. Singh, Semtex: A spectral element–Fourier solver for the incompressible Navier–Stokes equations in cylindrical or Cartesian coordinates, Comput. Phys. Commun. 245 (2019) 106804.  doi:10.1016/j.cpc.2019.05.015.
  * [50] R. J. Goldstein, R. J. Volino, Onset and development of natural convection above a suddenly heated horizontal surface, J. Heat Transf. 117 (4) (1995) 808–821.  doi:10.1115/1.2836296.
  * [51] J. Philippi, M. Berhanu, J. Derr, S. Courrech du Pont, Solutal convection induced by dissolution, Phys. Rev. Fluids 4 (10) (2019) 103801.  doi:10.1103/PhysRevFluids.4.103801.
  * [52] S. Backhaus, K. Turitsyn, R. E. Ecke, Convective Instability and Mass Transport of Diffusion Layers in a Hele-Shaw Geometry, Phys. Rev. Lett. 106 (10) (2011) 104501.  doi:10.1103/PhysRevLett.106.104501.
  * [53] A. Prakash, K. Yasuda, F. Otsubo, K. Kuwahara, T. Doi, Flow coupling mechanisms in two-layer Rayleigh-Benard convection, Experiments in Fluids 23 (3) (1997) 252–261.  doi:10.1007/s003480050108.
  * [54] Springer Handbook of Electrochemical Energy, 1st Edition, Springer Berlin Heidelberg, New York, NY, 2016.
  * [55] D. Goluskin, Internally Heated Convection and Rayleigh-Bénard Convection, Springer International Publishing, Cham, 2016.  doi:10.1007/978-3-319-23941-5.
  * [56] D. Johnson, R. Narayanan, Geometric effects on convective coupling and interfacial structures in bilayer convection, Phys. Rev. E 56 (5) (1997) 5462–5472.
  * [57] R. S. Peckover, I. H. Hutchinson, Convective rolls driven by internal heat sources, Phys. Fluids 17 (7) (1974) 1369–1371.
  * [58] G. Veronis, Penetrative convection, Astrophys. J. 137 (1963) 641–663.
  * [59] M. Jahn, H. H. Reineke, Free convection heat transfer with internal heat sources calculations and measurements, in: Proc. International Heat Transfer Conference 5, Tokyo, 1974, pp. 74–78.
  * [60] J. C. Ralph, R. McGreevy, R. S. Peckover, Experiments in turbulent thermal convection driven by internal heat sources, in: Heat Transfer and Turbulent Buoyant Convection, Vol. 2, McGraw-Hill, 1976, pp. 587–599.
  * [61] M. Tveitereid, Thermal convection in a horizontal fluid layer with internal heat sources, Int. J. Heat Mass Transf. 21 (1978) 335–339.
  * [62] F. A. Kulacki, D. E. Richards, Natural convection in plane layers and cavities with volumetric energy sources, in: S. Kakac, W. Aung, R. Viskanta (Eds.), Natural Convection: Fundamentals and Applications, Hemisphere, Washington, 1985, pp. 179–255.