Modeling discontinuous potential distributions using the finite volume method, and application to liquid metal batteries

TL;DR

This paper introduces a finite volume method-based model for simulating three-dimensional electrical potential distributions in batteries with internal voltage jumps, validated on a liquid metal battery system.

Contribution

It develops a novel finite volume discretization that accurately captures internal voltage jumps at electrode-electrolyte interfaces in battery models.

Findings

Successfully validated with analytical test case

Simulated current distribution in a Li||Bi liquid metal battery

Analyzed discharge curve behavior

Abstract

The electrical potential in a battery jumps at each electrode-electrolyte interface. We present a model for computing three-dimensional current and potential distributions, which accounts for such internal voltage jumps. Within the framework of the finite volume method we discretize the Laplace and gradient operators such that they account for internal jump boundary conditions. After implementing a simple battery model in OpenFOAM we validate it using an analytical test case, and show its capabilities by simulating the current distribution and discharge curve of a Li||Bi liquid metal battery.

Tables

Table 1. Table 1: Material properties at 460 ∘ C [ 69 , 70 , 71 , 72 , 73 , 74 ] .

property	Li	Bi	salt	Cu	Ta	Ni
$\rho$ in kg/m3	490	9831	2690	2800
$\sigma$ in S/m	2.7$\cdot {10}^6$	7.2$\cdot {10}^5$	271	$58\cdot {10}^6$	$2.9\cdot {10}^6$	$3\cdot {10}^6$

Equations

$$\nabla \cdot \sigma \nabla \phi =0,$$

$$\frac{\partial }{\partial t}c=\nabla \cdot (D\nabla c),$$

$${\nabla }_{n}c=-\frac{\mathbf{J}}{zFD}\mathbf{n},$$

$$\nabla \cdot \sigma \nabla \phi =0,$$

$$\mathbf{J}=\sigma \nabla \phi ,$$

$$\nabla \cdot \sigma \nabla \phi =\sum _f{\sigma }_f\mathbf{S}(\nabla \phi {)}_f,$$

$$\nabla \cdot \sigma \nabla \phi =\sum _f\mathbf{S}{\sigma }_f(\nabla \phi {)}_f=\sum _f|\mathbf{S}|{\sigma }_f\frac{{\phi }_{\mathrm{N}}-{\phi }_{\mathrm{P}}+\Delta \phi }{|\mathbf{d}|},$$

$$A_d+A_o=S$$

$$A_o=\frac{{\sigma }_f|\mathbf{S}|}{|\mathbf{d}|},$$

$$A_d=-\sum _n\frac{{\sigma }_f|\mathbf{S}|}{|\mathbf{d}|}=-\sum _r{A}_{r,o},$$

$$S=\sum _f-\frac{{\sigma }_f|\mathbf{S}|\Delta \phi }{|\mathbf{d}|}.$$

$$\nabla \phi =\frac{1}{V}\sum _f\mathbf{S}{\phi }_f,$$

$${\phi }_{\mathrm{wN}}={\phi }_{\mathrm{wP}}+\Delta \phi ,$$

$${\sigma }_{\mathrm{N}}\nabla {\phi }_{\mathrm{N}}\cdot {\mathbf{n}}_{\mathrm{N}}={\sigma }_{\mathrm{P}}\nabla {\phi }_{\mathrm{P}}\cdot {\mathbf{n}}_{\mathrm{P}},$$

$${\phi }_{\mathrm{wP}}=f\cdot ({\phi }_{\mathrm{N}}-\Delta \phi )+(1-f)\cdot {\phi }_{\mathrm{P}}$$

$$f=\frac{{\delta }_{\mathrm{P}}\cdot {\sigma }_{\mathrm{N}}}{{\delta }_{\mathrm{N}}{\sigma }_{\mathrm{P}}+{\delta }_{\mathrm{P}}{\sigma }_{\mathrm{N}}}.$$

$$a$$

$$x_0$$

$${\sigma }_1$$

$${\phi }_1(x)=\alpha x+\beta ,$$

$${\phi }_2(x)=\gamma x+\epsilon .$$

$$\alpha$$

$$\beta$$

$$E/\mathrm{V}=\frac{0.19}{x+0.41}+0.5$$

Full text

††journal: Electrochimica Acta

1 Introduction

The objective of this article is twofold: we want to model both, the current distribution in, and charge-discharge curves of electrochemical cells. For good overviews on the subject, see [1, 2, 3]. As a first step, we will discuss only modeling of the time-dependent cell potential.

1.1 Battery voltage, overpotentials and polarization curve

Most battery models do not describe the distribution of current and potential. Instead, they provide only the total cell voltage, which can be simply calculated as the difference of the open circuit potential and the sum of all overvoltages [4, 5, 6, 7, 8]. We are here especially interested in the ohmic overpotential. In many batteries, the electrolyte has by far the highest resistance. In such cases the ohmic losses can be calculated analytically or by potential theory, i.e., by solving a Laplace equation in the electrolyte [9, 10, 11]. More advanced models are used for porous electrodes, but rely on the same Laplace equation [12, 13, 14, 15, 16]. Even if the full potential distribution is computed in the electrolyte, all these models use in the end only a scalar voltage loss to compute the cell potential. They are therefore limited in many ways – and can not account for effects such as the potential drop in poorly conducting electrodes (e.g. Se, S [17]), complicated geometries where much current flows only through a part of the electrolyte, local activation overpotentials, or liquid electrodes where diverging current drives convection [18, 19].

In order to account for all these effects, we need to know the potential distribution in the battery. Fig. 1 shows a one-dimensional example – an example which is strongly simplified. It relies on the macroscopic approach, i.e., we do not resolve the electrochemical double layer, and consider all phases to be electrically neutral111For simulating alternating current response of electrochemical cells it will be necessary to consider the effects of the electrochemical double layer. For more information on the subject, see, e.g., [20]. [13, 21, 3]. Within these assumptions, the potential jumps at both electrode-electrolyte interfaces according to the Nernst equation. As illustrated in fig. 1a the cell voltage is simply the difference of both potential jumps. During operation of the electrochemical cell, the potential jumps are directly reduced by concentration and activation overpotentials. Additionally, the cell potential decreases by the ohmic overpotential, as illustrated in fig. 1b.

1.2 Modeling current and potential distributions

After describing one-dimensional potential profiles in the previous section, we will proceed with a discussion of 3-dimensional modeling. Within continuum mechanics, potential distributions can be obtained by solving a Laplace equation

[TABLE]

where $\varphi$ denotes the electric potential and $\sigma$ the ionic or electronic conductivity. This simple approach is only fully correct if there are no concentration gradients in the electrolyte [25, 6] and if double layer charging is neglected [3]. The challenge of solving equation 1 is the potential jump between electrolyte and electrodes. Generally, two solutions exist: the single- and the multi-potential approaches [26]. The latter relies on the idea of defining separate ionic and electronic potentials [27, 28], i.e., several Laplace equations are solved. The offset, i.e., the jump between the potentials, is typically defined by a volumetric source term in an envelope layer near the interface [29]. Such a model is ideally suited for porous electrodes, such as the catalyst layer in a fuel cell. The porous electrode and electrolyte are treated as superimposed continua [12, 27] of finite thickness [6]; a charge leaving the solid matrix must enter the pore liquid [30]. The superficial charge transfer – described by Faraday’s law – needs to be transferred into only a volumetric source term [12, 21]. Such multi-potential porous electrode models are extremely popular – for more details, see especially [31, 12, 13, 32, 30]. Please note that many researchers exclude the potential jump due to the Nernst potential, but model only the jump due to the activation overpotential.

Compared to the above described porous model, the single-potential approach is less common. It uses only one single potential field for the whole cell. The potential jump is not implemented as a source term, but as an internal jump boundary condition. Such models are also known as “interface models” because they assume the interface to be infinitely thin [3]. We will implement such a model and describe in the following how the equations are set up and how they are solved.

1.3 Solving coupled equation systems with internal jumps

As described in the previous section, the electric potential between electrodes, electrolyte and conductors needs to be coupled. We will denote this by the term “region coupling”. The voltage of electrochemical cells depends, among others, on temperature and concentration. Each of these fields is described by its own equation – and all of them are coupled. We will denote that as “equation coupling”. Newman’s original method for coupling different equations [33] is based on block-matrices, i.e., all equations are discretized in one matrix and solved in one step [34, 35, 36]. This method was slightly extended [37, 38, 39, 40] and used for many different applications. For a good overview about the method, see [41]. An alternative to block-coupling is a segregated solver: it solves each equation separately. As the equations are coupled, this process must be repeated iteratively until convergence is reached. On the one hand, the block-coupled single-matrix approach is surely better suited for highly coupled equations, and maybe for strongly non-linear equations [42]. On the other hand, the segregated approach is easier to implement and needs much less computer memory. Both approaches have advantages and drawbacks – there is not a single, and perfect solution for equation coupling.

Similarly, region coupling can be obtained by solving the potential on one global mesh, or by solving a separate equation for each conductor. The latter is well known in OpenFOAM. Using appropriate boundary conditions, a potential jump at the interface can easily be implemented [43, 44]. However, a relatively time consuming iteration between all regions is necessary. This can be extremely slow [45], because the regions are coupled only at the interface (and not in the volume). Solving the potential in all regions in a coupled way in one matrix is definitely the better way. For possible implementations of internal potential boundary conditions, see [46, 47, 48].

Finally, the equation system can be set-up and solved using the finite-difference (FDM), finite-element (FEM), boundary element or finite-volume method (FVM) [9, 2]. During the first years, the FDM was surely the most widely used approach for modeling potential problems [6]. Today, commercial codes such as COMSOL [49] (FEM), STAR-CCM+ [50] (FVM) and Ansys Fluent [51] (FVM) also allow for solving potential problems with internal jumps. However, their solution algorithm is not always published.

1.4 Our approach

We will model potential distributions with internal jumps using a single-potential approach. The set of equations will be solved in a segregated manner, while a block-matrix will be used for region coupling of the electric potential. We further use a multi-mesh approach [8, 45]. This means, we solve certain variables (as the potential) on a global mesh, but others (such as concentration) only in the electrodes. The equations are discretized using the finite volume method, and the model is implemented in the open source CFD library OpenFOAM v1806 [52].

2 Model

We present in this section the most simple model for the potential distribution in a concentration cell. With little effort, it can be extended to arbitrary electrochemical cells. We use the parent child-mesh technique, i.e., we provide one global mesh for the full geometry and a second mesh for the positive electrode only. Where necessary, we map variables between both meshes. The following simplifications apply:

fluid flow, heat transfer and variation of the layer thickness are neglected

• 2.

the charge transfer overpotential is neglected

• 3.

the electrochemical double layer is not resolved; we assume discrete potential jumps at the interfaces

• 4.

concentration gradients in the electrolyte are neglected [53]

In a first step we solve the concentration $c$ in the positive electrode as

[TABLE]

where $D$ denotes the diffusion coefficient. We apply zero flux boundary conditions at all interfaces except the electrode-electrolyte boundary. Here we set the concentration gradient corresponding to the current density $\boldsymbol{J}$ as

[TABLE]

where $\nabla_{n}$ denotes the interface-normal gradient, $z$ the number of electrons, $F$ the Faraday constant and $\boldsymbol{n}$ the interface normal vector.

The local concentration at the electrode-electrolyte interface determines the jump in potential. The latter could be computed from the Nernst equation using measured activities. However, for simplicity and to avoid numerical instabilities [42], we directly fit the measured open circuit voltage. We save this potential jump at the corresponding face centers and map it to the parent mesh. Please note that the potential jump is applied only at the interface between electrolyte and the positive electrode (and not at the negative electrode).

On the global mesh we solve a Laplace equation for the electric potential $\varphi$ as

[TABLE]

where $\sigma$ denotes the electrical conductivity. Finally, the current density is computed as

[TABLE]

and mapped to the electrode-electrolyte interface where it is needed to compute the boundary condition for the concentration. The potential jump at the electrolyte-positive electrode interface is accounted for when discretizing the gradient and Laplace operator as described in the next section.

We discretize the equations using the implicit Euler scheme for time and second order schemes for all spatial terms. The potential is solved using a PCG, and the concentration by a multigrid solver [54].

3 Discretization of jump condition

Solving the electrical potential on one single mesh, one needs to account for the jumping potential in two terms. The first is the Laplace equation 4 where the potential jump will appear in form of an additional source term. Secondly, when computing the current density by equation 5, both the potential jump and the discrete change of conductivity has to be observed.

3.1 Laplace operator

Within the finite volume method, the Laplace operator can be discretized using the Gauss theorem as [55]

[TABLE]

with $\boldsymbol{S}$ denoting the face area vector and $(\nabla\varphi)_{f}$ the gradient at the face. The face conductivity $\sigma_{f}$ is consistently discretized from the cell centered values222We use the word “cell” in section 3 in the sense of grid cell or control volume. using harmonic interpolation [45]. Denoting the potential in the owner and neighbor cell of face $f$ by $\varphi_{\mathrm{P}}$ and $\varphi_{\mathrm{N}}$ we can also write

[TABLE]

where $\Delta\varphi$ denotes the potential jump at the interface and $\boldsymbol{d}$ the vector connecting both cell centers. Writing the matrix equation as

[TABLE]

we find the off-diagonal coefficients as

[TABLE]

and the diagonal coefficients as

[TABLE]

where the index $r$ denotes related cells, i.e., cells which share a common face333For details, see https://openfoamwiki.net/index.php/OpenFOAM_guide/Matrices_in_OpenFOAM.. Finally, the source term is

[TABLE]

The potential jump $\Delta\varphi$ will appear only at the electrode-electrolyte interface; for all grid cells not touching the interface the source term will therefore be zero.

The above described procedure applies perfectly to orthogonal meshes. It can be easily extended to arbitrary polyhedral control volumes using the overrelaxed correction approach [56, 57, 55, 54, 58].

3.2 Gradient operator

Using the Gauss theorem, the gradient of the electric potential is discretized as [55]

[TABLE]

where $V$ denotes the cell volume and $\varphi_{f}$ the potential at a face. The latter is determined from the potential of the owner (P) and neighbor (N) cell using two continuity conditions. Firstly, we assume the potential itself to be continuous (or jumping) over the face as

[TABLE]

and secondly we ensure continuity of normal currents by

[TABLE]

where $\varphi_{\mathrm{wP}}$ and $\varphi_{\mathrm{wN}}$ denote the face potential in the owner and neighbor cell and $\Delta\varphi$ the jump between both; $\sigma_{\mathrm{P}}$ and $\sigma_{\mathrm{N}}$ denote (cell centered) conductivities, $\boldsymbol{n}_{\mathrm{P}}$ and $\boldsymbol{n}_{\mathrm{N}}$ the face normal vectors and $\varphi_{\mathrm{N}}$ and $\varphi_{\mathrm{P}}$ the cell values of the potential. Combining both above conditions leads to the face potential in the owner cell as

[TABLE]

with

[TABLE]

Here, $\delta_{\mathrm{P}}$ and $\delta_{\mathrm{N}}$ denote the distance between face and cell center for the owner and neighbor cell, respectively.

4 Validation

We use a simple analytical test case to validate our model. A 1-dimensional bar is made of two materials of different conductivity, and an artificial potential jump of 1 V is applied between both. Further, a constant potential is applied at both ends of the bar, such that a potential profile as illustrated in fig. 2 develops. The following parameters are used:

[TABLE]

The analytical solution predicts the potential profile as two lines:

[TABLE]

with

[TABLE]

The simulated electric potential in fig. 2 shows the jump as expected, and fits perfectly to the analytical solution. Moreover, the current density (not shown here) is continuous over the interface, and fits again to the analytical solution.

5 Application to a Li$||$Bi liquid metal battery

Liquid metal batteries are discussed as cheap stationary energy storage for fluctuating renewable energies [59, 60]. We simulate a Li$||$Bi liquid metal battery [61, 62] in 3D (for a 1D model, see [63]) – and compare with measured data. The concentration cell consists of a Bi positive electrode (0.1 mol), a LiF-LiCl-LiI molten salt electrolyte and a Li negative electrode [64]. The lithium is contained in a nickel spiral, and the vessel is made of tantalum to avoid as far as possible any reaction with the active materials and the electrolyte. Fig. 3 shows the setup, and tab. 1 gives the material properties. The cell is heated from below such that the temperature at the bismuth-electrolyte interface stays at approximately 460∘C. Ultra-dry LiI is vacuum dried for 12 h while being heated step wise from 100 to 300∘C. LiF-LiCl is vacuum dried for 12 h at 500∘. Finally, LiF-LiCl-LiI is mixed in eutectic composition [65] and filtered through a quartz frit. The cell is cycled at 1 A with a cycle length of 10 min.

The numerical model is simplified in three ways. Firstly, we insert a very thin gap artificially between electrolyte and vessel, because no current is allowed to flow there. Secondly, we model the negative electrode as it would consist of pure lithium. This is justified because the electric conductivities of nickel and lithium are very similar. Finally, we do not know the exact shape of the electrolyte-negative electrode interface – the exact immersion depth shown in fig. 3 is therefore an assumption. The numerical representation of the cell is shown in fig. 4a: 100 grid cells are used on the diameter of the FVM mesh, and the electrolyte-positive electrode interface is strongly refined. The potential is set to 0 V at the outer cable while a Neumann boundary condition corresponding to 1 A is applied at the negative contact. The open circuit potential at 460∘C is fitted using measurement values [66] as

[TABLE]

for the interval of molar fraction $0.1<x<0.3$. The density of the Li-Bi alloy is computed using Vegard’s law [67, 68].

Fig. 4c illustrates the electric potential field, and fig. 4d the profile along the axis. We clearly see the potential jump corresponding to the open circuit voltage. The current density, as illustrated in fig. 4a shows only an asymmetry at the bottom. Despite all current needs to flow to the lateral cable, the current distribution in the electrolyte is almost symmetric. Finally, fig. 4f shows the 10 min discharge cycle of 1 A, starting with a molar fraction of Li in Bi of $x=0.236$. Experiment and simulation fit very well. The corresponding Li molar fraction in Bi after 1, 5 and 10 minutes is shown in fig. 4e. As expected from the current distribution in fig. 4b, we find more Li in the center of the positive electrode. In reality, this could finally lead to the formation of intermetallic phases there.

6 Summary and outlook

We have discussed different approaches for modeling the cell voltage, as well as the macroscopic potential and current distribution in electrochemical cells. Special attention was paid to potential theory and the coupling of equations and different conductors. Thereafter, we have developed a three dimensional model for internal potential jumps in electrochemical cells, and implemented it into the finite volume method. An internal boundary condition, included into the Laplace and gradient operator allows for arbitrary floating potential jumps at electrode-electrolyte interfaces. As it can compute the full cell potential in 3D at one single mesh, the model is highly efficient, robust and universal. It is fully parallelized and can work in galvano- and potentiostatic mode. After validating the model by a simple analytical formula we have illustrated its capabilities by simulating a Li$||$Bi liquid metal battery and comparing with experimental data. We have correctly predicted the discharge cycle of the cell, and illustrated the current and potential distribution as well as the lithium concentration in bismuth.

Our discretization of the Laplace and gradient terms is quite universal: only a face field describing the jump needs to be provided. Our model can therefore be employed for modeling potential distributions in arbitrary electrochemical cells. Further, it has not escaped our notice that jumps can appear in electromagnetic fields, in temperature and potential fields due to contact resistance [75] or in concentrations fields at interfaces [29], as well. For all those applications, our model might potentially be useful.

In the future we plan to employ the developed model for studying the influence of mass transport on the cell potential of liquid metal batteries [76]. Specifically, we would like to extend the model by convection [18, 19], heat transfer [77, 73] and electrode kinetics.

Acknowledgment

This work was supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) by award number 338560565 as well as a postdoc fellowship of the German Academic Exchange Service (DAAD). The computations were performed on a PC-Cluster at the Center for Information Services and High Performance Computing (ZIH) at TU Dresden and on the cluster “Hemera” at Helmholtz-Zentrum Dresden – Rossendorf. Fruitful discussions with S. Beale, V. Galindo, G. Mutschke and U. Reimer are gratefully acknowledged. Kashif Mushtaq is grateful to the Portuguese Foundation for Science and Technology (FCT) for his PhD grant (reference: PD/BD/128041/2016) and support from the MIT Portugal Program.