Faster Lead-Acid Battery Simulations from Porous-Electrode Theory: I. Physical Model

TL;DR

This paper introduces an advanced porous-electrode model for lead-acid batteries that incorporates previously neglected physical phenomena, enabling more accurate predictions of battery behavior across different discharge rates.

Contribution

It extends existing models by including convection, pressure diffusion, and volume variation effects, providing a more comprehensive and physically accurate simulation framework.

Findings

Model accurately predicts transient voltages and concentration profiles.

Discharge capacity dependence on C-rate deviates from Peukert's law.

Good agreement between model predictions and experimental data.

Abstract

An isothermal porous-electrode model of a discharging lead-acid battery is presented, which includes an extension of concentrated-solution theory that accounts for excluded-volume effects, local pressure variation, and a detailed microscopic water balance. The approach accounts for three typically neglected physical phenomena: convection, pressure diffusion, and variation of liquid volume with state of charge. Rescaling of the governing equations uncovers a set of fundamental dimensionless parameters that control the battery's response. Total volume change during discharge and nonuniform pressure prove to be higher-order effects in cells where variations occur in just one spatial dimension. A numerical solution is developed and exploited to predict transient cell voltages and internal concentration profiles in response to a range of C-rates. The dependence of discharge capacity on C-rate deviates substantially from Peukert's simple power law: charge capacity is concentration-limited at low C-rates, and voltage-limited at high C-rates. The model is fit to experimental data, showing good agreement.

Figures

Tables

Table 1. Table 1: Dimensionless parameters, relative to the C-rate, 𝒞 = 8 ​ A cs ​ i ¯ / Q 𝒞 8 subscript 𝐴 cs ¯ 𝑖 𝑄 \mathcal{C}=8A_{\text{cs}}\bar{i}/Q . 𝒞 d subscript 𝒞 d \mathcal{C}_{\text{d}} is the diffusional C-rate.

Parameter	Definition	Interpretation	Value
			n	sep	p
${𝒞}_{\text{d}}$	$\frac{\overline{i}}{c^{\text{max}}{D}^{\text{max}}F/L}$	$\frac{\text{applied current density}}{\text{limiting current density}}$	$0.60𝒞$
${\iota }_{\text{s}}$	$\frac{{\sigma }^{\text{eff}}RT/FL}{\overline{i}}$	$\frac{\text{Ohmic current scale}}{\text{applied current density}}$	$3.8\times {10}^4/𝒞$	-	$55/𝒞$
${\beta }^{\text{surf}}$	$-\frac{c^{\text{max}}}{n_{\text{e}}}\sum _{k\in \{\text{ce}PbS{O}_4,\text{ce}Pb,\text{ce}Pb{O}_2\}}{\overline{V}}_k{s}_k$	Change in porosity associated with local half-reaction going to completion	$0.084$	-	$-0.064$
$\beta$	$-c^{\text{max}}\Delta \overline{V}$	Change in acid volume fraction associated with local half-reaction going to completion	$0.033$	-	$0.040$
${\gamma }_{\text{dl}}$	$\frac{C_{\text{dl}}RT𝒜L/F\overline{i}}{c^{\text{max}}FL/\overline{i}}$	$\frac{\text{capacitance time-scale}}{\text{discharge time-scale}}$	$2.1\times {10}^{-5}$	-	$1.7\times {10}^{-4}$
${\omega }_c$	$\frac{c^{\text{max}}{M}_{\text{e}}}{{\rho }^{\text{max}}}\left(1-\frac{M_{\text{w}}{\overline{V}}_{\text{e}}}{M_{\text{e}}{\overline{V}}_{\text{w}}}\right)$	Diffusive kinematic relationship coefficient	$0.70$
${\omega }_i$	$\frac{c^{\text{max}}{M}_{\text{e}}}{{\rho }^{\text{max}}}\left(t_{+}^{\text{w}}+\frac{M_{-}}{M_{\text{e}}}\right)$	Migrative kinematic relationship coefficient	$0.41$
${\pi }_{\text{os}}$	$\frac{{\mu }^{\text{max}}{v}^{\text{rxn}}L/d^2}{RT{c}^{\text{max}}}$	$\frac{\text{viscous pressure scale}}{\text{osmotic pressure scale}}$	$3.6\times {10}^{-5}𝒞$

Table 2. Table 1: Dimensional parameters from the literature. Parameters with several values indicate different values in negative electrode (n), separator (s) and positive electrode (p).

Parameter	Value			Units	Reference
	n	sep	p
$l$	$0.9\times {10}^{-3}$	$1.5\times {10}^{-3}$	$1.25\times {10}^{-3}$	m	Private communication
${\hat{\epsilon }}^{\text{max}}$	$0.53$	$0.92$	$0.57$	-	[21]
$H$	$11.4\times {10}^{-2}$			m	Measured
$W$	$6.5\times {10}^{-2}$			m	Measured
$s_{+}$	$-1$	-	$-3$	-	(2.1)
$s_{-}$	$1$	-	$-1$	-	(2.1)
$s_{\text{w}}$	$0$	-	$2$	-	(2.1)
$n_{\text{e}}$	$2$	-	$2$	-	(2.1)
${\overline{V}}_{\text{w}}$	$1.75\times {10}^{-5}$			m3 mol-1	[9]
${\overline{V}}_{\text{e}}$	$4.50\times {10}^{-5}$			m3 mol-1	[9]
${\overline{V}}_{+}$	$1.35\times {10}^{-5}$			m3 mol-1	[34]
${\overline{V}}_{-}$	$3.15\times {10}^{-5}$			m3 mol-1	[34]
${\overline{V}}_{\text{ce}Pb}$	$1.83\times {10}^{-5}$			m3 mol-1	[35]
${\overline{V}}_{\text{ce}Pb{O}_2}$	$2.55\times {10}^{-5}$			m3 mol-1	[35]
${\overline{V}}_{\text{ce}PbS{O}_4}$	$4.82\times {10}^{-5}$			m3 mol-1	[35]
$M_{\text{w}}$	$1.8\times {10}^{-2}$			kg mol-1	[35]
$M_{+}$	$0.1\times {10}^{-2}$			kg mol-1	[35]
$M_{-}$	$9.7\times {10}^{-2}$			kg mol-1	[35]
$F$	$96485$			C mol-1	Faraday constant
$R$	$8.314$			J mol-1 K-1	Ideal gas constant
$T$	$298.15$			K	Room temperature
$t_{+}^{\text{w}}$	$0.72$			-	[16, 13]
$\sigma$	$4.8\times {10}^6$	-	$8\times {10}^4$	S m-1	[14]
$j^{\text{ref}}$	$8\times {10}^{-2}$	-	$6\times {10}^{-3}$	A m-2	[36] (reported by [14])
$c^{\text{max}}$	$5.6\times {10}^3$			mol m-3	Private communication
$𝒜$	$2.6\times {10}^6$	-	$2.05\times {10}^7$	m-1	[36] (reported by [14])
$d$	${10}^{-7}$	-	${10}^{-7}$	m	[14]
$C_{\text{dl}}$	$0.2$	-	$0.2$	F m-2	[21]
$q^0$	$1$			-	Full initial SOC
$Q$	$17$			Ah	Manufacturer-specified

Table 3. Table 2: Dimensional functions of concentration, c ^ ^ 𝑐 \hat{c} (measured in mol/m 3 ) and ε ^ ^ 𝜀 \hat{\varepsilon} (dimensionless). (†) Empirical formulae are given by  [ 14 ] , citing  [ 36 ] , and agree with data of  [ 16 ] . (*) Our fit to data of given reference(s). (‡) m ​ ( c ^ ) = c ^ ​ V ¯ w / [ ( 1 − c ^ ​ V ¯ e ) ​ M w ] 𝑚 ^ 𝑐 ^ 𝑐 subscript ¯ 𝑉 w delimited-[] 1 ^ 𝑐 subscript ¯ 𝑉 e subscript 𝑀 w m(\hat{c})=\hat{c}\bar{V}_{\text{w}}/[(1-\hat{c}\bar{V}_{\text{e}})M_{\text{w}}] .

Function	Formula	Value at $\hat{c}=c^{\text{max}}$	Units	Reference
$\hat{D}$	$\left(1.75+2.6\times {10}^{-4}\hat{c}\right)\times {10}^{-9}$	$3.02\times {10}^{-9}$	m2 s-1	(†)
$\hat{\chi }$	$0.49+4.1\times {10}^{-4}\hat{c}$	$2.8$	-	[16, 37] (*)
$\hat{\kappa }$	$\hat{c}\exp (6.23-1.34\times {10}^{-4}\hat{c}$                    $-1.61\times {10}^{-8}{\hat{c}}^2)\times 10{}^{-4}$	$77$	S m-1	(†)
$\hat{\rho }$	$M_{\text{w}}/{\overline{V}}_{\text{w}}\left(1+\left(M_e{\overline{V}}_{\text{w}}/M_{\text{w}}-{\overline{V}}_e\right)\hat{c}\right)$	$1.32\times {10}^3$	kg m-3	Linear function of $\hat{c}$ [38]
$\widehat{𝒦}$	${\hat{\epsilon }}^3{d}^2/180{\left(1-\hat{\epsilon }\right)}^2$	-	m2	[14]
$\hat{\mu }$	$0.89\times {10}^{-3}+1.11\times {10}^{-7}\hat{c}+3.29\times {10}^{-11}{\hat{c}}^2$	$2.5\times {10}^{-3}$	Pa s	[16] (*)
${\hat{c}}_{\text{w}}$	$\left(1-\hat{c}{\overline{V}}_{\text{e}}\right)/{\overline{V}}_{\text{w}}$	$4.27\times {10}^4$	mol m-3	(2.3)
${\hat{j}}_0$	$j^{\text{ref}}{\left(\frac{\hat{c}}{c^{\text{ref}}}\right)}^{\left|\frac{s_{+}}{n_{\text{e}}}\right|+\left|\frac{s_{-}}{n_{\text{e}}}\right|}{\left(\frac{{\hat{c}}_{\text{w}}}{c_{\text{w}}^{\text{ref}}}\right)}^{\left|\frac{s_{\text{w}}}{n_{\text{e}}}\right|}$	$j^{\text{ref}}$	A m-2	Approach of [17]
${\hat{U}}_{\text{ce}Pb}$	$-0.295-0.074\log m-0.030{\log }^{2}m$           $-0.031{\log }^{3}m-0.012{\log }^{4}m$ (‡)	$-0.41$	V	[32]
${\hat{U}}_{\text{ce}Pb{O}_2}$	$1.628+0.074\log m+0.033{\log }^{2}m$           $+0.043{\log }^{3}m+0.022{\log }^{4}m$ (‡)	$1.76$	V	[32]

Equations

$$\begin{array}[]{l}\ce{Pb+HSO_{4}^{-}<=>[\text{\tiny discharge}][\text{\tiny charge}]${\ce{PbSO_4 + H^+ + 2 e^-}}$}~{}~{}\text{and}\\ \ce{PbO_{2}+HSO_{4}^{-}+3H^{+}+2e^{-}<=>[\text{\tiny discharge}][\text{\tiny charge}]${\ce{PbSO_4 + 2H_2O}}$,}\end{array}$$

$$\hat{c}={\hat{c}}_{+}={\hat{c}}_{-}.$$

$${\hat{c}}_{\text{w}}=\frac{1-{\bar{V}}_{\text{e}}\hat{c}}{{\bar{V}}_{\text{w}}},$$

$$\hat{\rho }={\hat{c}}_{\text{w}}{M}_{\text{w}}+\hat{c}M=\frac{M_{\text{w}}}{{\bar{V}}_{\text{w}}}+\left(\frac{M_{\text{e}}}{{\bar{V}}_{\text{e}}}-\frac{M_{\text{w}}}{{\bar{V}}_{\text{w}}}\right){\bar{V}}_{\text{e}}\hat{c},$$

$$\frac{1}{\mathcal{A}}\frac{\partial \hat{\epsilon }}{\partial \hat{t}}=-{\bar{V}}_{{\text{PbSO}}_4}{\hat{S}}_{{\text{PbSO}}_4}-{\bar{V}}_{\text{Pb}}{\hat{S}}_{\text{Pb}}-{\bar{V}}_{{\text{PbO}}_2}{\hat{S}}_{{\text{PbO}}_2},$$

$$\frac{\partial }{\partial \hat{t}}(\hat{\epsilon }{\hat{c}}_k)=-\hat{\nabla }\cdot {\hat{\mathbf{N}}}_k+\mathcal{A}{\hat{S}}_k,$$

$$\frac{\partial \hat{\epsilon }}{\partial \hat{t}}=-\hat{\nabla }\cdot {\hat{\mathbf{v}}}^{\square }+\mathcal{A}\left({\bar{V}}_{\text{w}}{\hat{S}}_{\text{w}}+{\bar{V}}_{-}{\hat{S}}_{-}+{\bar{V}}_{+}{\hat{S}}_{+}\right),$$

$$\hat{\nabla }\cdot \hat{\mathbf{i}}=\mathcal{A}\hat{j}.$$

$$\hat{\nabla }\cdot {\hat{\mathbf{i}}}_{\text{s}}=-\mathcal{A}\hat{j},$$

$$\begin{gathered} \hat{\mathbf{v}}-{\hat{\mathbf{v}}}^{\square }=\frac{{\bar{V}}_{\text{e}}}{\hat{\rho }}\left(\frac{M_{\text{e}}}{{\bar{V}}_{\text{e}}}-\frac{M_{\text{w}}}{{\bar{V}}_{\text{w}}}\right)\left({\hat{\mathbf{N}}}_{+}-\hat{c}{\hat{\mathbf{v}}}^{\square }\right) \\ +\frac{{\bar{V}}_{-}}{\hat{\rho }}\left(\frac{M_{\text{w}}}{{\bar{V}}_{\text{w}}}-\frac{M_{-}}{{\bar{V}}_{-}}\right)\frac{\hat{\mathbf{i}}}{F}, \end{gathered}$$

$${\hat{\mathbf{N}}}_{+}=-{\hat{D}}^{\text{eff}}\left(\hat{\nabla }\hat{c}-\frac{\hat{\psi }\hat{\nabla }\hat{p}}{RT\hat{\chi }}\right)+\frac{t_{+}^{\text{w}}\hat{\mathbf{i}}}{F}+\hat{c}{\hat{\mathbf{v}}}^{\square },$$

$$\hat{\psi }=\frac{{\bar{V}}_{\text{w}}{\hat{c}}_{\text{T}}{M}_{\text{w}}{\hat{c}}_{\text{w}}{M}_{\text{e}}\hat{c}}{2\hat{\rho }}\left(\frac{{\bar{V}}_{\text{w}}}{M_{\text{w}}}-\frac{{\bar{V}}_{\text{e}}}{M_{\text{e}}}\right).$$

$$\begin{gathered} \hat{\mathbf{i}}=-{\hat{\kappa }}^{\text{eff}}\hat{\nabla }\hat{\Phi }+\frac{2RT\hat{\chi }(1-t_{+}^{\text{w}}){\hat{\kappa }}^{\text{eff}}}{F{\bar{V}}_{\text{w}}{\hat{c}}_{\text{T}}}\hat{\nabla }\ln \hat{c} \\ +\frac{{\hat{\kappa }}^{\text{eff}}}{F\hat{c}}\left[\frac{M_{+}\hat{c}}{\hat{\rho }}-\frac{2\hat{\psi }(1-t_{+}^{\text{w}})}{{\bar{V}}_{\text{w}}{\hat{c}}_{\text{T}}}\right]\hat{\nabla }\hat{p}. \end{gathered}$$

$${\hat{D}}^{\text{eff}}=\hat{D}(\hat{c}){\hat{\epsilon }}^{3/2}\text{and}{\hat{\kappa }}^{\text{eff}}=\hat{\kappa }(\hat{c}){\hat{\epsilon }}^{3/2},$$

$${\hat{\mathbf{i}}}_{\text{s}}=-{\sigma }^{\text{eff}}\hat{\nabla }{\hat{\Phi }}_{\text{s}},$$

$${\sigma }^{\text{eff}}=\sigma {\hat{\epsilon }}_{\text{s}}^{3/2}.$$

$$\hat{\mathbf{v}}=-\frac{\hat{\mathcal{K}}}{\hat{\mu }}\hat{\nabla }\hat{p},$$

$${\hat{S}}_k=-\frac{s_k\hat{j}}{n_{\text{e}}F},$$

$$\hat{\nabla }\cdot {\hat{\mathbf{v}}}^{\square }=-\frac{\Delta \bar{V}\mathcal{A}j}{F}.$$

$$\begin{gathered} n_{\text{e}}\Delta \bar{V}={\bar{V}}_{\text{w}}{s}_{\text{w}}+{\bar{V}}_{-}{s}_{-}+{\bar{V}}_{+}{s}_{+}+{\bar{V}}_{{\text{PbSO}}_4}{s}_{{\text{PbSO}}_4} \\ +{\bar{V}}_{\text{Pb}}{s}_{\text{Pb}}+{\bar{V}}_{{\text{PbO}}_2}{s}_{{\text{PbO}}_2}. \end{gathered}$$

$$\begin{gathered} \frac{\partial }{\partial \hat{t}}\left(\hat{\epsilon }\hat{c}\right)+{\hat{\mathbf{v}}}^{\square }\cdot \hat{\nabla }\hat{c}=\hat{\nabla }\cdot \left[D^{\text{eff}}\left(\hat{\nabla }\hat{c}-\frac{\hat{\psi }\hat{\nabla }\hat{p}}{\hat{\chi }RT}\right)\right] \\ +\left(s+\hat{c}\Delta \bar{V}\right)\frac{\mathcal{A}\hat{j}}{F}. \end{gathered}$$

$$s=-\frac{s_{+}+n_{\text{e}}{t}_{+}^{\text{w}}}{n_{\text{e}}}.$$

$$\hat{j}=2{\hat{j}}_0\sinh \left(\frac{F\hat{\eta }}{RT}\right)+C_{\text{dl}}\frac{\partial ({\hat{\Phi }}_{\text{s}}-\hat{\Phi })}{\partial \hat{t}},$$

$${\hat{j}}_0=j^{\text{ref}}{\left(\frac{\hat{c}}{c^{\text{ref}}}\right)}^{|\frac{s_{+}}{n_{\text{e}}}|+|\frac{s_{-}}{n_{\text{e}}}|}{\left(\frac{{\hat{c}}_{\text{w}}}{c_{\text{w}}^{\text{ref}}}\right)}^{|\frac{s_{\text{w}}}{n_{\text{e}}}|},$$

$$\hat{\eta }={\hat{\Phi }}_{\text{s}}-\hat{\Phi }-\hat{U}.$$

$$\begin{gathered} {\hat{\mathbf{N}}}_{+}\cdot \mathbf{n}=\hat{\mathbf{i}}\cdot \mathbf{n}={\hat{\mathbf{v}}}^{\square }\cdot \mathbf{n}=0 \\ \text{ at x^=0,L, y^=0,W and z^=0.} \end{gathered}$$

$$\begin{gathered} \left({\hat{\mathbf{N}}}_{+}-{\hat{\mathbf{v}}}^{\text{head}}\hat{c}\right)\cdot \mathbf{n}=\hat{\mathbf{i}}\cdot \mathbf{n}=\left({\hat{\mathbf{v}}}^{\square }-{\hat{\mathbf{v}}}^{\text{head}}\right)\cdot \mathbf{n}=0, \\ \hat{p}={\hat{p}}_{\text{ext}},\text{ at z^=H+h^(t^)}. \end{gathered}$$

$${\hat{\Phi }}_{\text{s}}=\{\begin{array}{ll}0& \hat{\mathbf{x}}\in {\text{tab}}_{\text{n}},\\ \hat{V}(\hat{t})& \hat{\mathbf{x}}\in {\text{tab}}_{\text{p}}.\end{array}$$

$$-{\int }_{{\text{tab}}_{\text{n}}}{\hat{\mathbf{i}}}_{\text{s}}\cdot \mathbf{n}\mathrm{d}S={\int }_{{\text{tab}}_{\text{p}}}{\hat{\mathbf{i}}}_{\text{s}}\cdot \mathbf{n}\mathrm{d}S={\hat{\mathrm{I}}}_{\text{circuit}}(\hat{t})/8,$$

$$\begin{gathered} {\int }_0^H{\int }_0^W{\int }_0^{l_{\text{n}}}{\mathcal{A}}_{\text{n}}{\hat{j}}_{\text{n}}\mathrm{d}\hat{x}\mathrm{d}\hat{y}\mathrm{d}\hat{z} \\ =-{\int }_0^H{\int }_0^W{\int }_{L-l_{\text{p}}}^L{\mathcal{A}}_{\text{p}}{\hat{j}}_{\text{p}}\mathrm{d}\hat{x}\mathrm{d}\hat{y}\mathrm{d}\hat{z}={\hat{\mathrm{I}}}_{\text{circuit}}/8. \end{gathered}$$

Full text

1 Introduction

Lead-acid batteries are the most widely used electrochemical storage technology, with applications including car batteries and off-grid energy supply. Models can improve battery management—for example, by minimising overcharge to extend cycle life.

The most rigorous mechanistic approach to battery-cell modelling begins with a detailed microscopic description, wherein the electrolyte and electrodes occupy discrete spatial domains; volume averaging is then performed to produce a macroscopic model [1, 2, 3]. The details are beyond the present scope, but such a homogenization underpins the model presented below.

Our development of a detailed macrohomogeneous model of a typical lead-acid battery augments standard approaches by explicitly considering the balance of water, the variation of acid density with molarity, and the distribution of pressure. As well as ensuring thermodynamic consistency and retaining model closure when systems span multiple spatial dimensions [4], this reveals novel convective phenomena that may occur at high discharge currents.

Nondimensionalization of the model helps to assess the relative importances of different multiphysical phenomena within it. This facilitates numerical implementation because it greatly improves the conditioning of the system, and hence improves the speed of computations, by making most dependent variables close to unity. This paper focuses on the development of a detailed nonisobaric physical model of battery discharge. A novel nondimensionalization is presented, which helps to identify the key dimensionless parameters that control the battery’s transient response. Finally, a numerical procedure is developed to solve the nonlinear system of partial differential equations comprising the model. The results are fit against experimental discharge data and used to show how discharge capacity depends on C-rate. The nonlinear model shows that the battery response does not satisfy Peukert’s law. Instead, the capacity follows one power law at low C-rates, where average acid concentration controls the response, and a different power law at at high C-rates, where overpotentials control the response.

We show below that some of the excluded-volume and pressure effects can be neglected in the detailed physical model, on the basis that certain dimensionless factors are small for the particular materials used in the battery considered here. In cases where dimensionless parameters in a model are small or large, asymptotic analysis can be employed, producing simplified, more computationally efficient models that achieve the same level of physical accuracy. The full nonlinear model presented here serves as the baseline for testing a hierarchy of computationally efficient reduced-order models, which will be developed by perturbation expansion in part II.

2 Model

2.1 Unscaled governing system

Battery configuration and chemistry

Typically, a valve regulated lead-acid battery comprises six 2 V cells wired in series. Figure 1 depicts one such cell, which consists of five lead ($\ce{Pb}$) electrodes and four lead dioxide ($\ce{PbO_{2}}$) electrodes, sandwiched alternatingly around a porous, electrically insulating separator to produce eight electrode pairs, wired in parallel at the top edge of the electrode pile.

The pile has height $H$, depth $W$, and cross-sectional area $A_{\text{cs}}=HW$. The negative ($\ce{Pb}$) and positive ($\ce{PbO_{2}}$) electrodes have half-widths $l_{\text{n}}$ and $l_{\text{p}}$ respectively, and the separator has width $l_{\text{sep}}$, giving each internal electrode pair total width $L={l_{\text{n}}+l_{\text{sep}}+l_{\text{p}}}$. The layers repeat periodically, allowing the whole pile to be modelled via analysis of one electrode pair. Some thermal effects have been documented experimentally [5, 6] and considered theoretically [7, 8]; we assume an isothermal system for simplicity.

The electrodes are porous, permeated by an aqueous sulfuric acid ($\ce{H_{2}SO_{4}}$) electrolyte that also permeates the separator. The negative and positive half-reactions are

[TABLE]

respectively. In both reactions, solid lead sulfate ($\ce{PbSO_{4}}$) forms as bisulfate anions ($\ce{HSO_{4}^{-}}$) leave the liquid.

Without competitive chemistry, reactions (2.1) would simply reverse during recharge, as indicated. However, some authors [9, 10, 7, 11] have observed that gas evolution can become important when recharging. To avoid the need for considering such side reactions, the present analysis is limited to discharges.

Liquid thermodynamics

In water ($\ce{H_{2}O}$), $\ce{H_{2}SO_{4}}$ dissociates mainly to hydrogen cations ($\ce{H^{+}}$) and $\ce{HSO_{4}^{-}}$; in keeping with prior models [9, 12, 13, 14], speciation to sulfate is neglected. Thus the pore-filling liquid comprises $\ce{H_{2}O}$, $\ce{HSO_{4}^{-}}$, and $\ce{H^{+}}$, with partial molar volumes $\bar{V}_{\text{w}}$, $\bar{V}_{-}$ and $\bar{V}_{+}$, molar masses $M_{\text{w}}$, $M_{-}$, and $M_{+}$, and equivalent charges [math], $-1$, and $+1$, respectively. Unscaled variables are denoted with a hat: for example, the species molarities are $\hat{c}_{\text{w}}$, $\hat{c}_{-}$, and $\hat{c}_{+}$.

On any scale much larger than the Debye length, local electroneutrality holds in the liquid phase. Both ion concentrations then relate to the $\ce{H_{2}SO_{4}}$ molarity $\hat{c}$ through

[TABLE]

In an isothermal state, thermodynamic consistency of the liquid volume requires that $\hat{c}_{\text{w}}$ depends on $\hat{c}$ alone if the bulk modulus of the liquid is very large [15], through

[TABLE]

where $\bar{V}_{\text{e}}=\bar{V}_{+}+\bar{V}_{-}$ is the partial molar volume of the acid. Consequently the total liquid molarity $\hat{c}_{\text{T}}=2\hat{c}+\hat{c}_{\text{w}}$ also depends only on $\hat{c}$. Acid molarity further determines the mass density $\hat{\rho}$ of the liquid, because

[TABLE]

in which $M_{\text{e}}=M_{+}+M_{-}$ stands for the acid’s molar mass. Experiments show that $\hat{\rho}$ varies almost linearly with $\hat{c}$ [16], so the partial molar volumes can be assumed constant.

The mechanical state of the liquid is described by the external pressure $\hat{p}$, while its mixing free energy is parameterized by a thermodynamic Darken factor $\hat{\chi}\left(\hat{c}\right)$.

Balances

The position $\hat{x}$ within each electrode pair traverses three subdomains: the negative electrode ($0<\hat{x}<l_{\text{n}}$), separator ($l_{\text{n}}<\hat{x}<L-l_{\text{p}}$), and positive electrode ($L-l_{\text{p}}<\hat{x}<L$). The model structure is identical in each subdomain. Parameters describing the electrolyte are similar everywhere, while those describing pore geometry generally differ among subdomains; quantities that parameterise reactions differ in the negative- and positive-electrode subdomains, and vanish in the separator.

Liquid volume fraction $\hat{\varepsilon}$ and reactive solid area per volume $\mathcal{A}$ characterize the homogenized geometry within the electrode pair. In electrode subdomains, deposition (removal) of solid $\ce{PbSO_{4}}$ on pore surfaces is accompanied by removal (deposition) of solid $\ce{Pb}$ or $\ce{PbO_{2}}$, so these geometric parameters can generally vary locally. Solid $\ce{PbSO_{4}}$ does not tend to deposit as a compact thin film, so mechanistic lead-acid battery models usually let $\mathcal{A}$ vary with state of charge [9, 12, 13, 14, 11]. Since the functionality of this variation is disputed, we instead let the area be constant, following the approach of Liu et al., who showed this is sufficient to model Li/air-battery discharge [17]. From this perspective $\mathcal{A}$ describes an immobile Gibbs dividing surface that partitions the layers of Pb or $\ce{PbO_{2}}$ that contribute to the battery’s charge state from the current-collecting Pb layer beneath them, which does not. Thus the time change in pore volume relates to the molar volumes of the solid species in scheme (2.1) through

[TABLE]

in which $\bar{V}_{k}$ is the partial molar volume of species $k$ and $\hat{S}_{k}$ is the rate at which $k$ is generated by interfacial reactions per unit of pore area.

After homogenization, the local mass balance of species $k\in\left\{\text{w},-,+\right\}$ in the pore-filling liquid phase implies that

[TABLE]

where $\hat{\mathbf{N}}_{k}$ represents the molar flux of $k$. With electroneutrality condition (2.2) and equation of state (2.3), all three balances (2.6) combine to show liquid-volume continuity,

[TABLE]

in which $\hat{\mathbf{v}}^{\square}=\bar{V}_{\text{w}}\hat{\mathbf{N}}_{\text{w}}+\bar{V}_{-}\hat{\mathbf{N}}_{-}+\bar{V}_{+}\hat{\mathbf{N}}_{+}$ signifies the volume-averaged liquid velocity.

Under Faraday’s law, ion balances (2.6) also combine to demonstrate charge continuity in the liquid,

[TABLE]

Letting $F$ be Faraday’s constant, the liquid-phase current density is $\hat{\mathbf{i}}=F(\hat{\mathbf{N}}_{+}-\hat{\mathbf{N}}_{-})$ and the current density associated with interfacial charge exchange is $\hat{j}=F(\hat{S}_{+}-\hat{S}_{-})$. (Positive $\hat{j}$ flows into pores.) Note that any current leaving the liquid at a given location enters the solid there. Thus

[TABLE]

where $\hat{\mathbf{i}}_{\text{s}}$ is the solid-phase current density.

In a general isothermal setting, liquid convection is determined by a momentum balance, such as Cauchy’s equation [4, 15]. This governs the distribution of momentum density $\hat{\rho}\hat{\mathbf{v}}$, naturally expressed in terms of the mass-averaged velocity $\hat{\mathbf{v}}$. The kinematic relation

[TABLE]

specifies how $\hat{\mathbf{v}}$ must relate to $\hat{\mathbf{v}}^{\square}$.

Flux constitutive laws

Two Onsager–Stefan–Maxwell flux laws govern transport in the liquid phase. The law for the thermodynamic force on water [15] can be inverted [18] to give

[TABLE]

where $R$ is the gas constant, $T$ is the absolute temperature, $t^{\text{w}}_{+}$ is the cation transference number relative to the water velocity, and $\hat{D}^{\text{eff}}$ is the effective diffusivity of acid in water; the pressure-diffusion factor $\hat{\psi}$ depends on $\hat{c}$ through

[TABLE]

(To put equation (2.11) in the form given, one must have that $\bar{V}_{+}=(1-t^{\text{w}}_{+})\bar{V}_{\text{e}}$ and $\bar{V}_{-}=t^{\text{w}}_{+}\bar{V}_{\text{e}}$ [19].) The law for the thermodynamic force on cations can be linearly transformed to produce a current/voltage relation:

[TABLE]

Here $\hat{\kappa}^{\text{eff}}$ is the effective conductivity; $\hat{\Phi}$ is the potential measured by a reversible hydrogen electrode at $\hat{p}$ [20]. Effective properties appear in equations (2.11) and (2.13) because pore connectivity affects apparent transport rates. We let

[TABLE]

following Bruggeman’s tortuosity correlation.

The solid phase is electronically but not ionically conductive. Thus the current density there relates to the solid-phase potential, $\hat{\Phi}_{\text{s}}$, through Ohm’s law,

[TABLE]

where $\sigma^{\text{eff}}$ is the effective conductivity,

[TABLE]

Here $\sigma^{\text{eff}}$ is constant because $\hat{\varepsilon}_{\text{s}}$, the volume fraction of nonreacting Pb bounded within electron-exchange surface $\mathcal{A}$, does not vary. Note that $\hat{\varepsilon}_{\text{s}}$ does not count the volume fraction occupied by solid reactants, so $\hat{\varepsilon}_{\text{s}}\neq 1-\hat{\varepsilon}$.

A constitutive law for liquid stress closes the model [4]. Below, terms involving momentum will be shown to be negligible in the first approximation. To analyse the scale of these terms, it will suffice to assume that current density induces a flow with low Reynolds number, in which case the homogenization of Cauchy’s equation produces Darcy’s law,

[TABLE]

in which $\hat{\mu}$ is the liquid viscosity; pore geometry determines the Kozeny–Carman permeability $\hat{\mathcal{K}}$.

Interfacial constitutive laws

Electrons are the only solid-phase charge carrier in scheme (2.1), making interfacial electronic current a proxy for the half-reaction rates. No interfacial reactions occur in the separator domain, so $\hat{j}=0$ uniformly there and the reactive area $\mathcal{A}$ does not need to be defined there.

In an electrode subdomain where a single half-reaction involving $n_{\text{e}}$ electrons occurs, the molar flux of reactants is $\hat{j}/n_{\text{e}}F$ and hence every $\hat{S}_{k}$ is given by

[TABLE]

where $s_{k}$ is the signed stoichiometric coefficient of species $k$ in the half-reaction. (For a reduction half-reaction, $s_{k}$ is positive for a product and negative for a reactant.)

Relationship (2.18) permits simpler expressions to be used in place of the general material balances (2.6). First, equations (2.5), (2.6), (2.7) and (2.18) combine to give

[TABLE]

This liquid-phase volume balance introduces the volume of reaction $\Delta\bar{V}$, related to half-reaction stoichiometry by

[TABLE]

Second, one can combine equations (2.2), (2.5), (2.6), (2.8), (2.11), (2.18), and (2.19) to show that the acid is governed by a form of the convective diffusion equation:

[TABLE]

The generation term here includes a single additional parameter,

[TABLE]

Three aspects of the acid balance (2.21) are new. First, the convection term, and second, the volume of reaction $\Delta\bar{V}$, appear because state equation (2.3) imposes constraints on balances (2.6). Finally, a pressure-diffusion term appears because the flux laws are based on thermodynamic forces.

To complete the model, the current density $\hat{j}$ across the surface $\mathcal{A}$ is governed by a chemical-kinetic constitutive law. Generally such laws involve the voltage difference between the liquid and solid, the equilibrium potential of the half-reaction, and the chemical activities of the reactants. We assume the half-reactions in scheme (2.1) are elementary, following Butler–Volmer kinetics. Butler–Volmer laws naturally include a symmetry factor, which we take to be one half [13], yielding

[TABLE]

where $\hat{j}_{0}$ is the concentration-dependent exchange-current density,

[TABLE]

and $\hat{\eta}$ is the surface overpotential

[TABLE]

Here $\hat{U}$ stands for the half-reaction’s open-circuit potential (OCP) relative to a particular reference electrode. (The reference must be the same for all half-reactions.) Terms involving interfacial capacitance $C_{\textrm{dl}}$ help to smooth out numerics, but have a negligible effect on model predictions because the capacitive time-scale is very short [21].

Boundary conditions

Equations (2.5), (2.8) to (2.10), and (2.13) to (2.25) comprise a three-dimensional model with closure at every interior point in an electrode pair. Symmetry and insulating boundary conditions demand that no species in the liquid phase flows through the centers, sides, and bottom of the electrode pair, so

[TABLE]

Kinematic relation (2.10) shows that $\hat{\mathbf{v}}\cdot\mathbf{n}=0$ at these boundaries too. Flux laws (2.11), (2.13), and (2.17) further require the surface-normal gradients $\partial\hat{c}/\partial n$, $\partial\hat{p}/\partial n$ and $\partial\hat{\Phi}/\partial n$ to vanish at these boundaries.

Above the electrodes, there is a region of free electrolyte, with height $\hat{h}(\hat{t})$. At the top surface of this region, we impose a known external pressure $\hat{p}_{\text{ext}}$, and the absence of flux relative to the surface, which moves with velocity ${\mathbf{\hat{v}}^{\text{head}}=\left(\partial\hat{h}/\partial\hat{t}\right)\mathbf{e}_{z}}$:

[TABLE]

Note the final condition in (2.27) determines the a priori unknown height $\hat{h}(\hat{t})$.

Hereafter, subscripts n and p denote property values in the negative- and positive-electrode subdomains. We choose the negative electrode to be the ground state, and define the cell voltage to be the potential at the positive electrode tab:

[TABLE]

One can either control the voltage, or consider a current-controlled discharge where the voltage is determined by

[TABLE]

where $\hat{\mathrm{I}}_{\text{circuit}}$ is the current drawn from the battery, which is positive for a discharge; the factor of 8 appears because the cell comprises eight electrode pairs in parallel.

This paper focuses on experiments under ‘galvanic control’, following condition (2.29), which allow $\hat{\mathrm{I}}_{\text{circuit}}\left(\hat{t}\right)$ to be any function of time. Since six cells are connected in series, the voltage in the external circuit is $\hat{V}_{\text{circuit}}=6\hat{V}$.

Relationships between subdomains

The liquid phase permeates all three subdomains. Therefore scalar variables $\hat{c}$, $\hat{p}$, and $\hat{\Phi}$, as well as the normal components of all flux vectors, are continuous across electrode/separator boundaries.

There is no solid-phase current at either edge of the separator subdomain or pore-surface charge exchange within it, so $\hat{{i}}_{\text{s}}$ vanishes uniformly there. Since the separator subdomain electronically isolates the positive and negative electrodes, $\hat{\Phi}_{\text{s}}$ is not continuous across it.

Integrating the interfacial current distributions in equation (2.9) and applying the divergence theorem, boundary conditions (2.29) and the fact that $\hat{{i}}_{\text{s}}$ vanishes at the electrode/separator interfaces leads to integral constraints,

[TABLE]

In short, these say that the total current leaving the negative electrode domain must enter the positive electrode domain. Expressions (2.30) will help to analyze the scales of pressure and velocity in (2.42a).

Initial conditions

The initial electrolyte concentration and electrode porosities are spatially uniform, but depend on the state of charge, $q$, which we define as

[TABLE]

Let $c^{\text{max}}$, $\varepsilon^{\text{max}}_{\text{n}}$ and $\varepsilon^{\text{max}}_{\text{p}}$ be the values of electrolyte concentration, negative electrode porosity, and positive electrode porosity, respectively, at full state of charge. In a lead-acid battery, the state of charge is closely linked to the concentration of the electrolyte. Hence $q$ is chosen to be unity when the concentration of the electrolyte is at its maximum value, $c^{\text{max}}$, and zero when the concentration of the electrolyte is zero, so that the initial conditions are

[TABLE]

Parameters $q^{\text{max}}$ and $\varepsilon^{\Delta}$ are chosen to make (2.31) and (2.32) consistent with (2.21) and (2.5) (see A for details):

[TABLE]

Finally, because interfacial capacitance effects are included, initial conditions are needed for the potentials; we choose spatially homogeneous values such that $\hat{j}=0$ at $\hat{t}=0$:

[TABLE]

2.2 Dimensionless model

Nondimensionalization of the system presented in the Unscaled governing system section helps to determine the dominant effects in the system. If $\partial\hat{h}/\partial\hat{t}$ is sufficiently small and the electrode conductivity is sufficiently high, one can assume uniformity in the plane normal to $\hat{x}$, reducing the problem to one spatial dimension. In this case, the tabs cover the whole electrodes at $\hat{x}=0$ and $\hat{x}=L$, and so the boundary condition (2.29) becomes

[TABLE]

which introduces the variable $\hat{i}_{\text{cell}}$ to stand for the total current density at the cell level. Let $\hat{N}_{+}$, $\hat{i}$, $\hat{i}_{\text{s}}$, $\hat{v}^{\square}$ and $\hat{v}$ represent the $\hat{x}$-components of vectors $\hat{\mathbf{N}}_{+}$, $\hat{\mathbf{i}}$, $\hat{\mathbf{i}}_{\text{s}}$, $\hat{\mathbf{v}}^{\square}$ and $\hat{\mathbf{v}}$. A first set of dimensionless variables is formed by the natural rescalings

[TABLE]

Reduce the number of parameters by identifying the scale of interfacial current density, $\bar{i}/\mathcal{A}L$, and discharge time-scale, $c^{\text{max}}FL/\bar{i}$. Hence nondimensionalize interfacial current density, exchange-current density, and time as

[TABLE]

After nondimensionalization, ${j}_{0,\text{n}}$ and ${j}_{0,\text{p}}$ are $\mathcal{O}(1)$ functions. To nondimensionalize the potentials, note that the dominant exponents in the hyperbolic-sine terms of (2.23) should be $\mathcal{O}(1)$, since the terms multiplying the sinh functions are all $\mathcal{O}(1)$. Equation (A.1) can be exploited to define the dimensionless open-circuit potentials

[TABLE]

Then, to ensure that the exponents in equation (2.23) are $\mathcal{O}(1)$ and noting that $\hat{\Phi}_{\text{s}}=0$ at $\hat{x}=0$, introduce the dimensionless potentials

[TABLE]

Quantities $\hat{D}$ (and hence $\hat{D}^{\text{eff}}$) and $\hat{c}_{\text{w}}$ are appropriately scaled with their values at $c=c^{\text{max}}$:

[TABLE]

where $D^{\text{max}}=\hat{D}(c^{\text{max}})$ and $c_{\text{w}}^{\text{max}}=\hat{c}_{\text{w}}(c^{\text{max}})$. The conductivity and Darken thermodynamic factor rescale as

[TABLE]

where the quantity ${\alpha=-(2\bar{V}_{\text{w}}-\bar{V}_{\text{e}})c^{\text{max}}}$ is defined by Liu et al. [4] as the excluded-volume number.

In equation (2.24) for the exchange-current density, the reference concentrations are taken to be ${\hat{c}^{\text{ref}}=\hat{c}^{\text{max}}}$ and ${\hat{c}_{\text{w}}^{\text{ref}}=\hat{c}_{\text{w}}^{\text{max}}}$. The dimensionless widths of the negative electrode, separator and positive electrode become $\ell_{\text{n}}=l_{\text{n}}/L$, $\ell_{\text{sep}}=l_{\text{sep}}/L$, and $\ell_{\text{p}}=l_{\text{p}}/L$, respectively.

In the velocity equations (2.19), (2.17) and (2.10), $\hat{\rho}$ and $\hat{\mu}$ are scaled with their values at $c^{\text{max}}$, and $\hat{\mathcal{K}}$, with $d^{2}$, where $d$ is the characteristic pore size. The reaction velocity scale is $v^{\text{rxn}}=\bar{i}/c^{\text{max}}F$, and the Darcy pressure scale $\mu^{\text{max}}v^{\text{rxn}}L/d^{2}$, so that

[TABLE]

where $p^{\text{ref}}$ is a reference pressure (e.g. atmospheric). This scaling transforms the equations governing $v^{\square}$, $v$ and $p$ to

[TABLE]

in which the dimensionless parameters $\beta$, $\omega_{c}$ and $\omega_{i}$ are defined in Table 1. Further, the effect of pressure gradients in equations (2.13) and (2.21) is smaller than the effect of concentration gradients (ignoring the $\mathcal{O}(1)$ functions of concentration $\chi$ and $\psi$) by a factor of

[TABLE]

Since $\beta$ takes different values in the two electrodes, equation (2.43a) does not admit a solution where $v^{\square}$ vanishes at both electrode centers. However, as can be seen in Table 1, the dimensionless parameters $\beta$ and $\pi_{\text{os}}$ are small. We assume a limit where both of these parameters are zero (and hence $\partial\hat{h}/\partial\hat{t}=0$), in which case $v^{\square}$ vanishes everywhere, and $v$ and $p$ decouple from the other variables. The full model can then be solved to find the voltage without needing the velocity and pressure.111The limit of finite $\beta$ will be explored in a future, two-dimensional, model.

Having decoupled flow velocity and pressure from the rest of the model, the following dimensionless system for $c$, $j$, $\varepsilon$, ${i}$, $\Phi$, ${i}_{\text{s}}$ and $\Phi_{\text{s}}$ results:

[TABLE]

Composition dependences of the properties $D^{\textrm{eff}}$, $\chi$, $\kappa^{\textrm{eff}}$, $j_{0}$, $U_{\textrm{Pb}}$, and $U_{\textrm{PbO}_{2}}$ are established through the functions listed in table 2. The four dominant dimensionless parameters, $\mathcal{C}_{\text{d}}$, $\iota_{\text{s}}$, $\beta^{\text{surf}}$ and $\gamma_{\text{dl}}$, are defined in Table 1, which also states physical interpretations and typical values. In particular, the diffusional C-rate can be written as

[TABLE]

where ${8A_{\text{cs}}c^{\text{max}}FL}$ is the volumetric capacity of the battery (in Ah), $Q$ is the capacity of the battery (in Ah), and ${\mathcal{C}=8A_{\text{cs}}\bar{i}/Q}$ is the C-rate of operation. Alternatively, one can identify $\mathcal{C}_{\text{d}}$ to be the ratio between the applied current scale, $\bar{i}$, and the scale of the liquid-phase limiting current, $i_{\text{L}}=c^{\text{max}}D^{\text{max}}F/L$.

In the Results section, $q^{0}$ will be taken to be unity (the battery starts from a fully charged state) unless explicitly stated.

3 Numerical solution

The system of equations (2.45) was solved numerically. Code used to solve the model and generate the results below is available publicly on GitHub [22].

To facilitate numerical solution, the model was first manipulated into a form suitable for application of the Finite Volume Method. Letting $\phi=\Phi_{\text{s}}-\Phi$ and noting that $i_{\text{s}}=\mathrm{i}_{\text{cell}}-i$, one can replace equations (2.45d) to (2.45f) with

[TABLE]

We also eliminate $\partial\Phi/\partial x$ from equation (3.1) to find $i$ as a functional of $c$ and $\phi$:

[TABLE]

Equation system (3.2) is solved by discretising the spatial domain using Finite Volumes, choosing the spatial discretisation to be uniform within each subdomain and as uniform as possible across subdomains The results are robust to the total number of points chosen.

Having discretised in space, the resulting system of transient ordinary differential equations is solved using scipy.integrate in Python [23]. Finally, $\Phi$ is calculated as

[TABLE]

where the final term comes from the fact that $\Phi=-\phi$ at $x=0$. The voltage drop across the electrode pair is computed with

[TABLE]

Crucially, the system of partial differential equations (3.2) is much easier to solve than the differential-algebraic system (2.45).

4 Results

Figure 2 shows the (dimensional) voltage against capacity used for a range of C-rates. As is to be expected, the total capacity available decreases as the C-rate increases. This dependence is further elucidated by exploring the total capacity of the battery for constant-current discharges across a higher range of C-rates, summarized on Figure 3. By Peukert’s Law, one would expect the graph to be linear on this log-log axis. However, in this case deviations from linearity occur because there are two distinct capacity-limiting mechanisms. At low C-rates (below 1C), the battery capacity is concentration-limited, with full discharge occurring when the electrolyte concentration reaches zero. At high C-rates, the battery is voltage-limited, because the cut-off voltage of 10.5V is achieved before the bulk concentration falls to zero.

The model also allows internal variables to be explored, particularly the local water concentration. Electrolyte and water concentrations are depicted in Figure 4. At a very low C-rate of 0.1C (Figure 4a), both concentrations remain almost uniform throughout the discharge; electrolyte concentration decreases, and water concentration increases proportionally. At a higher C-rate of 0.5C (Figure 4b) and 2C (Figure 4c), the concentrations become spatially inhomogeneous, which leads to concentration overpotentials that limit the accessible capacity.

A Derivative-Free Gauss-Newton method [24] was also used to fit the model to data from a series of constant-current discharges of a 17 Ah BBOXX Solar Home battery at intervals of 0.5 A from 3 A to 0.5 A (Figure 5). Each constant-current discharge is followed by a two-hour rest period during which the current is zero. The fit is good, but this approach is slow since it requires solving the full PDE system at each iteration. A faster approach will be developed in part II.

5 Conclusions

Three novel phenomena were included in a porous-electrode model for lead-acid batteries. First, the mass-averaged and volume-averaged velocities of the electrolyte were both considered, the former associated with momentum transport, and the latter with kinematics. Due to density variation in the liquid, and volume changes associated with the electrode reactions, neither of these velocities remains solenoidal after homogenisation. Second, an extra convective term, associated with the volume-averaged velocity, drives cation transport. Third, a pressure-diffusion term drives cation transport, and appears also in the MacInnes equation (modified Ohm’s law) describing the liquid. Although these terms are small in magnitude for lead-acid batteries, they could be important for other chemistries where large volume changes occur during charge/discharge, such as lithium-ion batteries with silicon anodes [25, 26, 27, 28, 29, 30], or to understand the impedance signature of electromechanical/transport coupling [31]. Nondimensionalisation of this model allows us to identify key parameter groupings, and could also easily be extended to other models and chemistries.

Two distinct mechanisms determine the total cell capacity: concentration limitation at low C-rates, and voltage limitation at high C-rates. In addition to capacity-limitation and voltage-limitation, there may be additional capacity limitations from additional physical effects, such as pore occlusion at very high C-rates [17].

The model developed in this paper provides physical detail about the electrochemical processes occurring in a lead-acid battery during discharge, but is ultimately too computationally intensive to be used for advanced battery management systems. In part II, asymptotic analysis is used to derive three simplified models valid at low-to-moderate discharge rates. These can be solved much faster than the detailed model developed here, while giving additional physical insights.

Acknowledgements

This publication is based on work supported by the EPSRC center For Doctoral Training in Industrially Focused Mathematical Modelling (EP/L015803/1) in collaboration with BBOXX. JC, CP, DH and CM acknowledge funding from the Faraday Institution (EP/S003053/1).

List of symbols

Variables

$c$

concentration mol m*-3*

$\varepsilon$

porosity -

$j$

interfacial current density A m*-2*

$\mathbf{i}$

current density (3D) A m*-2*

$i$

current density in $x$-direction A m*-2*

$\mathbf{v}$

velocity (3D) m s*-1*

$v$

velocity in $x$-direction m s*-1*

$\mathbf{N}$

ion flux (3D) mol m*-2* s*-1*

$p$

pressure Pa

$\Phi$

potential V

Subscripts

n

in negative electrode

sep

in separator

p

in positive electrode

$+$

of cations

$-$

of anions

w

of solvent (water)

e

of electrolyte

s

of solid (electrodes)

Superscripts

[math]

initial

max

maximum

eff

effective

surf

surface

$\square$

convective

Accents

$\hat{}$

dimensional

Appendix A Parameters

The dimensional parameters are given in Table 1, and the concentration dependences of coefficients are laid out in Table 2. Formulae for open-circuit potentials were obtained empirically by Bode [32]. Note that these could be written as

[TABLE]

to reflect their relation to the Nernst equation (e.g. [33]), with $U^{\ominus}_{\ce{Pb}}=-0.295$ and $U^{\ominus}_{\ce{PbO_{2}}}=1.628$.

Consistent initial conditions

For consistent initial conditions, we take ‘starting at $x\%$ SOC’ to mean an internally equilibrated (i.e. spatially homogeneous) initial state at uniform concentration corresponding to this charge. The following determines initial conditions for $\hat{c}$ and $\hat{\varepsilon}$ consistent with this choice.

Considering a one-dimensional model and integrating (2.21) in $\hat{x}$ across the whole electrode pair and using the no-flux conditions (2.26), and the integral condition (2.30),

[TABLE]

Both sides of (A.2) should be zero when $q=0$; hence choose

[TABLE]

which, with the parameter values in Table 1, gives a maximum capacity of $26.1$ Ah for the battery. This compares favourably to the stated battery capacity of $17$ Ah.

Now, integrating (2.5) in $\hat{x}$ across the whole negative electrode and using the integral condition (2.30),

[TABLE]

Then, assuming that $\varepsilon^{0}_{\text{n}}$ is spatially uniform,

[TABLE]

Similarly,

[TABLE]

Appendix B Velocity

Integrating (2.43a) for the volume-averaged velocity from $x=0$ to $x=1$ and using the integral condition (2.30) for the interfacial current density and (2.45h) for $v^{\square}$ at $x=0$, one finds that

[TABLE]

which contradicts (2.45h) for $v^{\square}$ at $x=1$ since the sum of the volume changes across the whole cell is non-zero. It follows that the model can have no exact solution in a one-dimensional setting. This can be resolved by considering a multi-dimensional problem with a free electrolyte surface; this solution is beyond the present scope.