MHD stability of large scale liquid metal batteries
MHD stability of large scale liquid metal batteries
Email address for correspondence: V.Bojarevics@gre.ac.ukK. Pericleous
The aim of this paper is to develop a stability theory and a numerical model for the three density-stratified electrically conductive liquid layers. Using regular perturbation methods to reduce the full 3d problem to the shallow layer model, the coupled wave and electric current equations are derived. The problem set-up allows the weakly non-linear velocity field action and an arbitrary vertical magnetic field. Further linearisation of the coupled equations is used for the linear stability analysis in the case of uniform vertical magnetic field. New analytical stability criteria accounting for the viscous damping are derived for particular cases of practical interest and compared to the numerical solutions for variety of materials used in the batteries. The new criteria are equally applicable to the aluminium electrolysis cell MHD stability estimates.
Key words: Interfacial waves, MHD interaction, linear stability, liquid metal battery, shallow layer approximation, aluminium electrolysis cell
Liquid Metal Batteries (LMBs) have many important characteristics for efficient practical use in combination with renewable energy sources on a national energy grid scale in the future. A relatively high voltage efficiency at high current densities of this storage technology are due to liquid-liquid electrode-electrolyte interfaces that enable high speed charge transfer, high total current capability, low ohmic losses, as well as rapid mass transport of reactants and products to and from the electrode-electrolyte interfaces by means of liquid-state diffusion (Kim et al.2013).
The liquid state of the main components necessitates consideration of the fluid dynamics in LMBs. A number of recent publications are devoted to the problem, e.g. a possibility of the Tayler instability (Weber et al.2014; Herreman et al.2015), the thermal convection (Shen & Zikanov 2016), observation of vortical flow in a LMB model (Kelley & Sadoway 2014), and a simplified model of sloshing in a three layer system (Zikanov 2015). The main motivation of these investigations is to prevent the possibility of direct contact between the molten metal anode and cathode that may occur due to electro-magnetically driven destabilizing interface motion. On the other hand, the controlled mixing enhances mass transport improving the cell performance, preventing the accumulation of intermetallic compounds at the electrode-electrolyte interface.
LMBs are thought to be easily scalable on the cell level due to their simple construction using the natural density stratification of the liquid layers. Large cells of several cubic meters total volume have a potential to operate at very high power value (Bojarevics & Tucs 2017). High current densities coupled to the magnetic field (created by the currents in the cell, the supply bars and the neighbour cells) lead to significant electromagnetic forces. Such forces in stratified liquid layers with large surface areas may cause a long wave interfacial instability as it is well known in the case of Hall-Heroult cells (HHC), as first described by Sele (1977). The manifestation of this instability in LMB is the subject of this paper.
In a typical HHC the electric current, of total magnitude 150−800 kA, enters the cell from the carbon anodes, passes through the liquid electrolyte and aluminium layer, and exits via the carbon cathode blocks at the bottom of the cell. The liquid layers are relatively shallow, 4−30 cm in depth vs 4−20 m in horizontal dimension. The small relative depth of the layers and the small difference of the liquid densities facilitates the instability development.
The ratio of electrical conductivities of the cell materials is another significant parameter. The liquid metal is a better conductor (∼106 S/m) than the carbon (∼104 S/m), while the electrolyte is about two orders of magnitude less conductive (∼102 S/m). The significantly lower conductivity of the electrolyte means that this layer is responsible for the majority of electrical losses in the cell. Joule heating is necessary to heat the cell and to keep the metal liquid, however the total voltage drop must be as low as possible in order to achieve a better electrical efficiency. A small perturbation of the interface between liquid layers may cause a substantial redistribution of the current in the cell.
First attempts to explain the interfacial instabilities were made by Sele (1977), Urata (1985), Sneyd (1985) and Moreau & Ziegler (1986). A more involved understanding of the physical mechanism was provided by Sneyd & Wang (1994), Bojarevics & Romerio (1994), and Davidson & Lindsay (1998). The mechanism is based on the standing gravity wave modification due to the electric current redistribution. The electric current density in the electrolyte increases above the wave crests, resulting in a high density horizontal current in the shallow liquid metal layer. In the presence of a vertical magnetic field the electromagnetic force excites another standing wave mode orthogonal to the initial perturbation. The new wave mode is coupled to the original mode, and the oscillation frequency is shifted. The frequency shift increases with the rise of the magnetic field until at a critical value, when the two wave frequencies coincide, an exponential growth of the amplitude indicates the onset of instability. In general, the above process is described by the following set of equations:
where H\boldmathk is a vector which represents the amplitudes of the original gravitational modes \boldmathk\boldmathk=(kx,ky), \mathsfbiω2\boldmathk is the matrix of the gravitational frequencies, \mathsfbiG\boldmathk\boldmathk′ is the interaction matrix, E is the dimensionless parameter characterizing the electromagnetic forces. The mode coupling is included in \mathsfbiG\boldmathk\boldmathk′, where each column represents the Lorentz force (Fourier decomposed) in response to the gravitational wave modes. These coupled equations represent an eigenvalue problem for the square of the new complex frequencies μ (H\boldmathk∼eμt). The matrix \mathsfbiG\boldmathk\boldmathk′ is real anti-symmetric, and in a general case the eigenvalues are shifted increasing the magnetic field (Bojarevics & Romerio 1994). Onset of the instability starts at a critical value of E at which the exponentially growing part of the complex eigenvalue μ appears. See also the paper by Antille & von Kaenel (2002) where the numerical eigenvalue solution is analysed when approaching the instability threshold. A key point noted in the papers, is that the dominant contribution to the perturbed Lorentz force arises from the interaction between a horizontal current in the aluminium layer and the vertical component of the background magnetic field.
Davidson & Lindsay (1998) derived a simple mechanical analogue which captures the basic features of the instability. The liquid aluminium layer is represented by a compound pendulum that consists of a large flat aluminium plate attached to a top surface by a light, rigid strut. The strut is pivoted at its top end so that the plate is free to swing along two horizontal axes x and y. The fluid system of infinite motion freedom is reduced to only two degrees of freedom. Zikanov (2015) constructed a similar mechanical model for instability description in the LMB taking into account an additional top liquid metal layer. The metal layers of the battery are represented by solid metal slabs rigidly attached to weightless rigid struts pivoted at the top. The free oscillations of the slabs imitate the sloshing motion of the liquid layers. The slabs are separated from each other by a layer of a poorly conducting electrolyte. Two destabilization mechanisms were considered: 1) interaction of a purely vertical magnetic field and horizontal currents, similar to HHC, 2) interaction between the current perturbations and the azimuthal self-magnetic field from the total vertical current. The first mechanism will occur in real batteries if a sufficiently strong vertical magnetic field is present due to the presence of external current supply. The batteries of a square or a circular horizontal cross section will be always unstable if even a small field is present. The second mechanism appears to be more challenging since the azimuthal magnetic field, unlike the vertical magnetic field, cannot be reduced via optimization of the current supply lines unless they cross the liquid layer (Weber et al.2014). The existence of the second instability type was predicted by Munger & Vincent (2008) for HHC case, yet needs more clarification for the LMB case. The approach developed by Davidson & Lindsay (1998) and Zikanov (2015) is purely mechanical. However, the principal physical mechanism could be valid, due to the fact that sloshing motions generated in the shallow liquid layers are inherently large scale, and so their qualitative behaviour can be approximately described using the coupled pair of long wave modes approach.
More realistic fluid dynamic description can be achieved starting from the full set of Navier-Stokes equations by means of the shallow layer approximation and a systematic derivation of a set of coupled wave equations governing the three fluid layers. The hydrodynamic coupling is realised by pressure continuity at the common interfaces. The continuity of the electric potential and the supplied electric current will introduce the electromagnetic coupling of waves. In the following sections the linear stability of coupled modes will be investigated in the presence of a purely vertical magnetic field, accounting for the continuous electric current in the LMB model. The role of dissipation rate will be analysed using both analytical tools and numerical solutions. Analytical criteria for the cell stability will be established using approximations suitable for a practical cell design, including the solid bottom and top friction effects on the shallow layers.
2 Interfacial dynamics
Hydrodynamics of the three density stratified electrically conductive liquid layer system, schematically represented in figure 1, in the presence of electro-magnetic fields, are described by the following equations
where the indices i,j=1,2,3 correspond to the coordinates (x,y,z), the velocity components are given as (u1,u2,u3), the summation over repeated indices is implied, ρ represents density, ν - effective viscosity, g - the gravitational acceleration, p - pressure and the vector of Lorentz force \boldmathf\boldmathf=\boldmathj\boldmathj×\boldmathB is computed in each of the 3 layers, j is the current density and B the magnetic field. In this paper the horizontal dimensions of the cell are assumed to be much larger compared to the vertical depth, so that the description can be based on a systematically derived shallow layer approximation. The velocity components in each layer can be represented as an expansion in a small aspect ratio parameter δ=maxh/minL, where h is a typical depth, for instance the unperturbed metal layer, and L is the characteristic horizontal dimension (width of the cell):
where a stretched vertical coordinate ¯¯¯z=z/δ is introduced. The u3 expansion starts with the δ-order due to (2). If all three components of the electromagnetic force density are of the same order of magnitude: fx∼fy∼fz, and the horizontal pressure gradient components are of the same order as the corresponding force components: ∂ip≈fi, then the vertical component of the gradient ∂3p∼−ρg≫fz. According to these estimates, the leading horizontal (i=1,2) components of (2) are
The vertical component of the equation (2) gives the leading order terms as:
The hydrostatic pressure in the liquid layers adjacent to the interface H1(x,y,t), see the figure 1, can be expressed by
where the index m=1,2 stands for the layer number and p(1)p is the reference pressure at the moving interface z=H1(x,y,t). Similarly the same pressure can be referenced to the interface H2(x,y,t), the corresponding pressure given by
The last equation can be linearised if an additional approximation of a small wave amplitude is introduced: hk(x,y,t)=h0k+εh′k(x,y,t) for the layer thickness or equivalently Hm(x,y,t)=H0m+εζm(x,y,t) for the interface position, where the additional small parameter ε=maxA/h is introduced. A is a typical wave amplitude and h0k, H0m are the unperturbed values, ζm are the interfacial perturbations. For each particular layer the depth average horizontal velocity divergence can be expressed as:
The volume conservation requires:
In order to estimate the leading order of terms in the depth averaged momentum equations (2) dimensionless variables of order O(1) are introduced using the following scaling: L - for the coordinates x,y; h for a typical layer thickness; ε√gh - for the wave velocity; L/√gh - for the time; ρkgh for the pressure (k=1,2,3). For typical geometries considered in this paper δ=h/L≈ 0.2 m/8 m = 0.0250≪1, whereas, ε=A/h≈ 0.005 m/0.2 m = 0.0250≪1.
The horizontal pressure gradient from the expressions (2) and (2), neglecting terms of the δ and higher order, can be substituted in the depth averaged, nondimensionalized horizontal momentum equation (2). For the layers adjacent to the lower interface H1(x,y,t) the respective momentum equations are
The equations (2.0) and (2.0) formally give the connection between the reference pressures p(1)p and p(2)p defined on the two interfaces, however being valid in the same fluid layer k=2 (electrolyte). The alternative representations are required for the wave equation derivation. After the integration over depth the dissipative terms in (2) are replaced by empirical expressions used for the shallow layer approximation (Moreau & Evans 1984; Rodi 1987) using a linear in velocity friction law with the coefficients kfk. The electromagnetic interaction parameter (the ratio of electromagnetic force to the gravity force perturbation) is introduced as Ek=IB0/(L2ρkgεδ), where I is the total electric current, B0 is a typical magnitude of magnetic field. The corresponding magnitude of E can be estimated, using typical values for I=105 A, B0=10−3 T, L=8 m (width of cell), ρ=1.6×103 kg m−3 (liquid magnesium for the top metal), g=9.8 m s−2, ε=δ=0.025: E=0.32=O(1). The electromagnetic term is of the same order of magnitude as the leading terms, while the nonlinear wave motion terms are of lower order (∼ε) and will be neglected later in the linear theory, but retained for a numerical solution.
The wave equations for the coupled interfaces can be derived following the procedure described in Bojarevics (1992):
take time derivative of the non-dimensional linearised equations (2.0),(2.0):
substitute (2.0), (2.0) into the horizontal divergence of (2.0),(2.0),
take the difference of the resulting equations.
This procedure eliminates the common unknown pressure p(1)p on the interface ζ1:
The corresponding boundary conditions for the normal velocity un=0 at the side walls can be obtained by taking the difference between (2.0) and (2.0) to eliminate the common pressure at the interface ζ1:
In a similar manner the wave equation for the upper interface ζ2 can be obtained:
and the corresponding boundary conditions are
The new constants introduced in the above equations are defined as:
Note, that in (2) and (2) summation over repeated indices is limited to the two horizontal dimensions. As it can be seen from (2) and (2), both interfaces can not be considered independently due to the presence of coupling terms. In the following we will assume that the interfacial friction kf2 is negligible. The set of the equations in the electrically nonconductive limit and in the absence of viscous dissipation is in correspondence with the set of the equations obtained by Robino et al. (2001) derived for the dynamics of internal solitary waves in stratified 3 layer ocean. The aluminium electrolysis cell MHD wave model can be recovered if ζ2=0 in (2).
3 Electric current flow
For energy storage and supply the LMBs must operate in two regimes: charge and discharge, resulting in the current flowing (upwards or downwards). In this paper only the charging process is considered due to the physical symmetry of both operational regimes. The current flow in the layered structure is illustrated in figure 1. Similarly to Davidson & Lindsay (1998), we assume that the characteristic time-scale for the wave motion is much larger than the diffusion time of the magnetic field to satisfy the low magnetic Reynolds number approximation, leading to Rm=μσUh≪1, where μ is the magnetic permeability, σ the electrical conductivity, U is a typical velocity. Using typical values for the liquid metal μ=4π×10−7 H m−1, σ=3.65×106 S m−1, U=0.01 m s−1, h=0.1 m the estimated value of Rm≈0.004. A similar estimate can be obtained for the wave motion magnetically induced electric current ration to the basic current density: σUB/(I/L2)∼σε√ghB0/(I/L2)∼0.001.
In the low Rm approximation, when the flow effect is neglected, the electric current can be expressed as
and is described by a set of coupled Laplace equations for the electric potential φk(x,y,z):
where k=1,2,3 corresponds to the layer number. The continuity conditions for the electric potential and the normal current component \boldmathj\boldmathj⋅\boldmathn at the interfaces z=Hm(m=1,2) are
The normal derivatives at the deformed interfaces are defined as (assuming the summation over the repeated index i only):
With (3.0) the current continuity (3.0) at the interfaces H1 and H2 can be written explicitly in the nondimensional form in order to estimate the leading order terms:
where the four orders of magnitude difference in the electrical conductivities permit to define σ2/σ1=s1δ2, σ2/σ3=s3δ2 and the stretched ¯¯¯ζi=ζi/δ. These definitions allow us to compare numerically the electrical conductivities in the poorly conducting electrolyte relative to the well conducting liquid metals, and the effect of the small depth (∼δ) of the layers. The side walls of the domain are considered to be electrically insulating:
In this paper we assume, that the applied current distributions at the top and the bottom are uniform, equal and the interfacial perturbation do not influence the electric current distributions in the collectors:
where ¯¯¯¯¯Hk=Hk/δ. In principle, j(x,y,t) could be used, however requiring an external circuit solution.
The set of Laplace equations (3.0) can be rewritten in a nondimensional form
The shallow layer approximation requires that the potential is expanded in terms of the parameter δ:
where the expansion terms are expressed in a similar manner as in Bojarevics & Romerio (1994):
where a, b, c, d, e, g are the coordinate x and y dependent functions that correspond to the unperturbed interfaces. The functions A, B, C, D, E, G are x, y and time t dependent, corresponding to the perturbed interfaces. Taking into account the previously described boundary conditions and neglecting the higher order terms, the unknown coefficients can be determined as shown in the Appendix A by equalising the similar order of magnitude terms.
Finally, the resulting set of the equations governing the electric current distribution in the system is obtained by introducing the perturbed potentials in both metal layers: Φ1=εB1, Φ3=εB3. For application in the LMB the dimensional equations for the electric potential perturbations are linearly correlated to the respective interface perturbations:
where h0k is the unperturbed layer thickness. The current distribution in the electrolyte is almost purely vertical due to fact that σ2≪σ1∼σ3 (similarly to Zikanov (2017)). According to (3.0) and (3.0), the current flow is perturbed by the electrolyte thickness perturbations. Finally, the corresponding dimensional current components can be expressed as
In the following section it will be shown that some of these perturbations may become electromagnetically coupled due to the presence of magnetic field and may lead to an instability resulting in a short circuit state at the extreme case.
4 Linear stability analysis
4.1 Coupled 3 layer problem
The wave equations (2)-(2.0) can be linearised neglecting the ε-order nonlinear terms and assuming a given magnetic field distribution, which can be expanded in terms of the small parameter δ: \boldmathB\boldmathB(x,y,z)=\boldmathB0(x,y)+δ\boldmathB1(x,y,z)+O(δ2). This expansion follows the same principle as the electric current representation \boldmathj\boldmathj(x,y,z)=\boldmathj0(x,y)+δ\boldmathj1(x,y,z)+O(δ2). At this stage we assume that the leading part of the magnetic field is caused by external sources (connectors, supply lines, neighbouring batteries etc.), therefore \boldmath∇×\boldmathB0=0 in the liquid zone. This allows us to neglect the EM force components resulting from the horizontal magnetic field interaction with the vertical unperturbed current: −j0zB0y%
\boldmathey due to the zero contribution to the horizontal divergence term in (2), (2.0): j0z(∂yB0x−∂xB0y)\boldmathez=0. The leading horizontal force components contain the horizontal electric current and the vertical magnetic field: j1yB0z%
\boldmathey, noting that j1z is δ-order lower than the horizontal perturbation current. This confirms to the assumptions made in the Introduction and implied in the previous studies on MHD stability of HHC, that the magnetic field is purely-vertical \boldmathB\boldmathB=\boldmathB\boldmathB0=B0z(x,y)\boldmathez, and it is caused by external sources.
The set of the wave equations for this case has the following form, after assuming that the friction at the electrolyte top and bottom is negligible in comparison to the friction at the solid top and bottom:
with the boundary conditions at the side-walls:
The electric potential distribution is governed by the set of equations (3.0), (3.0). Note, that the coefficients on the left hand side of (3.0) and (3.0) contain only the constant parts of the layer thickness.
The problem can be rewritten in a weak form by means of integrating the equations on the horizontal interface Γ∈[0,Lx;0,Ly] against a set of regular test functions (see the integral representation in Appendix B). The solution can be constructed in a Sobolev space W1,2(Γ), so that ζ and Φ are satisfying the corresponding equations for all test-functions ψ and q that belong to W1,2(Γ). The following set of functions is introduced:
where ˆζ\boldmathk(t), ˆΦ\boldmathk(t) are the spectral wave amplitudes and the perturbed potentials in Fourier space, whereas \boldmathk\boldmathk=(kx,ky). Note that the boundary conditions are satisfied in the weak sense only when using the functions (4.0)-(4.0). Taking into account the orthogonality properties of the cosine functions the set of wave equations including the boundary conditions can be rewritten (see the Appendix B) in the spectral coefficient space for k and \boldmathk\boldmathk′ mode interactions: