Thermal Rayleigh-Marangoni convection in a three-layer liquid-metal-battery model

TL;DR

This study investigates the combined effects of Rayleigh-Bénard and Marangoni convection in a three-layer liquid metal battery model, revealing how these instabilities interact and influence heat transfer and flow patterns.

Contribution

It provides a linear stability analysis and direct numerical simulations of convection modes in a liquid metal battery model, highlighting the role of Marangoni effects and layer thickness.

Findings

Marangoni effects increase instability growth rates.

Critical Grashof and Marangoni numbers decrease with thicker middle layer.

Convection thresholds are realistic for laboratory-sized batteries.

Abstract

The combined effects of buoyancy-driven Rayleigh-B\'{e}nard convection (RC) and surface tension-driven Marangoni convection (MC) are studied in a triple-layer configuration which serves as a simplified model for a liquid metal battery (LMB). The three-layer model consists of a liquid metal alloy cathode, a molten salt separation layer, and a liquid metal anode at the top. Convection is triggered by the temperature gradient between the hot electrolyte and the colder electrodes, which is a consequence of the release of resistive heat during operation. We present a linear stability analysis of the state of pure thermal conduction in combination with three-dimensional direct numerical simulations of the nonlinear turbulent evolution on the basis of a pseudospectral method. Five different modes of convection are identified in the configuration, which are partly coupled to each other: RC in the upper electrode, RC with internal heating in the molten salt layer, MC at both interfaces between molten salt and electrode as well as anti-convectionin the middle layer and lower electrode. The linear stability analysis confirms that the additional Marangoni effect in the present setup increases the growth rates of the linearly unstable modes, i.e. Marangoni and Rayleigh-B\'{e}nard instability act together in the molten salt layer. The critical Grashof and Marangoni numbers decrease with increasing middle layer thickness. The calculated thresholds for the onset of convection are found for realistic current densities of laboratory-sized LMBs. The global turbulent heat transfer follows scaling predictions for internally heated RC. The global turbulent momentum transfer is comparable with turbulent convection in the classical Rayleigh-B\'{e}nard case.

Tables

Table 1. Table 1: List of all material parameters of a three-layer battery system Li || LiCl-KCl || Pb-Bi at an operating temperature of 500 ∘ C. All references for the numerical values are given in the text.

Physical quantity	Symbol	SI-Unit	Value
mass density	${\rho }_{ref}^{\left(1\right)}$	kg/m3	$1.0065\times {10}^4$
	${\rho }_{ref}^{\left(2\right)}$	kg/m3	$1597.9$
	${\rho }_{ref}^{\left(3\right)}$	kg/m3	$484.7$
kinematic viscosity	${\nu }^{\left(1\right)}$	m${}_{}{}^2/$s	$1.29\times {10}^{-7}$
	${\nu }^{\left(2\right)}$	m2/s	$1.38\times {10}^{-6}$
	${\nu }^{\left(3\right)}$	m2/s	$6.64\times {10}^{-7}$
thermal diffusivity	${\kappa }^{\left(1\right)}$	m2/s	$1.015\times {10}^{-5}$
	${\kappa }^{\left(2\right)}$	m2/s	$1.90\times {10}^{-7}$
	${\kappa }^{\left(3\right)}$	m2/s	$2.48\times {10}^{-5}$
thermal conductivity	${\lambda }^{\left(1\right)}$	W/(mK)	14.41
	${\lambda }^{\left(2\right)}$	W/(mK)	0.365
	${\lambda }^{\left(3\right)}$	W/(mK)	50.12
specific heat	$C_p^{\left(1\right)}$	J/(kgK)	141.05
	$C_p^{\left(2\right)}$	J/(kgK)	1201.6
	$C_p^{\left(3\right)}$	J/(kgK)	4169
electrical conductivity	${\sigma }_e^{\left(1\right)}$	S/m	$7.81\times {10}^5$
	${\sigma }_e^{\left(2\right)}$	S/m	$187.1$
	${\sigma }_e^{\left(3\right)}$	S/m	$3\times {10}^6$
expansion coefficient	${\beta }_T^{\left(1\right)}$	1/K	$-1.28\times {10}^{-4}$
	${\beta }_T^{\left(2\right)}$	1/K	$-3.32\times {10}^{-4}$
	${\beta }_T^{\left(3\right)}$	1/K	$-2.08\times {10}^{-4}$
change in interfacial	${\sigma }_{ref}^I{\alpha }_T^I$	N/(mK)	$3.1\times {10}^{-5}$
tension	${\sigma }_{ref}^{II}{\alpha }_T^{II}$	N/(mK)	$7.98\times {10}^{-5}$

Table 2. Table 2: Derived non-dimensional parameters and parameter ratios which are used for the DNS of the nonlinear evolution. All ratios follow from ( 4 ). Parameter ξ 𝜉 \xi is given by ( 31 ). Parameters that depend on layer heights and electrical current density are computed with d ( 1 ) = d ( 2 ) = d ( 3 ) superscript 𝑑 1 superscript 𝑑 2 superscript 𝑑 3 d^{(1)}=d^{(2)}=d^{(3)} =20 mm and j 0 = 3 subscript 𝑗 0 3 j_{0}=3 kA/m 2 .

Quantity	Symbol	Value
Grashof number	$G$	-3.97$\times {10}^6$
Marangoni number	Ma	-310.02
Prandtl number	Pr(1)	0.0127
	$P{r}^{\left(2\right)}$	7.24
	$P{r}^{\left(3\right)}$	0.0268
Layer height	$d_{21}$	1
	$d_{31}$	1
Thermal conductivity	${\lambda }_{21}$	0.0253
	${\lambda }_{31}$	3.48
Mass density	${\rho }_{21}$	0.159
	${\rho }_{31}$	0.048
Kinematic viscosity	${\nu }_{21}$	10.66
	${\nu }_{31}$	5.14
Thermal diffusivity	${\kappa }_{21}$	0.019
	${\kappa }_{31}$	2.44
Thermal expansion	${\beta }_{21}$	2.59
	${\beta }_{31}$	1.63
Interfacial tension	$\xi$	2.57
Electrical conductivity	${\sigma }_{e,21}$	2.396$\times {10}^{-4}$
	${\sigma }_{e,31}$	3.841
Velocity unit	$U_{vis}$	6.46$\times {10}^{-6}$m/s
Time unit	${\tau }_{vis}$	3097 s
Length unit	$L_{vis}$	20 mm
Temperature unit	$\Theta$	6.59K

Table 3. Table 3: List of full nonlinear DNS runs to verify the eigenvalues and eigenvectors of the numerical solution of the linear stability analysis. Grashof and Marangoni numbers as well as spectral resolutions are listed. All other parameters not listed here correspond with those listed in Tab. 2 .

Case	$d_{21}$	$G$	$\xi$	$Ma$	$l_x$	$N_x$	$N_y$	$N_z^{\left(1\right)}$	$N_z^{\left(2\right)}$	$N_z^{\left(3\right)}$
LS01	1	0	2.57	-76	12	128	128	32	64	32
LS02	1	0	2.57	-77	12	128	128	32	64	32
LS03	1	-1.29$\times {10}^4$	0	0	12	128	128	32	64	32
LS04	1	-1.30$\times {10}^4$	0	0	12	128	128	32	64	32
LS05	0.3	0	2.57	-229	12	128	128	32	64	32
LS06	0.3	0	2.57	-230	12	128	128	32	64	32
LS07	0.3	-3.75$\times {10}^5$	0	0	6	128	128	32	64	32
LS08	0.3	-3.85$\times {10}^5$	0	0	6	128	128	32	64	32
LS09	0.1	-2.30$\times {10}^6$	0	0	6	128	128	32	64	32
LS10	0.1	-2.40$\times {10}^6$	0	0	6	128	128	32	64	32
LS11	1	-1.3314$\times {10}^4$	2.57	-10.16	6	128	128	32	64	32
LS12	1	-1.3314$\times {10}^4$	0	0	6	128	128	32	64	32

Table 4. Table 4: List of full nonlinear simulation runs with different electrical current densities j 0 subscript 𝑗 0 j_{0} and/or layer heights d ( 1 ) superscript 𝑑 1 d^{(1)} . For each case displayed, an additional simulation with the Marangoni effect disabled is performed. Again N x = N y subscript 𝑁 𝑥 subscript 𝑁 𝑦 N_{x}=N_{y} and l x = l y subscript 𝑙 𝑥 subscript 𝑙 𝑦 l_{x}=l_{y} . The runs with superscript “+ ” are found to be linearly stable.

Case	$d_{21}$	$G$	${\tilde{G}}_c$	$k_c^{\mathrm{RC}}$	$\mathrm{Ma}$	$\stackrel{~}{}M{a}_c$	$k_c^{\mathrm{MC}}$	$N_x\left(=N_y\right)$	$l_x\left(=l_y\right)$	$N_z^{\left(1\right)}$	$N_z^{\left(2\right)}$	$N_z^{\left(3\right)}$	$d^{\left(1\right)}$ [mm]	$j_0$ [kA/m2]
H1+	1	-3.88$\times {10}^3$	-1.29$\times {10}^4$	3	-4.84	-76.4	2.10	128	6	32	64	32	5	3
H2	1	-1.24$\times {10}^5$	-1.29$\times {10}^4$	3	-38.75	-76.4	2.10	128	6	32	64	32	10	3
H3	1	-9.42$\times {10}^5$	-1.29$\times {10}^4$	3	-130.79	-76.4	2.10	128	6	32	64	32	15	3
H4	1	-3.97$\times {10}^6$	-1.29$\times {10}^4$	3	-310.02	-76.4	2.10	256	6	32	64	32	20	3
H5	1	-3.01$\times {10}^7$	-1.29$\times {10}^4$	3	-1.05e3	-76.4	2.10	256	6	64	128	64	30	3
H6	1	-1.27$\times {10}^8$	-1.29$\times {10}^4$	3	-2.48e3	-76.4	2.10	512	6	64	128	64	40	3
JM1	0.5	-1.10$\times {10}^5$	-9.23$\times {10}^4$	5.75	-8.61	-138.8	3.75	128	3	32	64	32	20	1
JM2	0.5	-9.92$\times {10}^5$	-9.23$\times {10}^4$	5.75	-77.50	-138.8	3.75	128	3	32	64	32	20	3
JM3	0.5	-2.76$\times {10}^6$	-9.23$\times {10}^4$	5.75	-215.29	-138.8	3.75	256	3	64	128	64	20	5
JM4	0.5	-1.10$\times {10}^7$	-9.23$\times {10}^4$	5.75	-861.16	-138.8	3.75	256	3	64	128	64	20	10
JS1+	0.3	-3.97$\times {10}^4$	-3.81$\times {10}^5$	9.30	-3.1	-229.2	6.25	128	2	32	32	32	20	1
JS2+	0.3	-3.57$\times {10}^5$	-3.81$\times {10}^5$	9.30	-27.91	-229.2	6.25	128	2	32	32	32	20	3
JS3	0.3	-9.92$\times {10}^5$	-3.81$\times {10}^5$	9.30	-77.50	-229.2	6.25	128	2	32	32	32	20	5
JS4	0.3	-3.97$\times {10}^6$	-3.81$\times {10}^5$	9.30	-310.01	-229.2	6.25	256	2	64	64	64	20	10
JS5	0.3	-1.59$\times {10}^7$	-3.81$\times {10}^5$	9.30	-1240.1	-229.2	6.25	256	2	64	64	64	20	20
JXS1+	0.1	-1.76$\times {10}^6$	-2.37$\times {10}^6$	2.55	-137.78	-647.9	19.25	128	6	64	32	64	20	20
JXS2	0.1	-4.81$\times {10}^6$	-2.37$\times {10}^6$	2.55	-24.80	-647.9	19.25	128	6	64	32	64	40	3

Table 5. Table 5: Statistical quantities which are averaged over the time interval τ a ​ v subscript 𝜏 𝑎 𝑣 \tau_{av} . We list root mean square velocities in all three layers as well as the Nusselt number. The values in column RMC and RC (Ma=0) belong to the DNS which are listed in Tab.  4 . The last column gives the velocity unit for a system with dimension listed in Tab.  4

Case	$u_{rms}^{\left(1\right)}$		$u_{rms}^{\left(2\right)}$		$u_{rms}^{\left(3\right)}$		$\mathrm{Nu}$		$U_{vis}$
	RMC	RC	RMC	RC	RMC	RC	RMC	RC	in ${10}^{-6}$ m/s
H2	4.26	3.80	32.37	31.46	10.92	9.97	1.49	1.45	12.9
H3	15.32	15.26	87.14	86.56	31.20	31.05	2.22	2.12	8.6
H4	31.40	37.14	163.23	165.05	56.13	62.26	2.88	2.67	6.5
H5	38.86	35.81	380.43	385.30	226.43	170.53	4.45	4.61	4.3
H6	81.07	92.44	704.04	709.52	1180.94	1359.30	6.13	6.04	3.2
JM1	2.60	2.06	14.18	13.26	4.81	4.19	1.07	1.08	6.5
JM2	7.44	6.67	67.12	65.17	17.32	15.84	1.46	1.42	6.5
JM3	13.94	13.86	111.51	108.37	31.27	30.25	1.79	1.69	6.5
JM4	19.09	29.47	208.71	210.72	58.32	141.46	2.48	2.15	6.5
JS3	6.88	6.30	51.79	48.57	13.19	11.56	1.19	1.16	6.5
JS4	11.72	8.78	113.57	104.54	26.20	20.48	1.48	1.39	6.5
JS5	43.63	26.71	281.00	218.57	478.55	582.28	2.10	1.79	6.5

Equations

$$L_{vis}=d^{(1)},{\tau }_{vis}=\frac{{\left(d^{(1)}\right)}^2}{{\nu }^{(1)}},U_{vis}=\frac{{\nu }^{(1)}}{d^{(1)}},$$

$$\Theta =\frac{Q^{(2)}(d^{(2)}{)}^2}{8{\lambda }^{(2)}}=\frac{j_0^2(d^{(2)}{)}^2}{8{\lambda }^{(2)}{\sigma }_e^{(2)}},$$

$$P_{vis}={\rho }_{ref}^{(1)}{U}_{vis}^2=\frac{{\rho }_{ref}^{(1)}{\left({\nu }^{(1)}\right)}^2}{{\left(d^{(1)}\right)}^2}.$$

$${\varphi }_{ij}=\frac{{\varphi }^{(i)}}{{\varphi }^{(j)}}$$

$${\partial }_t{\mathbf{u}}^{(1)}=$$

$$-G{T}^{(1)}{\mathbf{e}}_z,$$

$${\partial }_t{\mathbf{u}}^{(2)}=$$

$$-G{\beta }_{21}{T}^{(2)}{\mathbf{e}}_z,$$

$${\partial }_t{\mathbf{u}}^{(3)}=$$

$$-G{\beta }_{31}{T}^{(3)}{\mathbf{e}}_z,$$

$$G=\frac{g{\beta }_T^{(1)}{\left(d^{(1)}\right)}^3\Theta }{{\left({\nu }^{(1)}\right)}^2}=\frac{g{\beta }_T^{(1)}{j}_0^2{\left(d^{(1)}\right)}^3{\left(d^{(2)}\right)}^2}{8{\left({\nu }^{(1)}\right)}^2{\lambda }^{(2)}{\sigma }_e^{(2)}}.$$

$$\nabla \cdot {\mathbf{u}}^{(i)}$$

$${\partial }_t{T}^{(1)}=$$

$$\displaystyle+\frac{1}{\Pr^{(1)}}\Big{[}\nabla^{2}T^{(1)}+\dfrac{8\lambda_{21}\sigma_{e,21}}{d_{21}^{2}}\Big{]}\,,$$

$${\partial }_t{T}^{(2)}=$$

$$\displaystyle+\frac{\kappa_{21}}{\Pr^{(1)}}\Big{[}\nabla^{2}T^{(2)}+\dfrac{8}{d_{21}^{2}}\Big{]}\,,$$

$${\partial }_t{T}^{(3)}=$$

$$\displaystyle+\frac{\kappa_{31}}{\Pr^{(1)}}\Big{[}\nabla^{2}T^{(3)}+\dfrac{8\lambda_{23}\sigma_{e,23}}{d_{21}^{2}}\Big{]}\,.$$

$${\mathrm{Pr}}^{(i)}=\frac{{\nu }^{(i)}}{{\kappa }^{(i)}}.$$

$$T^{(3)}$$

$$T^{(1)}$$

$$u_x^{(1)}$$

$$u_z^{(1)}$$

$${\partial }_z{T}^{(1)}$$

$$T^{(1)}$$

$$\frac{\mathrm{Ma}}{{\mathrm{Pr}}^{(1)}}{\partial }_x{T}^{(1)}$$

$$\frac{\mathrm{Ma}}{{\mathrm{Pr}}^{(1)}}{\partial }_y{T}^{(1)}$$

$$u_x^{(2)}$$

$$u_z^{(2)}$$

$${\partial }_z{T}^{(2)}$$

$$T^{(2)}$$

$$\frac{\mathrm{Ma}\xi }{{\mathrm{Pr}}^{(1)}{\mu }_{21}}{\partial }_x{T}^{(2)}$$

$$\frac{\mathrm{Ma}\xi }{{\mathrm{Pr}}^{(1)}{\mu }_{21}}{\partial }_y{T}^{(2)}$$

$${\sigma }^I$$

$${\sigma }^{II}$$

$$\mathrm{Ma}=\frac{{\alpha }_{\mathrm{T}}^{\mathrm{I}}{\mathrm{j}}_0^2{\sigma }_{\mathrm{ref}}^{\mathrm{I}}{\left({\mathrm{d}}^{(2)}\right)}^2{\mathrm{d}}^{(1)}}{8{\mu }^{(1)}{\kappa }^{(1)}{\lambda }^{(2)}{\sigma }_{\mathrm{e}}^{(2)}}$$

$$\xi =\frac{{\sigma }_{\text{mbox}ref}^I{\alpha }_T^I}{{\sigma }_{\text{mbox}ref}^{II}{\alpha }_T^{II}},$$

$${\sigma }^I$$

$${\sigma }^{II}$$

$$\left(\frac{{\mathrm{d}}^2}{\mathrm{d}{z}^2}-{\gamma }^{(i)}\right)q^{(i)}(z)=f^{(i)}(z)$$

$$T_{cond}^{(1)}(z)$$

$$T_{cond}^{(2)}(z)$$

Full text

I Introduction

In the wake of the rapid growth of renewable energies, the intermediate storage of energy is of central importance. Besides thermal and mechanical methods, chemical energy storage in liquid metal batteries (LMB) is a promising way that received increasing attention in recent years Kim et al. (2012); Wang et al. (2014). A liquid metal battery consists of three stratified liquid layers. The heaviest layer is the cathode, which consists of a liquid metal alloy. It is separated from the liquid-metal anode at the top by a molten salt layer rather than by an ion-permeable solid separator as in standard batteries Newman and Thomas-Alyea (2012). This implies that the operating temperature of present prototypes is several hundreds degrees Celsius compared to room temperature in standard devices. The feasibility of such configurations was demonstrated recently in the laboratory. Among them are electrode combinations of lithium and lead-antimony alloys (Li||Pb-Sb) (Wang et al., 2014) or magnesium and antimony (Mg||Sb) (Bradwell et al., 2012), respectively.

LMBs are multi-physics fluid systems that couple thermal effects, magnetic fields and electrochemical reactions to turbulent flows in each of the three layers. Different aspects of LMB systems have been studied in the last years. The current-driven Tayler instability, which may cause an electrical short-circuit between both electrodes for large systems with diameters of the order of a meter, has been investigated by Weber et al. Weber et al. (2013, 2015) and Herreman et al. (2015). The metal pad instability, which is already known from aluminum reduction cells, has been studied by Zikanov Zikanov (2015). Turbulent mixing processes in the liquid metal electrodes due to heating have been investigated experimentally by Kelley and Sadoway (2014). Very recently, thermal convection in the three-layer liquid metal battery was studied numerically by Shen and Zikanov Shen and Zikanov (2016). The authors showed in this work that in the presence of Joule heating, convective motion is always triggered for the typical Rayleigh or Grashof numbers in practical configurations. It can be expected that a significant part of the open physical questions and problems arises at the upper and lower interfaces between the electrodes and the molten-salt layer. This provides the motivation for our study of interfacial convection in LMBs. Moreover, the temperature distribution is an important optimization parameter in the operation of LMBs since a high temperature leads to increased energy losses while a low temperature increases the likelihood of short circuits by solidification.

In the present work, we want to investigate the effects of interfacial tension in a LMB. We ignore other processes in LMBs such as mass transport and chemical reactions due to modeling uncertainties and in order to keep the number of equations and parameters manageable. More precisely, we will extend the three-layer thermal convection setting by Marangoni effects and study the linear stability and the full nonlinear evolution of such a LMB model. Our three-dimensional numerical simulations (DNS) will show that thermal Marangoni convection is an important factor for the flow inside a LMB. In this respect, we extend the model of Shen and Zikanov (2016) to account for the differences of all transport coefficients and the gradient of interfacial tension. By doing this, interfaces are coupled by the continuity of all three velocities components across the interface and the balance of tangential viscous stresses with interfacial tension gradients. Effects of the electrical current and the resulting Ohmic heating are included by different heating rates in the three layers, which reflect the differences in electrical conductivity. As a consequence, temperature gradients are produced that potentially cause interfacial tension driven flows (i.e. Marangoni convection) and buoyancy driven flows (i.e. Rayleigh-Bénard convection). Our numerical simulations will allow us to study these effects either in combination or independently, whereby their relative importance can be estimated.

The deformations of layers, which have been incorporated, e.g., in magnetohydrodynamic numerical simulations of Herreman et al. (2015), are neglected. The azimuthal magnetic field, which is caused by the charge current across the cell is also neglected. This simplification is in line with the results of Shen and Zikanov (2016): The azimuthal magnetic field of a homogeneous current remains negligible for laboratory sized configurations which are in the focus of the present work, particularly in view to the nonlinear evolution of the three-layer model. The typical layer height is then of the order of centimeters.

Also, we will consider an internal fraction of the triple-layer only and exclude side walls. This will allow us to use the simpler Cartesian geometry in combination with periodic boundary conditions in the horizontal directions. The latter step opens the possibility to apply exponentially fast converging pseudospectral simulations Köllner (2016). Finally, as in all previous studies, the model is considerably simplified by disregarding the mass transport of metal from the top layer to the bottom layer and the resulting change of density and interfacial tension with the composition.

Thermophysical properties, especially for molten salts, are difficult to obtain from the literature for currently investigated LMBs. Thus we select the following representative substances for the LMB configuration:

• •

For the upper electrode, we take lithium with a melting point at 181∘C. Lithium is a typical metal with favorable properties to construct LMBs (Kim et al., 2012). Material properties of lithium are taken from ref. Maroni et al. (1973).

• •

For the layer that separates both liquid metal electrodes, a lithium chloride (59 mol%)-potassium chloride (41 mol%) eutectic mixture is taken which is denoted as LiCl-KCl. The melting point is at 355∘C. Most of the required transport properties of eutectic LiCl-KCl are collected in Williams et al. (2006). The surface tension value can be found in ref. Smirnov and Stepanov (1982) and the electrical conductivity in Van Artsdalen and Yaffe (1955). The mass density as a function of temperature is given in ref. Smirnov and Stepanov (1982). It is noted that in Shen and Zikanov (2016) the same salt used.

• •

For the lower electrode, we take eutectic lead-bismuth (Pb-Bi), which has been also used in the mixing experiments by Kelley and Sadoway (2014). The material properties of eutectic Pb-Bi are also collected in ref. Fazio et al. (2015).

As a typical temperature in an experiment (Wang et al., 2014) we take $T^{\infty}=500^{\circ}$C (773.15K), which is far above the melting temperature of all three battery components. The material properties at this operating temperature are collected in Tab.1. The interfacial tension between layers is estimated by the method of Girifalco and Good (1957), see Sec. II.3.

Our paper is organized as follows. Section II presents the non-dimensional model equations and the numerical method. Sec. III presents the linear stability analysis for the triple-layer configuration. Four convection regimes are identified: Rayleigh-Bénard convection (RC) in the middle and top layers, Marangoni convection (MC) at each of the two interfaces. Based on linear stability investigations, Sec. IV discusses three-dimensional direct numerical simulation (DNS) of the full nonlinear evolution at given material parameters, which represent typical laboratory-scale configurations in terms of currents and layer heights. In Sec. V, we derive a practical Rayleigh number criterion to decide when each of the four different regimes can appear and/or will dominate. The anti-convection mode (AC) in the bottom electrode is also considered. Moreover, the relevance of interfacial deformations caused by the convection in the layers is estimated. Finally, we summarize the present work in Sec. VII and give a brief outlook.

II Three-layer model and numerical methods

II.1 Characteristic units

The model describes three immiscible Newtonian liquids that are stably stratified due to their density difference. Figure 1 sketches this three-layer model in the non-dimensional formulation. The lower layer, representing a dense liquid-metal alloy, is located in the range $-d^{(1)}<z<0$ and related quantities will be denoted with a superscript $(1)$. The molten salt layer occupies the interval $0<z<d^{(2)}$. On top, a layer of liquid metal is located with $d^{(2)}<z<d^{(2)}+d^{(3)}$. Layer $\Omega^{(1)}$ and $\Omega^{(2)}$ are separated by the planar interface I, layer $\Omega^{(2)}$ and $\Omega^{(3)}$ by the interface II, respectively.

The three layers are bounded by solid walls at the bottom and top, which are held at a uniform temperature $T^{\infty}$. A homogenous internal heating is applied in each of the three layers, which mimics a prescribed homogeneous current density $j_{0}\mathbf{e}_{z}$. The differences in the electrical conductivity $\sigma_{e}^{(i)}$ are also taken into account. Convection is driven by gradients in mass density $\rho^{(i)}(T)$ as well as the interfacial tensions $\sigma^{I}(T)$ and $\sigma^{II}(T)$, all of which depend on the temperature field $T(\mathbf{x},t)$.

The characteristic units for time, $\tau_{vis}$, length, $L_{vis}$ and velocity, $U_{vis}$ are based on the height of the bottom layer and the characteristic time of viscous equilibration across this layer. They are given by

[TABLE]

with $\nu^{(1)}$ being the kinematic viscosity in layer $\Omega^{(1)}$. The lateral aspect ratios in $x$– and $y$–directions are always equal, i.e. $l_{x}=l_{y}$. The appropriate temperature unit $\Theta$ is chosen to represent the maximum temperature appearing at pure conduction in the middle layer (Kulacki and Goldstein, 1972). It is given by

[TABLE]

where $Q^{(2)}=j_{0}^{2}/\sigma_{e}^{(2)}$ is the volumetric Joule dissipation rate (measured in J/(s m3)) in the middle layer and $\lambda$ a thermal conductivity. Also, $j_{0}$ is the constant current density and $\sigma_{e}$ is the electrical conductivity. Dimensionless temperatures are thus given by $T=(\tilde{T}-T^{\infty})/\Theta$ where $\tilde{T}$ is the physical temperature. In Sec. III.1, the case of pure conduction is solved for the present three-layer problem. It will turn out that a zero heating rate in the top and bottom layer provides a good approximation to the exact conduction solution since the electrical conductivity of the middle molten salt layer is considerably lower than in the outer liquid metal layers. Furthermore, the unit of temperature in eq. (2) will be a good estimate for the maximum temperature in the cell. Finally, the unit for pressure is given by

[TABLE]

From now on, we will consider dimensionless units only. The ratios of material parameters $\phi\in\{\nu,\kappa,\sigma,\lambda,\beta_{T},...\}$ are denoted by the following abbreviation

[TABLE]

They are summarized in Tab. 2. This will simplify the notation of the set of equations which are discussed next.

II.2 Equations and boundary conditions

In the following, we list the balance equations of momentum, mass and energy in nondimensional form. These equations have been successfully used to simulate Rayleigh-Bénard or Marangoni convection in other configurations (Nepomnyashchy et al., 2006; Chillà and Schumacher, 2012). Since the magnitude of all material parameters are related to those of the bottom layer, material parameter ratios appear additionally. The transport of momentum is given by the incompressible Navier-Stokes equations in the Boussinesq approximation (Nepomnyashchy et al., 2006; Chillà and Schumacher, 2012)

[TABLE]

with the Grashof number

[TABLE]

Here, $\mathbf{g}=(0,0,-g)$ is the vector of acceleration due to gravity and $\beta_{T}$ is the thermal expansion coefficient. The mass balance is given by

[TABLE]

with $i=1,2,3$. The variable $\mathbf{u}$ denotes the velocity field. The energy balance reduces to a transport equation for the temperature field in each layer

[TABLE]

The last term in each of the previous equations describes the non-dimensional volumetric heating rate. The Prandtl number of the layer (i) is given by

[TABLE]

The Prandtl number of the middle layer is comparable to water whereas the metal layers have very small Prandtl numbers of order $10^{-2}$, see Tab. 2.

No-slip and isothermal boundary conditions are imposed for the solid walls at the bottom and top:

[TABLE]

The matching conditions at the lower planar interface I at $z=0$ are as follows

[TABLE]

with the dynamic viscosities $\mu^{(i)}=\rho^{(i)}\nu^{(i)}$. At the upper planar interface II at $z=d_{21}$, the matching conditions are as follows

[TABLE]

These conditions (16)-(27) enforce the continuity of velocity components across the interfaces I and II, the balance of the tangential stresses with interfacial tension gradients, and a planar interface, which is an approximation to a more general normal-stress balance. The variation of interfacial tension with temperature is given in dimensional form by

[TABLE]

with $\alpha_{T}^{I}$ and $\alpha_{T}^{II}$ being the interfacial tension coefficients at both interfaces. The Marangoni number is defined with respect to interface I and given by

[TABLE]

The ratio of interfacial tension change with temperature between the upper and the lower interface is quantified by

[TABLE]

such that $Ma^{II}=\xi Ma$. Table 2 summarizes all parameter values and ratios. The lateral boundary conditions are periodic, which is required by expanding the solution into a sum of Fourier modes. Simulations are typically started with random velocity and zero temperature fields, representing the sudden switch-on of an electrical current.

II.3 Estimation of interfacial tension

Most of the material properties required for our model are well documented. The appropriate references are given in Sec. I. However, the interfacial tension, a property that represents the molecular interaction between particles in the electrode and the molten salt (Rowlinson and Widom, 1989), is unknown for the particular case. Laboratory measurements of interfacial tension between liquid metals and molten salts are to the best of our knowledge not existent for the present configuration.

To obtain interfacial tension, the method of Girifalco and Good (1957) is employed here. It calculates interfacial tension as the sum of the well-known surface tensions of each component and subtracts the energy of the intermolecular bonds across the interface. Formally, it reads

[TABLE]

with the exchange parameter $\Phi\leq 1$. Values of $\Phi$ smaller than one account for molecular interactions that are different from dispersion type forces (Berg, 2010). Especially, for the present combination of ionic and metallic bonds, this parameter can be considerably different from unity. For the interface of aluminum and a salt mixture $\Phi=0.42$ has been noted by Roy and Utigard (1998)). In the following, we take a value of $\Phi=0.7$.

The following surface tension values $\sigma^{(i)}$ of each phase at an operating temperature of 500 ∘C are extracted from the references in Sec. I. The particular values are $\sigma^{(1)}=0.387$ N/m [$\partial_{T}\sigma^{(1)}=7.99\times 10^{-5}N/(mK)$] for Pb-Bi, $\sigma^{(2)}=0.1327$ N/m [$\partial_{T}\sigma^{(2)}=8.26\times 10^{-5}N/(mK)$] for LiCl-KCl and $\sigma^{(3)}=0.3493$ N/m [$\partial_{T}\sigma^{(3)}=1.6\times 10^{-4}N/(mK)$] for Li in the top electrode. With these values and Eqns. (32,33) all free parameters in Eqns. (28) and (29) can be calculated (see also Tab. 1). Note also that the decrease of interfacial tension with temperature is in line with the fact that the mutual solubility of salt and liquid metal typically increases with temperature Kim et al. (2012).

II.4 Numerical method

The numerical method is based on a pseudospectral algorithm using a Fourier expansion of all fields in $x$– and $y$–directions with $N_{x}$ and $N_{y}$ modes, respectively. No-slip boundary conditions at the top and bottom as well as nonpenetrative boundary conditions at the interfaces require a Chebyshev polynomial expansion with respect to $z$–direction with a polynomial degree of $N_{z}^{(i)}$ in each layer $\Omega^{(i)}$. For the polynomial degrees, we only use powers of two because of the fast Fourier transformations. The time-stepping scheme is a combination of the implicit Euler backward formula for linear terms and the Adams-Bashforth formula for nonlinear terms. The algorithm is a straightforward extension of the solver developed by Boeck et al. (2002), Köllner et al. (2013) and Köllner (2016). The time step size is adapted such that the Courant-Friedrichs-Lewy number is in the range from 0.1 to 0.2.

The main difference to our former simulations (Köllner et al., 2013, 2016) is that now three instead of two layers are coupled. After discretization with respect to the horizontal $x$– and $y$–directions as well as with respect to time, each individual Fourier mode of a hydrodynamic field satisfies a one-dimensional inhomogeneous Helmholtz equation

[TABLE]

with $\gamma^{(i)}$ being a real number, $f^{(i)}(z)$ a given complex-valued given function, and $q^{(i)}(z)$ an unknown complex-valued function. These three coupled equations (for $i\in\{1,2,3\}$) are solved directly with the Chebyshev tau method Canuto et al. (1988) in the same way as described for two layers by Boeck et al. (2002). The solver is written in the C language and parallelized using the Message Passing Interface library. The actual resolutions used are specified below.

The solver is validated by comparing with simple test cases, i.e., steady conduction and the Poiseuille flow that is driven by a homogeneous volume force in $x$-direction. Further checks were done by reproducing the stability threshold for the case treated by Géoris et al. (1999). We also reproduced the simulation results of Shen and Zikanov (2016), which however could be compared only qualitatively, because of different lateral boundary conditions. Nevertheless, flow structures were found to be rather similar. Furthermore, we checked that the kinetic energy balance holds and successfully reproduced the linear stability thresholds of convection, which is discussed in more detail in Sec. III.

III Linear stability analysis

In the following, we will investigate the linear stability of the stationary pure conduction state. The closed form of the stationary temperature distribution is detailed in the following Sec. III.1. Thereafter, the linear stability problem is formulated (Sec. III.2), its solution procedure is discussed, and critical Marangoni and Grashof numbers for the present reference system are calculated (Sec.III.3).

III.1 Temperature profile of pure conduction

The case of pure conduction ($\mathbf{u}=0$) with initial conditions $T^{(i)}(t=0)=0$ converges to a time-independent solution for $t\rightarrow\infty$ and can be found analytically. To further simplify the present problem, we neglect the Joule dissipation in the liquid metal electrodes and proceed with the assumption $\sigma_{e}^{(1)}=\sigma_{e}^{(3)}=\infty$, which will be justified later by comparison of the stability predictions with DNS of the fully nonlinear equations of motion. The steady state temperature distribution is a solution of Eqns. (10)–(12) and given by

[TABLE]

where we introduced the abbreviations $A=d_{21}\lambda_{32}+\lambda_{31}+d_{31}$, $B=d_{21}\lambda_{32}+2d_{31}$. This temperature distribution is plotted together with the numerical solution that accounts for Joule dissipation in the liquid metal layers in Fig. 2. One can readily observe that the top layer has higher conductivity than the lower one, resulting in a smaller temperature gradient in $\Omega^{(3)}$. Numerical and analytical solutions agree perfectly.

The temperature at both interfaces is given by

[TABLE]

In the reference case, which is displayed in Fig. 2, the interfacial temperatures are $T^{I}_{cond}=0.0995$ and $T^{II}_{cond}=0.0296$. The maximum temperature appears in the middle layer. It is not exactly at $z=0.5$, but shifted slightly towards the layer with lower thermal conductivity. It occurs at

[TABLE]

with a value of

[TABLE]

For the reference case, this results in $z_{max}=0.49$ and $T_{max}=1.06$.

III.2 Linearized equations and solution method

We proceed now with an extension of the linear stability analysis of Rayleigh-Marangoni convection to a three-layer model. Former studies with a temperature difference applied between the boundaries can be found in Géoris et al. (1993); Nepomnyashchy et al. (2005, 2006). The vertical velocity component and the temperature perturbation are expanded into normal modes Chandrasekhar (1961)

[TABLE]

Quantities $k_{x},k_{y}$ are components of the horizontal wavenumber vector with magnitude $k=(k_{x}^{2}+k_{y}^{2})^{1/2}$. The governing temperature equations are linearized around the basic state $T_{cond}(z)$. Furthermore, by applying the curl two times to the linearized Navier-Stokes equations, the following system of linear equations (Nepomnyashchy et al., 2006) is derived

[TABLE]

Here $D=d/dz$. For the temperature it follows that

[TABLE]

The matching and boundary conditions at $z=0$ are given by

[TABLE]

and at $z=d_{21}$ by

[TABLE]

The boundary conditions at the no-slip walls are as follows. For $z=-1$

[TABLE]

and at $z=d_{21}+d_{31}$

[TABLE]

Furthermore, we set $\rm{Ma}^{I}=\rm{Ma}$ and $\rm{Ma}^{II}=\xi\rm{Ma}$ in Eqns. (52) and (57), in order to study the effect of Marangoni convection on both interfaces separately.

The resulting linear stability problem takes the form

[TABLE]

where the column vector $\mathbf{q}$ consists of the independent fields $\{w^{(3)},w^{(2)},w^{(1)},\theta^{(3)},\theta^{(2)},\theta^{(1)}\}$ and the linear differential operators $\bm{M}$ and $\bm{L}$ encode the bulk equations and the boundary conditions. A general solution can be expressed by the exponential of the operator $(\bm{M}^{-1}\bm{L})$ (Trefethen and Embree, 2005), which can be further analyzed to derive several properties of the solution. Also, one can expand the $\mathbf{q}(z,t)$ into the eigenfunctions $\hat{\mathbf{q}_{i}}(z)$ of $(\bm{M}^{-1}\bm{L})$ by $\mathbf{q}=\sum_{i}\hat{\mathbf{q}_{i}}(z)\exp(ts_{i})$ with complex growth rates $s_{i}$. We are interested in the marginal stability properties only where $s_{i}=0$. Then, the system can be rewritten as an eigenvalue problem with one control parameter.

For the analysis of the onset of Marangoni convection, i.e. marginal stability, eq. (62) with $\partial_{t}=0$ is discretized with the Chebyshev collocation method Schmid and Henningson (2012); Trefethen (2000) and rearranged thereafter to an eigenvalue problem with respect to the Marangoni number Ma. The result is a generalized eigenvalue problem with quadratic matrices $\bm{A}$ and $\bm{B}$ of dimension $6(N+1)\times 6(N+1)$, which follows to

[TABLE]

Here, $N$ is the degree of the Chebyshev polynomials which are used to discretize the vertical direction in each layer. It is usually set to $N=32$.

For a given set of parameters, the critical Marangoni number $\rm{Ma}_{c}$ is the eigenvalue to Eq. (63) with the minimal magnitude over all wavenumbers $k$ of the perturbations. Note that we only consider situations for which density and interfacial tension decrease with increasing temperature. This is equivalent to $\rm{Ma}\leq 0$ and $G\leq 0$. Eigenvalues are calculated with Matlab.

III.3 Results of linear stability analysis

III.3.1 Pure Marangoni instability ($\rm{Ma}\neq 0$, $G=0$)

The neutral stability curve for the reference case with the parameters taken from Tab. 2 is displayed in Fig. 3(a). We find a critical value of $\rm{Ma}^{I}_{c}=-76.39$ at a wavenumber of $k_{c}=2.10$. The marginal values of $\rm{Ma}^{I}$ (note that the coupling $\rm{Ma}^{II}=\xi\rm{Ma}^{I}$ applies here implicitly) were calculated for different polynomial degrees $N$. The critical value does not change up to the fourth digit.

According to the critical value $\rm{Ma}^{I}_{c}$, two nonlinear simulation are performed at $\rm{Ma}^{I}=-76$ (LS01) and $\rm{Ma}^{I}=-77$ (LS02), respectively (see Tab. 3). Their root mean square (rms) velocity which is defined by

[TABLE]

is plotted in the bottom panel of Fig. 3 as function of time, successfully verifying the prediction from the eigenvalue calculation. In Fig. 4, the eigenfunctions (crosses) at the critical parameters $k_{c},\rm{Ma}_{c}^{I}$ are plotted together with the respective simulation profiles (solid lines). Pointwise differences between two snapshots are compared. Eigenfunctions have been rescaled to have a range of values which is comparable with the simulation results. Convection is caused at both interfaces since the temperature perturbation amplitude is non-zero there. At a given horizontal position, the temperature may be locally increased at the upper interface ($\theta(z=d_{21})>0$) causing a lower interfacial tension which drives a divergent flow at this point, i.e., a flow directed towards the interface with $\partial_{z}w(d_{21})<0$. At the same horizontal position, interfacial tension is then increased on the lower interface, and the flow is directed away from the interface. At both interfaces work is performed on the liquid metal phases.

However, the upper interface II induces a stronger flow than the lower interface I as seen in the bottom panel of Fig. 4. To understand the impact of each interface separately, we calculated the critical Marangoni number for each interface by disabling the effect of the other. When the lower interface is disabled, i.e., $\rm{Ma}^{I}=0$, the value of $\rm{Ma}^{II}_{c}=-233.55$ is obtained which corresponds with $\rm{Ma}^{II}_{c}/\xi=-90.87$, at a critical wavenumber $k_{c}=2.2$ which is very close to the value for both interfaces being active. When the upper interface is disabled, i.e., $\rm{Ma}^{II}=0$, one obtains $\rm{Ma}^{I}_{c}=-106.88$ at $k_{c}=2.25$. The latter case of $\rm{Ma}^{II}=0$, is also verified by two DNS, one at $\rm{Ma}^{I}=-105$, showing decay and $\rm{Ma}^{II}=-110$, showing a growth of perturbations. Consequently, both interfaces collectively contribute to the instability. Moreover, the eigenfunctions with one active interface (not displayed) look similar to the eigenfunctions that result from two active interfaces.

III.3.2 Rayleigh-Marangoni instability ($\rm{Ma}\neq 0$, $G\neq 0$)

Next, a closer look at the impact of buoyancy effects with $G\neq 0$ is taken. In order to do so, we calculate the critical Marangoni number $\rm{Ma}^{I}_{c}$ as a function of the Grashof number $G$. Figure. 5(a) shows this computed relation as crosses, while Fig. 5(b) displays the corresponding critical wavenumber. The magnitude of the critical Marangoni number decreases as |$G$| is increased. In what follows, the relation $\rm{Ma}^{II}=\xi\rm{Ma}^{I}$ will always hold. Thus, we will continue with the single Marangoni number $\rm{Ma}\equiv\rm{Ma}^{I}$ only. An approximately linear trend of the form,

[TABLE]

is observed, where $\tilde{\rm}{Ma}_{c}$ and $\tilde{G}_{c}$ are the critical numbers in absence of the correspondingly other physical process. Such linear scaling has been already observed in the one-layer case (Nield, 1964) when both, the Marangoni and the Rayleigh instability, destabilize the system.

The particular value of the critical Grashof number $\tilde{G}_{c}$ is calculated from an eigenvalue problem similar to (63) but now for Grashof instead of Marangoni number. We find a critical value $\tilde{G}_{c}=-1.2937\times 10^{4}$ with a critical wavenumber of $k_{c}=3$. The eigenfunctions of the pure Rayleigh-Bénard case are displayed in Fig. 6. They show a similar flow structure as in the Marangoni case, though now, caused by the downwelling of colder fluid in the upper half of the mid layer. The flow in the top layer is caused by viscous friction at the interface. In this layer, temperature perturbations are negative which result in a stabilization.

Again we compared these findings to two DNS simulations, namely DNS runs LS03 and LS04. After a decay of initial perturbations, the unstable modes start to develop in the supercritical case LS04 while perturbations decays in the subcritical case LS03.

The middle layer Rayleigh-Bénard mode leads to eigenfunctions that are qualitatively equal to those of Marangoni convection. As a consequence, Rayleigh-Bénard (RC) and Marangoni convection (MC) act together to destabilize the system at hand. The case when both sources for convection are present will be denoted by Rayleigh-Marangoni convection (RMC). The Rayleigh and Marangoni modes act together in the middle layer of the present three-layer system as in the single layer case. The actual thresholds can be calculated by Eq. (65).

III.3.3 Middle layer thickness

While a change in the electrical current density magnitude and the overall system height at fixed ratios $d_{21}$ and $d_{31}$ is accounted for in the Marangoni and the Grashof numbers, respectively, the stability threshold is changed when the mid-layer thickness $d_{21}$ becomes smaller compared to the outer layer, i.e., $d_{21}<1$ in our notation. Therefore, we also computed the stability thresholds for both convection mechanisms as a function of the middle layer height $d_{21}$.

Figure 7 shows how the magnitude of $\tilde{\rm}{Ma}_{c}$ grows with decreasing layer height, which is depicted with the compensated plot of $\tilde{\rm}{Ma}_{c}\,d_{21}$ as a function of $d_{21}$. In this case, $G=0$ was taken. For example, the threshold values for $d_{21}=0.3$ are $\rm{Ma}_{c}=229.2$ and $k_{c}=6.3$. For this case, we verified again by two DNS, namely, LS05 and LS06, that this Marangoni number is indeed critical.

The magnitude of the critical Grashof number ($\rm{Ma}=0$) increases as well with decreasing middle layer height. The critical Grashof number which is compensated with the third power of $d_{21}$ is shown in Fig. 8. In this compensated representation, one observes that the Grashof number increases less than the third power of $d_{21}$. For the shallowest layer $d_{21}=0.1$ one observes a particularly strong decrease of $\tilde{G}_{c}\,(d_{21})^{3}$.

An inspection of the corresponding eigenfunctions reveals how the velocity perturbation changes with $d_{21}$. The results are shown in Fig. 9. For a shallow middle layer with $d_{21}=0.3$, the upper layer, which is initially coupled by viscous shear only, gets active and develops internal convection rolls. This is because the coupling of layers by viscous stresses and heat transport are not mutually ’compatible’ due to the much higher thermal diffusivity in the liquid metal top layer. For $d_{21}=0.1$, convection in the upper-layer layer dominates. A corresponding upper layer Rayleigh number can be calculated by

[TABLE]

with $T^{II}=(\tilde{T}^{II}-T^{\infty})/\Theta$ (see Eq. (8)). For the calculated threshold ($d_{21}=0.1$) of $G_{c}=-2.37\times 10^{6}$, $k_{c}=2.55$, this number equals to $\rm{Ra}^{(3)}=-1290.4$. This should be compared to the classical linear stability analysis in one layer where the critical Rayleigh number ranges from $\rm{Ra}_{c}=120$ with adiabatic temperature ($D\theta=0$) and free-slip velocity ($w=D^{2}w=0$) boundary conditions to a value of $\rm{Ra}_{c}=1707.76$ with isothermal temperature ($\theta=0$) and no-slip velocity ($w=Dw=0$) boundary conditions (Colinet et al., 2001). Hence for shallow middle layers a fourth regime of convection can be expected that is similar to the classical Rayleigh-Bénard convection, caused by the linear temperature gradient in the upper layer. Note also that we probed successfully the stability threshold for $d_{21}=0.1,0.3$ with DNS simulation runs LS07 to LS10 (see again Tab. 3).

Finally, it is investigated how both effects work together, in the same way as in the previous Sec. III.3.2. It is found that relation (65) can be applied up to $d_{21}\geq 0.2$. However, for $d_{21}=0.1$, we found two different behaviors of $G_{c}(\rm{Ma})$: first, for $|\rm{Ma}|<400$, the Marangoni effect counteracts the upper layer RC, such that the critical Grashof number grows in magnitude with increasing $|\rm{Ma}|$. Second, for |Ma|>400, the middle layer RC is dominant, thus, $|G_{c}|$ is decreasing with increasing $|\rm{Ma}|$.

III.3.4 Convection patterns near the threshold values

The characteristic structures of convection near the onset threshold are shown in two DNS which are denoted LS11 and LS12, respectively. Here, we used the reference current $j_{0}=3$kA/m2, but adapted the height to the threshold of the onset of Rayleigh-Bénard convection. Therefore, all three layer heights are set to $d^{(1)}=6.4$ mm, which results in an overcritical Grashof number (see DNS run LS11 in Tab. 3). Moreover, a simulation LS12 with $\rm{Ma}=0$, but otherwise the same parameters as LS11 was conducted.

The rms velocity $u_{rms}^{(2)}(t)$ for the two simulations are shown in Fig. 10. The solid line represents Rayleigh-Bénard convection and the dashed line Rayleigh-Marangoni convection. We observe that the Marangoni effect accelerates the growth of the rms velocity amplitude and thus of convection. The larger growth rate is expected from our former results. The flow structures of these simulations are visualized in Fig. 11(a-d). We show vertical cuts at $y=0$ with the temperature and velocity for RC in panel (a) of the figure and RMC in panel (c). Furthermore, the interfacial temperature and velocity vectors at the same time instant are displayed in panels (b,d). A nearly hexagonal and stationary convection pattern appears with an outflow in the center of cells. This is similar to observations for one layer with internal convection Kulacki and Goldstein (1972); Tveitereid and Palm (1976). In conclusion, the former results that the Rayleigh and Marangoni effect act together in the exemplary configuration are hereby again confirmed.

IV Nonlinear evolution

In the following section we will present a series of DNS, listed in Tab. 4, that investigate the full nonlinear and turbulent evolution of the RMC in the three-layer battery model. Several parameters will be varied for this purpose. Runs H1 to H6 are conducted at $d_{21}=1$ and fixed $j_{0}$, but different dimensional heights $d^{(1)}$. An increase of $d^{(1)}$ is in line with an increasing Grashof number. In the series JM1 to JM4, we fix the physical height of the bottom layer $d^{(1)}=20$ mm as well as $d_{21}=0.5$ while varying $j_{0}$ (J = varying current density, M = medium middle-layer height). The same holds for the series JS1 to JS5 where $d_{21}$ is further decreased to 0.3 in comparison to the JM series (S = small). Finally, the runs JXS1 and JXS2 set $d_{21}$ to 0.1 (XS = extra small).

IV.1 Dynamics of the reference case H4

As the reference case, we take the configuration with each layer having a height of 20 mm and a current density of 3kA/m2, (see case H4 in Tab. 4). All other parameters are given again in Tab. 2. Although for commercial applications shallower middle layers are preferable to decrease Ohmic losses in the separator, a thicker middle layer helps us here to reveal flow properties in the cell center and at the interfaces. To separate the different convection mechanisms from each other, we performed three simulations (1) RC at $\rm{Ma}=0$, (2) MC at $G=0$, and (3) RMC. Table 4 also lists the critical Marangoni and Grashof numbers as well as the corresponding critical wavenumbers which allow us to estimate the relative contribution of the effects. One finds a factor of supercriticality of $G/\tilde{G}_{c}=308$ for buoyancy–driven convection compared to $\rm{Ma}/\tilde{\rm}{Ma}_{c}=4$ for interfacial tension–driven convection.

The rms velocity in each of the three layers is plotted in Fig. 12 for RC (top panel), MC (middle panel), and RMC (bottom panel). After the onset of convection, all systems run into a statistically stationary state. For the MC case the flow velocity is weakest. The differences in the rms velocity between RMC and RC remain small, but become better visible when plotting the $x,y$-averaged vertical profiles, $\sqrt{\langle\mathbf{u}^{2}\rangle_{x,y,\tau_{av}}}$, which is done in the top panel of Fig. 13. The profiles are averaged over a time interval $\tau_{av}$, which is at least 10 convective time units long.

For RC and RMC, the primary flow structure is caused by buoyant convection in the middle layer. Interestingly for the RMC case, the rms velocities are smaller in the outer layers than for RC. For MC, the flow is driven by interfacial tension gradients and therefore rms velocity amplitudes are largest at the interfaces. The lower interface contributes less because of the ratio $\xi=2.57=\rm{Ma}^{II}/\rm{Ma}^{I}$. For these three cases, the middle layer velocity has on average always the highest amplitude while the rms velocity in the bottom electrode is always smallest. In cases RMC and RC, the small rms velocity magnitude is most probably caused by the stable density stratification in the lower layer. In case of MC, the upper-layer and middle-layer rms velocity are of similar magnitude, which is also in line with the specific driving mechanism, namely that interface II between molten salt and upper electrode is the main source for Marangoni convection.

The impact of convection on the $x,y$- and time-averaged, vertical temperature profiles is displayed in the bottom panel Fig. 13. The temperature is smallest in the case of RMC and highest in the case of MC. In connection to this observation, we measure the global transport of heat by the Nusselt number in the following way. Following the works on convection with internal heating in a single layer by Goluskin and van der Poel (2016), we define a Nusselt number Nu by the maximum of the $x,y$-averaged temperature, $T_{max}$. It is given by

[TABLE]

Furthermore, one takes the maximum temperature for pure conduction $T_{max}^{con}$ which is given by Eq. (41) and arrives at the definition

[TABLE]

of the Nusselt number, which will be used in the following. Figure 14 shows the Nusselt number as a function of time for the three convection cases. The difference between the three cases supports our observation for the temperature profiles in Fig. 13. The initial drop is caused by the initialization with zero temperature everywhere. The three-layer system requires a finite time to heat up. This heating time can be estimated from the balance between the time derivative of temperature and the heating rate in Eq. (11), which provides a heating period of about $\Pr^{(1)}{d_{21}}^{2}/(8\kappa_{21})\approx 0.1$ The reason of why RMC transports the heat more efficiently than RC can be answered by an examination of the typical flow structures, which is done next.

In Figs. 15, 16 and 17, we display snapshots of the flow and temperature from different perspectives for cases RC, MC and RMC, respectively. All figures are taken at time $t=0.3$. In the contour plots at the upper interface II as well as in the vertical cuts at $y=0$, we kept the same temperature range in all three cases. The differences were too large for the lower interface I such that we had to take different contour level intervals here.

Case RC is organized in irregular cells with an upwelling hotter fluid in the center and a downwelling colder fluid at the cell boundaries, as seen in all three two-dimensional cuts in Fig. 15. The downwelling colder fluid appears in form of sheets and jets, which sometimes carry fluid all the way down to the lower interface I across the middle layer $\Omega^{(2)}$. In the vertical cut at $y=0$, one also observes that the upper and lower layers are coupled by the viscous stresses to the mid layer since the velocity fields are aligned. The isosurface of $T^{(2)}=0.37$ in the bottom right panel completes the description of the convection structures. Grooves are formed by downwelling cold fluid. They resemble the typical patterns which are visible in the temperature at the upper interface, $T^{II}(x,y,d_{21},t=0.3)$.

Case MC which is displayed in Fig. 16 shows slightly different flow structures at the same time. The temperature distribution at the upper interface II depicts convection cells that are driven by the increase of interfacial tension from the cell center (hot) to the cell boundary (cold). The Marangoni effect also determines the flow behavior at the lower interface, i.e., the fluid streams from the hotter to colder regions. However, the structures across the middle layer are strongly influenced by the convection cell patterns that are formed at the upper interface. The isosurface of $T=0.5$ in the bottom right panel Fig. 16 encloses hotter fluid. It is seen how fluid is transported from the center of the middle layer to the upper interface where it cools down. The vertical cut of the figure complements the view on the convection structure in this case.

Case RMC resembles mostly case RC, which can be seen by the similar composition of temperature distributions in the top left panel of Fig. 17. However, there are some slight differences visible. The convection cells in the top left panel have a smaller typical length scale. This is the reason for the enhanced heat transport and the higher Nu which was shown in Fig. 14. Moreover, downwelling fluid appears preferentially in form of stronger jets rather than sheets. These structures couple the top and bottom electrode directly. The more frequent jets are also visible by the hot spots in the lower interface in the top right panel and in the isosurface in the lower right panel. This is because cold regions – the source of downwelling jets – have a tendency now to merge by the Marangoni effect. In conclusion, the flow pattern of RMC is of smaller length scale than the one for pure RC. This decrease of length scales is to our view the reason for the enhanced Nusselt number since the heat can be carried along a larger number of up- and downwelling structures.

IV.2 Variation of the layer height in H1 to H6

After the detailed description of the reference case H4, we turn in the following to the variation of the layer heights. In the simulation series H the ratios $d_{21}$ and $d_{31}$ will be kept to unity and the electrical current density remains fixed to 3kA/$m^{2}$ (see cases H1 - H6 in Tab. 4). The physical height $d^{(1)}$ will be changed. For case H1 with $d^{(1)}=5$ mm layer, the system is linearly stable and all initial infinitesimal perturbations decay with respect to time. As already discussed in Sec. IV.1, the simulations run always at least 10 convective time units for any statistical analysis. In addition to each of the RMC cases, which are given in Tab. 4, corresponding RC runs are also conducted.

Fig. 18 displays the rms velocity in each of the three layers. As already described for the reference case, the strongest flow is observed in the middle layer, except in case H6. For the tallest system, the upper layer shows the highest magnitude of rms velocity. The reason lies in the unstable stratification in the upper layer. It is strong enough to develop an additional large-scale convection cell (in the sense of Rayleigh-Bénard convection) that contributes dominantly to the rms magnitude. This behavior can also be rationalized by the Rayleigh number $\rm{Ra}^{(3)}$ for the upper layer $\Omega^{(3)}$ which takes values of $-193$ for H4, $-1469$ for H5 and $-6192$ for case H6, respectively. As visible from Eq. (66), case H6 is far beyond any linear instability threshold. Thus, a significant convective fluid motion can be expected to exist for this case only.

We also indicate a power law scaling of $G^{0.45}$ in Fig. 18 that matches the data for the rms velocities in the molten salt layer very well. All computed values are also listed in Tab. 5 in the appendix. The Reynolds number is a measure for the global transport of momentum in turbulent convection flows. The definition for each of the three layers in the present liquid metal battery model is given by

[TABLE]

Thus, the Reynolds numbers will exhibit exactly the same trend as the rms velocities in Fig. 18 and are therefore not shown. Only the relative amplitudes with respect to each other vary since the liquid metal has lower kinematic viscosity than the molten salt. In classical turbulent Rayleigh-Bénard convection with no-slip and isothermal walls for air Schumacher et al. (2014) or liquid metals Scheel and Schumacher (2016) as working fluids the Reynolds numbers obeys nearly the same scaling with increasing Grashof number. The scaling deviates slightly from the exponent of 1/2. The trends for the Reynolds numbers in both liquid metal electrodes vary stronger and cannot be fitted by a power law. The reason for this behavior can most probably be attributed to the more complex boundary conditions at the interfaces which couple the motion between separator and electrode.

The scaling of the Nusselt number – the measure for the global transport of heat – with respect to the Grashof number is displayed in Fig. 19. A power law exponent of $0.2$ as suggested for Rayleigh-Bénard convection with internal heat sources by Goluskin and van der Poel (2016) fits well to the present simulation data. RMC has generally a slightly higher Nusselt number in comparison to RC, which is a result of the additional Marangoni effects. The values for the Nusselt numbers are also listed in Tab. 5 in the appendix. Note also that the Marangoni number varies as $\rm{Ma}\propto G^{3/5}$ for the present cases of variable layer height. Figure 20 shows a vertical cut through the three-layer system of case H6. The temperature contours indicate a strong temperature gradient at the interface II. It can also be seen that the thermal boundary layers at both interfaces have different thicknesses.

The variation of the convection pattern with height is represented by the vertical gradient of the temperature $\partial_{z}T^{(2)}(z=d_{21})$ at the upper interface as seen in Fig. 21. We use this quantity because it displays a larger contrast between inflow and outflow regions in comparison to a plot of the interfacial temperature. At the bright temperature contour areas one observes inflow and at the darker temperature contour areas outflow of the fluid from the interface. With increasing Grashof number, the cellular pattern is fragmented into ever finer substructures such that a larger amount of heat can be carried across the layer.

IV.3 Joint variation of the middle-layer size and the electrical current density

In the present subsection, we will vary the relative thickness of the middle layer, $d_{21}$. It will be reduced to 10mm, which corresponds to $d_{21}=0.5$ in series JM, to 6mm ($d_{21}=0.3$) in series JS, and to 2mm ($d_{21}=0.1$) for series JXS. The outer electrode heights will remain unchanged at a value of 20mm. Moreover, for each of these three geometries, different current densities $j_{0}$ were taken as summarized in Tab. 4.

Let us start with the medium sized middle layer at $d_{21}=0.5$. All four cases JM1 - JM4 are unstable to Rayleigh-Bénard convection which is readily predicted by a comparison of $G$ with the corresponding critical Grashof numbers $\tilde{G}_{c}$ in Tab. 4. Again, we simulated the cases of RMC as listed in the table together with corresponding RC simulations with $Ma=0$ in order to emphasize the effects of the interfacial tension better.

With a further reduction of the thickness of the middle layer $\Omega^{(2)}$ the corresponding Rayleigh-Bénard and Marangoni modes are damped to the subcritical range as they rely on convection in this layer. Thus a higher current density of $j_{0}$=5kA/$m^{2}$ and consequently a higher heating rate are required to sustain convection (see also Tab. 4). The case JXS1 with the shallowest middle layer remains stable although the current density is increased to a value as high as $j_{0}=20$kA/m2. Only when the outer layer sizes are increased to $d^{(1)}=40$mm, the system becomes unstable again and shows the upper layer Rayleigh mode (case JXS2).

The observations from these three additional series of simulation runs are as follows. Figure 22 (top panel for JM and bottom panel for JS) shows the Reynolds numbers $\rm{Re}^{(i)}$ as a function of the Grashof number magnitude $|G|$. The data for both series follow again a similar power law trend as in the previous parameter study for H1 to H6 with an exponent that is close to 1/2. The Reynolds numbers $\rm{Re}^{(3)}$ for cases with the largest Grashof numbers (JM4 and JS5) show a significant enhancement. This large amplitude of the Reynolds number is caused by the strong Rayleigh-Bénard mode which is established in the upper liquid metal electrode. The corresponding Rayleigh number $\rm{Ra}^{(3)}$, which is defined in Eq. (66), is $\rm{Ra}^{(3)}=-1092$ for JM4 and $\rm{Ra}^{(3)}=-2673$ for JS5.

Figure 23 shows the Nusselt numbers for simulation series JM and JS. It is seen that $\rm{Nu}$ is systematically larger for the shallower middle layer at $d_{21}=0.3$ when plotted against the Grashof number. Also, RMC has a consistently higher value of $\rm{Nu}$ than RC. An approximately universal scaling relation for the turbulent heat transfer, for which all data points collapse to a single line in a double-log plot, can be obtained when the original Grashof number is substituted by the relative Grashof number $G_{r}$. It is defined as the actual Grashof number divided by the corresponding critical value $G_{r}=G/\tilde{G}_{c}$. This relation is shown in the bottom panel of Fig. 23 for all simulation cases in series JM and JS. The data are fitted well again by the power law of $\rm{Nu}\propto G^{0.2}$ that was discussed already in section IV B. This particular scaling of the turbulent heat transfer seems thus to be very robust in our parametric studies.

V Different convection regimes

As the full nonlinear simulations in the present three-layer model showed, there are four different modes for thermal convection that we have already discussed, and an additional one not discussed so far. These are (1) the Rayleigh-Bénard mode due to internal heating in the middle layer, (2) the standard Rayleigh-Bénard mode in the upper layer due to the hot interface II and the cold top, (3) the Marangoni mode at the upper interface II, (4) the Marangoni mode at the lower interface I, and (5) an anti-convection mode in the middle layer and bottom electrode that has not been discussed so far. On the one hand, these modes are partly coupled to each other, one the other hand, they are partly insignificant for particular parameter sets. They contribute differently to the global transfer of heat and momentum. The following section will discuss and summarize these different mechanisms and how their appearance can be predicted.

V.1 Rayleigh-Bénard mode in molten salt layer

The primary mechanism for convection in the investigated system is the unstable temperature distribution between the maximum temperature close to $z=d_{21}/2$ and the interfacial temperature at the upper interface II, causing buoyant convection. In a classical study, Sparrow et al. (1964) investigated the stability of a single fluid layer under internal heating, with no-slip and isothermal boundary conditions at the top and bottom. These results of Sparrow et al. (1964) for one layer can be translated to the present three-layer system. Aside from the formal linear stability of a layer with no-slip and isothermal boundaries, they found that the threshold of convection can be described by a single Rayleigh number criterion if the interface temperatures do not differ too much. The relevant Rayleigh number $\rm{Ra}^{(2)}$ relates the difference of the temperature maximum and the upper interfacial temperature and their distance. Formally one can thus define the Rayleigh number as follows

[TABLE]

where the distance of the temperature maximum to the interface is expressed with the help of (40). In terms of this Rayleigh number (72), Sparrow et al. (1964) showed that the threshold for the linear instability falls in a range of $\rm{Ra}^{(2)}=560$ – 595 (depending on the outer temperature difference which they denoted by a parameter $N_{s}$). Later, Char and Chiang (1994) studied the linear stability of a single layer subject to internal heating with different boundary conditions for the linearized Boussinesq equations. The authors derived a critical Rayleigh number for a free-slip, isothermal and free-slip, thermally insulated upper boundary, respectively. They derived a critical Rayleigh number of $\rm{Ra}^{(2)}$=296.20 (112.64) for the isothermal (thermally insulated) case. Which boundary conditions are better suited to account for layer (3) is unclear. A hint for the best thermal boundary condition might be given by the ratio of thermal conductivities $\lambda_{32}=137.54$. This ratio suggests an isothermal rather than an thermally insulating boundary condition.

To check the relevance of this newly defined Rayleigh number $\rm{Ra}^{(2)}$, the calculated critical Grashof numbers can be transformed to critical Rayleigh number for a fixed set of material parameters and geometry. For the four different geometries which are given in Tab. 4, one obtains a critical Rayleigh number of $\tilde{\rm}{Ra}^{(2)}_{c}=-290.7$ for $d_{21}=1$, $\tilde{\rm}{Ra}^{(2)}_{c}=-281.7$ for $d_{21}=0.5$, $\tilde{\rm}{Ra}^{(2)}_{c}=-276.9$ for $d_{21}=0.3$ and $\tilde{\rm}{Ra}^{(2)}_{c}=-92.7$ for $d_{21}=0.1$. Obviously, for the thick middle layer the critical Rayleigh number fits thus well to the prediction of Char and Chiang (1994) of an isothermal upper interface. When the middle layer thickness is reduced, the upper-layer Rayleigh-Bénard mode gets more involved and reduces the stability threshold. For a shallow middle layer, the internal convection mode gets insignificant. At this point we recall that the separator layer should be thin to keep the Ohmic losses at a minimum

Finally, let us note that in the present case, the number of parameters in the marginal stability analysis, which can be varied independently of each other in the nondimensional equations, is less than 18. The stability problem in Sec. III.2 can be rescaled to a new velocity $w^{+}=w\Pr$. By doing this, one observes that only the product $\rm{Ra}=G\,\Pr$ appears in the linearized equations for perturbations with a growth rate of zero. Furthermore, the conductivities of the outer layer can be regarded as infinite as is done in the solution for pure conduction. Overall, this results in a reduction to 15 essential parameters that govern linear stability threshold.

V.2 Rayleigh-Bénard mode in upper liquid metal electrode

The top electrode has a linear temperature profile, and thus its stability to buoyant convection can be estimated by the Rayleigh number $\rm{Ra}^{(3)}$ which was introduced in Eq. (66). For the simulated geometries, this Rayleigh-Bénard mode was the primary source for convection in the series of the very shallow cases, JXS. According to our stability calculation the critical Rayleigh number for $d_{21}=0.1$ is $\rm{Ra}^{(3)}_{c}=-1290.4$. For thicker middle layers the convection mode is observed for higher electrical current densities only. The Rayleigh number for which convection is observed is $\rm{Ra}^{(3)}=-2673$ for JS5, $-1092$ for JM4 and $-6192$ for H6. In summary, if shallow middle layers are combined with thicker layers for the upper electrode, this specific Rayleigh-Bénard mode can become the primary source for convection in the three-layer system. If the internal convection mode in the molten salt layer is active additionally, this Rayleigh-Bénard mode leads to a considerable increase of the velocity fluctuations in the top layer. This parameter regime would be an example where both RC flows dominate the global transport.

V.3 Marangoni modes at the liquid metal-molten salt interface

The general requirement for the existence of stationary MC can be inferred from the stability theory of two layers by Sternling and Scriven (1959). The respective interface is prone to MC if heat is transferred from the phase with lower thermal diffusivity into the phase with higher one, given that interfacial tension decreases with temperature. Formally, the requirement for the respective interfaces in the present system is that $\kappa^{(2)}<\kappa^{(1)}$ and $\kappa^{(2)}<\kappa^{(3)}$. For the present liquid metal battery model both inequalities hold. Also, we can expect that this ratio of diffusivities is typical for the compound system of a liquid metal and a molten salt in general since the diffusivity of the molten salt is by two orders of magnitude lower than that of the liquid metals in the electrodes.

In some of our cases, linearly unstable MC modes could grow, but the Rayleigh-Bénard mode due to internal heat sources was usually amplified stronger. We demonstrated that the main reason of the appearance of MC, as a part of the RMC regime, is a reduction of the length scales of the cellular flow patterns which is in line with an increase of the Nusselt number.

The stability threshold of MC depends crucially on the ratio of transport coefficients of the adjacent phases. Thus, a threshold based on simple parameters, as in the case of the RC modes in the middle and upper layers, respectively, cannot be provided. It might be possible that MC is the primary mechanism for other sets of material parameters. This is stated in view of the missing experimental data for the interfacial tension of liquid metal–molten salt interfaces as a function of the temperature. What can be certainly stated here is that MC is damped by a reduction of the middle layer thickness (cf. Fig. 7), by a reduction of the temperature difference and by an increased viscous friction, respectively. Thus a relevant Marangoni number should remain proportional to $(d^{(2)})^{3}$ as suggested by Fig. 7(a). The scaling with the third power of thickness is due to our choice of temperature scale. Finally, we note that a closed expression for the stability threshold of pure MC with the generalized temperature profile of Eqns. (35)–(37) might be derived similar to the case of linear temperature profiles (Géoris et al., 1993).

V.4 Anti-convection mode in the lower liquid metal alloy electrode

The fifth source of convection concerns the bottom and the middle layer and was not discussed so far. Owing to the low Prandtl number $\Pr^{(1)}$=0.012 in the bottom layer, the diffusive transport of the temperature field dominates in this layer over advection for almost all simulated cases, which is readily estimated by the Péclet number $\rm{Pe}=\rm{Re}^{(1)}\,\Pr^{(1)}<1$. This circumstance and the rather low heat capacity of the Pb-Bi layer will generate a diffusive transport into the bottom layer that is initiated at the localized hotspots due RC mode in the middle layer. This can be seen in Fig. 20 around $x=2$. The resulting horizontal temperature gradients cause a flow where denser colder fluid tends to displace hotter fluid. And indeed, the work that is done by gravity on the lower layer and which is proportional to $\langle-\beta_{T}^{(1)}T^{(1)}u_{z}^{(1)}\rangle_{xyz}$ is found to be positive for all simulated cases.

Anti-convection has been classically described for two layers heated from above (Welander, 1964; Boeck et al., 2008; Merkt and Bestehorn, 2012). In this regime, which is noted as "stably stratified", buoyant convection is nevertheless possible when the fluid mass density decreases with temperature, as it is usually the case. The requirement for the anti-convection regime is that material parameters of both layers are in a certain proportion to each other as discussed by Welander (1964). His prediction for a system with two layers heated from above can be transferred to the lower and middle layers in our case. We found that anti-convection would be possible if the thermal expansion coefficient of the middle layer is smaller than the one used for simulations. For a value of $\beta_{T,21}=0.25$ instead of 2.59 anti-convection should become possible.

To demonstrate that this mode can indeed play a role in our case of internal heating in the middle layer, we performed a simulation with the parameters of run JM4 but suppress all sources for convection, i.e., we set $\rm{Ma}=0$ as well as $\beta_{21}=\beta_{31}=0$. The only physical mechanism which is kept is the change of fluid density and the resulting buoyancy effect in the lower layer. Figure 24 demonstrates that indeed – though the lower layer is stably stratified – convective fluid motion sets in. In the middle layer, temperature perturbations are generated which heat the interface and produce an upwelling flow in the bottom layer. However, in the real system, the buoyancy in the middle layer counteracts while the Marangoni effect concurs with anti-convection since the interfacial flow is from the hot to cold regions. For the chosen parameters, the present system develops a state of chaotic pulsating and translating convection cells. The time-averaged flow amplitude is considerable. One gets values of $u_{rms}^{(1)}=56.8$, $u_{rms}^{(2)}=36.97$ and $u_{rms}^{(3)}=9.86$, respectively.

VI Interfacial deformation

The present model neglects effects of the deformation of the interfaces. We, therefore, need to verify the consistency of if this assumption. To this end, we estimate the flow-induced deformation of the interface from the equilibrium shape.

The shape of the interface is essentially governed by the balance of normal stresses at the interface. This condition was approximated by assuming zero vertical velocity. For a material interface, the relation by Edwards et al. (1991) holds at any point on the interface. It is given by

[TABLE]

where all quantities carry a physical dimension here. Three new quantities have to be defined. First, the interface normal $\mathbf{n}$ is pointing from phase (2) to phase (3) which is required by sign convention. This surface normal can be determined from a height function $h(x,y,t)$ that represents the distance of the interface from the unperturbed level, $z=0$, namely

[TABLE]

Second, the mean curvature $H$ encodes the deformation of the interface. It can be calculated by means of the height function $h(x,y,t)$ and is given by

[TABLE]

where $h_{x}$ and $h_{y}$ denote partial derivatives of $h$ with respect to $x$ and $y$.

Third, the rate of strain tensor in the respective layer, $\hat{\mathbf{E}}^{(i)}$, is given by

[TABLE]

where $\tilde{\mathbf{u}}^{(i)}=\mathbf{u}^{(i)}U_{vis}$. The rate of strain tensor accounts for the normal viscous stresses at the interface. Without loss of generality, we restrict the discussion to the upper interface since velocities have a higher amplitude there and density stratification is less significant. Further calculations are done using dimensional quantities.

In the present framework, we assumed a planar interface, i.e., $h=0$ results in $H=0$, which is appropriate if the interface is pinned at the edges and if the right hand side of Eq. (73) remains small. Small should be understood in the sense that only deformations $\delta h$ of the interface are considered that are small compared to any length scale of the flow. A formal way to estimate the deformation is to perform a perturbation expansion of Eq. (73) with respect to $h$ as done in Nepomnyashchy et al. (2006) (see pp. 13 therein). This requires the calculation of the pressure from the zero-order solution of the Navier-Stokes equations.

Here, we will take a simpler approach by estimating a typical interface deformation as follows. Let us consider a fluid parcel that moves downwards into the bulk of the lower layer away from the interface with velocity $U$. Simultaneously a column of denser liquid is present due to a temperature reduction by a magnitude of $\delta T$ in the top layer at this position. This causes a pressure difference which can be estimated by

[TABLE]

where $L_{V}$ is the typical vertical variation scale. This process forces the interface to bulge downwards by $\delta h$ thus causing a negative curvature $\delta H$. Simultaneously, the "hydrostatic contribution" to the pressure at the interface would rise by $g\delta h(\rho^{(2)}-\rho^{(3)})$. Now, one can relate these contributions via Eq. (73) and use the scaling relation for the viscous contribution, which is given by Eq. (76), and the curvature from Eq.(75) to obtain

[TABLE]

where the curvature term Eq. (75) is estimated with the help of a horizontal length scale $L_{H}$ and the viscous contribution as well as the density contribution by means of a vertical length scale $L_{V}$. Finally, the interfacial deformation is estimated by rearranging (78). One obtains

[TABLE]

In order to estimate the deformation, we will adopt typical values of case H6. The velocity in the middle layer is $U\sim 3$mm/s (see Tab. 5), the length scales $L_{V}\sim L_{H}\sim 10$mm as seen in Fig. 21. A typical temperature perturbation is $T^{II}\Theta\approx$ 1K. We note that the length unit for H6 is 40 mm. Inserting those values in Eq. (79) yields a $\delta h\approx$ -1.22 $\mu$m which is considerably smaller than the smallest flow structures. In this view, the plane interface approximation is indeed appropriate and its errors should be insignificant compared to the other physical processes which have been neglected. For instance, vibrations acting on the system could trigger gravity waves.

Finally, we employ the estimate of the interface deformation by Herreman et al. (2015) (see their Eq. (5.7)) that equates the simulated kinetic energy with an increase in potential energy. This leads in our case to an estimate of

[TABLE]

Thus we can conclude that no relevant dynamic interface deformation is triggered by thermal convection. The density difference between the layers is strong enough that any deformation is limited to a tiny magnitude. Furthermore, for an LMB application, the size of the middle layer for H6 is comparatively large and would cause a high voltage drop. Thus no higher velocities from the internal convection modes are to be expected. However, in systems with a tall upper layer and a strongly amplified RC mode, the deformation might become considerable again. In this case, we do not have reliable estimates of the velocity. Similar conclusions have been drawn by (Zikanov, 2015) for thermal convection and Herreman et al. (2015) for magnetohydrodynamic effects.

A curved equilibrium shape, which appears if the contact angle differs from 90 degrees, can have a particular influence on thermal convection. Such wall effects should be studied together with the particular vessel geometry since the thermal properties of the battery vessel may induce horizontal temperature gradients, which in turn will induce Marangoni convection.

VII Summary and Conclusions

We have numerically studied thermal convection, induced by a uniform resistive heating, in a system of three liquid layers which serve as a simplified model for a liquid metal battery. Our model comprised the flow induced by interfacial tension gradients and gradients in density both of which are due to temperature gradients. The present LMB model has been simplified by disregarding electrochemical effects, the transport of mass across the separator layer, and the interaction with a vessel wall, respectively. Furthermore we exclude magnetohydrodynamic effects in connection with the current density in the cell. Our main conclusions can be summarized as follows.

(1) In the three-layer system Li||LiCl–KCl||Pb–Bi, driven by resistive heating, four main modes can drive convection. A fifth mode of anti-convection can appear as a secondary effect. The primary source for convection is the Rayleigh-Bénard convection mode due to internal heating in the molten salt layer. It obeys similar properties as a single convection layer subject to internal heating Goluskin and van der Poel (2016). This process acts together with the second potential source of convective motion, the Marangoni convection effects at the upper interface. The typical thermal diffusivity ratios between a molten salt and a liquid metal lead to the appearance of Marangoni convection at the upper as well as at the lower interface of the separator to the electrodes since heat is transported from the middle layer with significantly lower thermal diffusivity to the electrodes. The fourth source for convection is the unstable stratification in the upper electrode which causes a turbulent convection flow that is comparable to classical Rayleigh-Bénard convection. Finally, anti-convection can increase the convective motion in the lower layer. For a fixed temperature drop, all modes will be damped when the thickness of the molten-salt layer is decreased, except the Rayleigh-Bénard convection mode in the top electrode.

(2) The exact linear stability thresholds for the onset of convection have been calculated for different system geometries and electrical current densities. By comparison with a full nonlinear DNS, it was found that those linear stability predictions indeed describe the conditions for the onset of convection well. The most critical mode in almost all of the simulated cases is the internal- heating Rayleigh-Bénard mode. The stability threshold of this mode can be estimated by the middle layer Rayleigh number $\rm{Ra}^{(2)}$. Only for a particularly shallow middle layer, specifically for $d^{(2)}/d^{(1)}=0.1$ in our study, the classical Rayleigh-Bénard mode in the upper layer can become the main source for the onset of convection. The threshold for this mode can be determined by the Rayleigh number $\rm{Ra}^{(3)}$.

(3) Several of our simulated configurations are also unstable with respect to Marangoni convection. Although this process is not the dominant mechanism of instability, it contributes significantly when both physical effects – Rayleigh-Bénard and Marangoni effects – are jointly enabled. The characteristic scale of the cellular flow patterns changes to a smaller length scale and the Nusselt number is increased simultaneously. The impact of the Marangoni effect on flow velocities is twofold. On the one hand, it is found that MC leads to a decrease of the velocity in the upper layer once the upper RC mode is active. On the other hand, for shallower middle layers, MC leads to a considerable increase of flow velocity in the lower and middle layers, respectively. To conclude, although thermal MC turns out not to be the main source for convection, it affects the flow structures considerably. Further studies of convection in LMB should to our view include MC effects.

We also expect that a number of changes will result from the presence of lateral walls, in particular when the aspect ratio of the cell (diameter/height) is small. Small aspect ratios will affect the cellular pattern formation and can thus suppress the heat transfer. This is known from RC Bailon-Cuba et al. (2010). A finite thermal conductivity of the side walls will amplify the heat transfer, e.g. by a more pronounced thermal plume formation at the side walls Verzicco (2002). Resulting horizontal temperature gradients can also drive buoyant as well as Marangoni convection. Finally, the wetting behavior of the liquid phase governs the interface geometry near the vessel walls. Additional corner modes could thus be amplified in connection with the heat transport through the walls.

In future studies, the foundations of the physical model as well as the practical setup of prototypes will require significant advances. The coupling between electrochemical reactions and the related ion transport has to be combined with the thermal convection processes in order to receive a more precise and complete picture of the dynamical turbulent processes in LMBs. Furthermore, the assumption of a homogeneous current has to be improved since particular electrode geometries (Ning et al., 2015) are employed in LMB prototypes where the coupling between to the reaction rates of the electrochemical reactions and convection processes can lead to heterogeneous heating. Moreover, the gradients in the composition of the lower electrode could lead to additional gradients in the surface tension since lithium is transferred across the interfaces. This might be an additional source for a stable density stratification. As we have discussed, some unknown material parameters, such as the interfacial tension, or resulting uncertainties limit the predictability of the importance of different convection processes. They ask for further experimental studies in the near future.

Acknowledgements.

Financial support by the Deutsche Forschungsgemeinschaft in the framework of the Priority Programme SPP1506 is gratefully acknowledged. Furthermore, we thank the Computing Centre (UniRZ) of TU Ilmenau for access to its parallel computing resources. We also thank Oleg Zikanov and Norbert Weber for helpful discussions.

Appendix A Root mean squared velocity and Nusselt numbers

In this appendix, we list the root mean square velocities and Nusselt numbers of the series H, JM and JS, respectively. They have been discussed in section IV. We exclude those runs which were linearly stable as well as the series JXS.