Faster Lead-Acid Battery Simulations from Porous-Electrode Theory: II. Asymptotic Analysis

TL;DR

This paper develops and analyzes simplified models for lead-acid batteries using asymptotic methods, balancing physical accuracy and computational efficiency for various discharge rates.

Contribution

It introduces three reduced-order models derived from a detailed porous-electrode theory, enabling faster simulations with physical interpretability.

Findings

Lumped model accurate below 0.1C discharge rate

Higher-order model retains accuracy up to 5C

Parameter estimation achieved through experimental data fitting

Abstract

Electrochemical and equivalent-circuit modelling are the two most popular approaches to battery simulation, but the former is computationally expensive and the latter provides limited physical insight. A theoretical middle ground would be useful to support battery management, on-line diagnostics, and cell design. We analyse a thermodynamically consistent, isothermal porous-electrode model of a discharging lead-acid battery. Asymptotic analysis of this full model produces three reduced-order models, which relate the electrical behaviour to microscopic material properties, but simulate discharge at speeds approaching an equivalent circuit. A lumped-parameter model, which neglects spatial property variations, proves accurate for C-rates below 0.1C, while a spatially resolved higher-order solution retains accuracy up to 5C. The problem of parameter estimation is addressed by fitting experimental data with the reduced-order models.

Figures

Tables

Table 1. Table 1: Dimensionless parameters, relative to the C-rate, 𝒞 = I circuit / Q 𝒞 subscript 𝐼 circuit 𝑄 {\mathcal{C}=I_{\text{circuit}}/Q} . Further details and interpretations can be found in part I.

Parameter	Value
	n	sep	p
${𝒞}_{\text{d}}$	$0.60𝒞$
${\iota }_{\text{s}}$	$3.8\times {10}^4/𝒞$	-	$55/𝒞$
${\beta }^{\text{surf}}$	$0.084$	-	$-0.064$
${\gamma }_{\text{dl}}$	$2.1\times {10}^{-5}$	-	$1.7\times {10}^{-4}$
$\ell$	$0.25$	$0.41$	$0.34$
$s$	$-0.2$	-	$0.8$
$q^0$	$1$
${\epsilon }^{\text{max}}$	$0.53$	$0.92$	$0.57$
${\epsilon }^{\Delta }$	$0.24$	-	$-0.13$

Table 2. Table 2: Speed comparison for the four models for a constant-current discharge using the parameters from literature (Tables  1 and 2 ). Time is CPU time in seconds, obtained on an AMD FX(tm)-4350 Quad-Core Processor, and averaged over 100 runs; speed-up is ratio of numerical time to own time.

	0.1C		0.5C		2C		5C
Solution Method	Time	Speed-up	Time	Speed-up	Time	Speed-up	Time	Speed-up
Numerical	0.643	-	0.265	-	1.3	-	1.19	-
Composite	0.050	13	0.048	5	0.051	26	0.292	4
FOQS	0.005	126	0.005	55	0.005	261	0.005	239
LOQS	0.002	407	0.001	183	0.001	887	0.001	805

Table 3. Table 3: Performance of fitting algorithms for each model. Error is the sum-of-squares of voltage prediction errors, in Volts; time is CPU time in seconds, obtained on an AMD FX(tm)-4350 Quad-Core Processor.

	LOQS		FOQS		Composite		Numerical
Algorithm	Error	Time	Error	Time	Error	Time	Error	Time
SciPy leastsquares	14.03	36	13.43	48	328.72	122	351.52	1443
DFO-GN	14.03	74	13.43	184	17.04	96	9.93	1144

Table 4. Table 4: Parameters obtained from SciPy with the FOQS model, and DFO-GN with the Numerical model. All algorithms given Tables  1 and 2 as initial data. A full version of this table is given in D .

	${\epsilon }_{\text{n,p}}^{\text{max}}$	${\epsilon }_{\text{sep}}^{\text{max}}$	$j_{\text{n}}^{\text{ref}}$	$R_{\text{circuit}}$	$q_{\text{2.5A}}^0$	$q_{\text{2A}}^0$	$q_{\text{1.5A}}^0$	$q_{\text{1A}}^0$	$q_{\text{0.5A}}^0$
SciPy + FOQS	0.55	0.81	0.19	0.08	1.00	0.98	0.95	0.90	0.89
DFO-GN + Numerical	0.74	0.55	0.24	0.07	1.00	0.98	0.95	0.90	0.89

Table 5. Table 1: Functions of concentration, c 𝑐 {c} . References available in part I, and relevant parameters in Table  2 . (†) α = − ( 2 ​ V ¯ w − V ¯ e ) ​ c max 𝛼 2 subscript ¯ 𝑉 w subscript ¯ 𝑉 e superscript 𝑐 max {\alpha=-(2\bar{V}_{\text{w}}-\bar{V}_{\text{e}})c^{\text{max}}} . (‡) 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}}]} .

Dimensional		Dimensionless
$\hat{D}\left(\hat{c}\right)$	$\left(1.75+2.6\times {10}^{-4}\hat{c}\right)\times {10}^{-9}$	$D\left(c\right)$	$\hat{D}\left(c^{\text{max}}c\right)/\hat{D}\left(c^{\text{max}}\right)$
$\hat{\chi }\left(\hat{c}\right)$	$0.49+4.1\times {10}^{-4}\hat{c}$	$\chi \left(c\right)$	$2\left(1-t_{+}^{\text{w}}\right)\hat{\chi }\left(c^{\text{max}}c\right)/\left(1-\alpha c\right)$ (†)
$\hat{\kappa }\left(\hat{c}\right)$	$\hat{c}\exp (6.23-1.34\times {10}^{-4}\hat{c}$                    $-1.61\times {10}^{-8}{\hat{c}}^2)\times 10{}^{-4}$	$\kappa \left(c\right)$	$RT\hat{\kappa }\left(c^{\text{max}}c\right)/F^2{D}^{\text{max}}{c}^{\text{max}}$
${\hat{c}}_{\text{w}}\left(\hat{c}\right)$	$\left(1-\hat{c}{\overline{V}}_{\text{e}}\right)/{\overline{V}}_{\text{w}}$	$c_{\text{w}}\left(c\right)$	${\hat{c}}_{\text{w}}\left(c^{\text{max}}c\right)/{\hat{c}}_{\text{w}}\left(c^{\text{max}}\right)$
${\hat{j}}_0\left(\hat{c}\right)$	$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_0\left(c\right)$	$𝒜L{\hat{j}}_0\left(c^{\text{max}}c\right)/287𝒞$
${\hat{U}}_{\text{ce}Pb}\left(\hat{c}\right)$	$U_{\text{ce}Pb}^{\ominus }-0.074\log m-0.030{\log }^{2}m$           $-0.031{\log }^{3}m-0.012{\log }^{4}m$ (‡)	$U_{\text{ce}Pb}\left(c\right)$	$\frac{F}{RT}\left({\hat{U}}_{\text{ce}Pb}\left(c^{\text{max}}\hat{c}\right)-U_{\text{ce}Pb}^{\ominus }\right)$
${\hat{U}}_{\text{ce}Pb{O}_2}\left(\hat{c}\right)$	$U_{\text{ce}Pb{O}_2}^{\ominus }+0.074\log m+0.033{\log }^{2}m$           $+0.043{\log }^{3}m+0.022{\log }^{4}m$ (‡)	$U_{\text{ce}Pb{O}_2}\left(c\right)$	$\frac{F}{RT}\left({\hat{U}}_{\text{ce}Pb{O}_2}\left(c^{\text{max}}\hat{c}\right)-U_{\text{ce}Pb{O}_2}^{\ominus }\right)$

Table 6. Table 2: Relevant dimensional parameters from the literature, for Table  1 . Parameters with several values indicate different values in negative electrode (n), separator (sep) and positive electrode (p). References are available in part I.

Parameter	Value			Units
	n	sep	p
$L$	$3.65\times {10}^{-3}$			m
$A_{\text{cs}}$	$7.4\times {10}^{-3}$			m2
$s_{+}$	$-1$	-	$-3$	-
$s_{-}$	$1$	-	$-1$	-
$s_{\text{w}}$	$0$	-	$2$	-
$n_{\text{e}}$	$2$	-	$2$	-
${\overline{V}}_{\text{w}}$	$1.75\times {10}^{-5}$			m3 mol-1
${\overline{V}}_{\text{e}}$	$4.50\times {10}^{-5}$			m3 mol-1
$M_{\text{w}}$	$1.8\times {10}^{-2}$			kg mol-1
$F$	$96485$			C mol-1
$R$	$8.314$			J mol-1 K-1
$T$	$298.15$			K
$t_{+}^{\text{w}}$	$0.72$			-
$j^{\text{ref}}$	$8\times {10}^{-2}$	-	$6\times {10}^{-3}$	A m-2
$c^{\text{max}}$	$5.6\times {10}^3$			mol m-3
$𝒜$	$2.6\times {10}^6$	-	$2.05\times {10}^7$	m-1
$Q$	$17$			Ah
$U^{\ominus }$	-0.295	-	1.628	V

Table 7. Table 1: Parameters obtained from optimisation algorithms/model combinations. All algorithms were given the values in Tables  1 and 2 as initial data for the parameters.

	${\epsilon }_{\text{n},\text{p}}^{\text{max}}$	${\epsilon }_{\text{sep}}^{\text{max}}$	$j_{\text{n}}^{\text{ref}}$	$R_{\text{circuit}}$	$q_{\text{2.5A}}^0$	$q_{\text{2A}}^0$	$q_{\text{1.5A}}^0$	$q_{\text{1A}}^0$	$q_{\text{0.5A}}^0$
SciPy + LOQS	0.55	0.81	0.19	0.08	1.00	0.98	0.95	0.90	0.89
SciPy + FOQS	0.55	0.81	0.19	0.08	1.00	0.98	0.95	0.90	0.89
SciPy + Composite	0.60	0.90	0.07	0.13	1.00	1.00	1.00	0.98	1.00
SciPy + Numerical	0.60	0.90	0.08	0.15	1.00	1.00	1.00	1.00	1.00
DFO-GN + LOQS	0.75	0.53	0.19	0.08	1.00	0.98	0.95	0.90	0.89
DFO-GN + FOQS	0.51	0.88	0.32	0.07	1.00	0.98	0.95	0.90	0.89
DFO-GN + Composite	0.60	0.78	0.12	0.06	1.00	0.98	0.96	0.90	0.89
DFO-GN + Numerical	0.74	0.55	0.24	0.07	1.00	0.98	0.95	0.90	0.89

Equations

$$\frac{\partial }{\partial t}(\epsilon c)$$

$$\frac{\partial \epsilon }{\partial t}$$

$$\frac{\partial i}{\partial x}$$

$${\mathcal{C}}_{\text{d}}i$$

$$\frac{\partial i_{\text{s}}}{\partial x}$$

$$i_{\text{s}}$$

$$j$$

$${\Phi }_{\text{s}}=\frac{\partial c}{\partial x}=i=0,i_{\text{s}}={\mathrm{i}}_{\text{cell}}$$

$$i_{\text{s}}=0$$

$$c$$

$$\epsilon$$

$$\Phi$$

$${\Phi }_{\text{s}}$$

$${\int }_0^{{\ell }_{\text{n}}}{j}_{\text{n}}\mathrm{d}\hat{x}=-{\int }_{1-{\ell }_{\text{p}}}^1{j}_{\text{p}}\mathrm{d}\hat{x}=i_{\text{cell}},$$

$$\frac{\partial {\Phi }_{\text{s}}}{\partial x}=0$$

$${\Phi }_{\text{s},\text{n}}=0,{\Phi }_{\text{s},\text{p}}=V(t).$$

$$\frac{\partial }{\partial t}(\epsilon c)$$

$$\frac{\partial \epsilon }{\partial t}$$

$${\mathcal{C}}_{\text{d}}j$$

$$j_{\text{n}}$$

$$j_{\text{p}}$$

$$\frac{\partial c}{\partial x}=\frac{\partial \Phi }{\partial x}=0\text{ at }x=0,1,$$

$${\int }_0^{{\ell }_{\text{n}}}{j}_{\text{n}}\mathrm{d}x=-{\int }_{1-{\ell }_{\text{p}}}^1{j}_{\text{p}}\mathrm{d}x={\mathrm{i}}_{\text{cell}},$$

$$f(x,t)=f^{(0)}(x,t)+{\mathcal{C}}_{\text{d}}{f}^{(1)}(x,t)+{\mathcal{C}}_{\text{d}}^2{f}^{(2)}(x,t)+\mathcal{O}({\mathcal{C}}_{\text{d}}^3),$$

$$0$$

$$\frac{\partial {\epsilon }^{(0)}}{\partial t}$$

$$0$$

$$j_{\text{n}}^{(0)}$$

$$j_{\text{p}}^{(0)}$$

$$\frac{\partial c^{(0)}}{\partial x}=\frac{\partial {\Phi }^{(0)}}{\partial x}=0\text{ at }x=0,1,$$

$${\int }_0^{{\ell }_{\text{n}}}{j}_{\text{n}}^{(0)}\mathrm{d}x=-{\int }_{1-{\ell }_{\text{p}}}^1{j}_{\text{p}}^{(0)}\mathrm{d}x={\mathrm{i}}_{\text{cell}},$$

$$c^{(0)}=c^0,{\epsilon }^{(0)}={\epsilon }^0.$$

$$\frac{\partial }{\partial t}$$

$$\frac{\partial {\epsilon }^{(1)}}{\partial t}$$

$$j^{(0)}$$

$$j_{\text{n}}^{(1)}$$

$$j_{\text{p}}^{(1)}$$

$${\eta }_{\text{n}}^{(0)}=-\left({\Phi }^{(0)}+U_{\text{ce}Pb}(c^{(0)})\right),$$

$${\eta }_{\text{p}}^{(0)}=V^{(0)}-{\Phi }^{(0)}-U_{\text{ce}Pb{O}_2}\left(c^{(0)}\right),$$

$${\eta }_{\text{n}}^{(1)}=-\left({\Phi }^{(1)}+c^{(1)}{U}_{\text{ce}Pb}^{\prime }\left(c^{(0)}\right)\right),$$

$${\eta }_{\text{p}}^{(1)}=V^{(1)}-{\Phi }^{(1)}-c^{(1)}{U}_{\text{ce}Pb{O}_2}^{\prime }\left(c^{(0)}\right),$$

$$\frac{\partial c^{(1)}}{\partial x}=\frac{\partial {\Phi }^{(1)}}{\partial x}=0\text{ at }x=0,1,$$

$${\int }_0^{{\ell }_{\text{n}}}{j}^{(1)}\mathrm{d}x=-{\int }_{1-{\ell }_{\text{p}}}^1{j}^{(1)}\mathrm{d}x=0,$$

$$c^{(1)}={\epsilon }^{(1)}=0$$

$$j_{\text{n}}^{(0)}={\mathrm{i}}_{\text{cell}}/{\ell }_{\text{n}},j_{\text{p}}^{(0)}=-{\mathrm{i}}_{\text{cell}}/{\ell }_{\text{p}}.$$

$$\frac{\mathrm{d}{c}^{(0)}}{\mathrm{d}t}$$

$$\times [(s_{\text{n}}-s_{\text{p}}){\mathrm{i}}_{\text{cell}}$$

$$-c^{(0)}\frac{\mathrm{d}}{\mathrm{d}t}\left({\ell }_{\text{n}}{\epsilon }_{\text{n}}^{(0)}+{\ell }_{\text{sep}}{\epsilon }_{\text{sep}}^{\text{max}}+{\ell }_{\text{p}}{\epsilon }_{\text{p}}^{(0)}\right)],$$

$$\frac{\mathrm{d}{\epsilon }_{\text{n}}^{(0)}}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\epsilon }_{\text{p}}^{(0)}}{\mathrm{d}t}$$

$${\mathrm{i}}_{\text{cell}}/{\ell }_{\text{n}}$$

Full text

1 Introduction

The popular equivalent-circuit approach to battery modelling [1] is efficient, but has limited physical detail and extrapolates poorly. Electrochemical models [2, 3, 4, 5, 6, 7, 8, 9, 10] require far more computational power, but include detailed descriptions of physical mechanisms, which presumably enhances predictive capability. Battery management could be improved if there existed easily-solved models with greater mechanistic detail. To that end, this paper puts forward several reduced-order models of lead-acid battery discharge, each derived from a mechanistic description based on an extension of Newman’s porous-electrode theory [11], which we developed in part I.

Several authors have simplified mechanistic lead-acid-battery models to improve their computational efficiency. Newman and Tiedemann [7] recognise that spatial gradients can be ignored at low current; they state a ‘lumped parameter model’ (LPM) that depends only on time, but do not show how it derives from a porous-electrode model. Gandhi et al. [12] propose a LPM to underpin an analytical current/voltage relation. Knauff [13] simplifies a porous-electrode model by assuming, without justification, that current is linear in space, and acid molarity, quadratic.

We deploy perturbation methods [14] to produce a hierarchy of increasingly complex models. After nondimensionalization, a diffusional C-rate, $\mathcal{C}_{\text{d}}$–the C-rate scaled with the diffusion time-scale—is found to control how simply the full model can be approximated. Three reduced-order models are derived, validated against the full model, and applied to experiments for parameter estimation.

A leading-order expansion in the diffusional C-rate produces a LPM of the Newman–Tiedemann type, found to be accurate for C-rates below $0.1C$. The first-order expansion accounts for quasi-static spatial heterogeneity within the electrode sandwich. As well as improving the fit of the full model, this correction has a computationally efficient closed-form expression. Finally, the first-order solution is improved by accounting for diffusion transients. This composite model includes just one linear partial differential equation, but matches the full model well up to $5C$.

2 Dimensionless model

In part I, we proposed a general three-dimensional, thermodynamically consistent, isothermal porous-electrode model of a discharging lead-acid battery. The detailed model was simplified slightly on the basis of dimensional analysis to allow solution in a one-dimensional setting.

After nondimensionalization, we obtained the following dimensionless system governing the electrolyte concentration $c$, porosity $\varepsilon$, current density ${i}$ and potential $\Phi$, electrode current density ${i}_{\text{s}}$ and potential $\Phi_{\text{s}}$, and interfacial current density $j$:

[TABLE]

where property values in the negative and positive electrode are designated with subscripts n and p, respectively. Typical values of the dimensionless parameters $\mathcal{C}_{\text{d}}$, $\iota_{\text{s}}$, $\beta^{\text{surf}}$, $\gamma_{\text{dl}}$, $\ell$, $s$, $q^{0}$, $\varepsilon^{\text{max}}$ and $\varepsilon^{\Delta}$ are given in Table 1, while concentration-dependent functions $D$, $\kappa$, $\chi$, $j_{0}$ and $U$ are given in Table 1. The dimensionless applied current is $i_{\text{cell}}(t)=I_{\text{circuit}}(t)/8A_{\text{cs}}$, where $I_{\text{circuit}}(t)$ is the applied current in the external circuit and $A_{\text{cs}}$ is the electrode cross-sectional area. We define $\bar{i}$ to be the maximum value of $i_{\text{cell}}(t)$ with respect to time.

The key parameter is the diffusional C-rate, $\mathcal{C}_{\text{d}}$, which is the C-rate as measured on the diffusion time-scale (or alternatively, the ratio of the applied current scale to the scale of the limiting current).

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

3 Solutions

We now derive three analytical, approximate solutions to the model system (2.1), and compare these to the numerical solution of the full model computed in part I, which we treat as ‘ground truth’. To do this, we note that the diffusional C-rate, $\mathcal{C}_{\text{d}}$, is small for most practical (low C-rate) applications, and perform an asymptotic analysis near the limit of small $\mathcal{C}_{\text{d}}$.

3.1 Leading-order quasi-static solution

In this section, we will derive the quasi-static solution in the limit of small $\mathcal{C}_{\text{d}}$, $\gamma_{\text{dl}}$ and $1/\iota_{\text{s}}$. Since $\gamma_{\text{dl}}$ and $1/\iota_{\text{s}}$ are much smaller than one (Table 1), we only take the leading-order terms in their expansions. In contrast, $\mathcal{C}_{\text{d}}$ can sometimes be close to one, so we will consider both the leading order and first order in $\mathcal{C}_{\text{d}}$.

To leading order in $1/\iota_{\text{s}}$, (2.1f) becomes

[TABLE]

in each electrode, and so $\Phi_{\text{s}}$ can be approximated as a function of time only. Applying the boundary condition (2.1h) and defining $V(t)=\left.\Phi_{\text{s}}\right\rvert_{x=1}$, we can now replace (2.1f) with

[TABLE]

We use the integral condition (2.1m) so that we do not need to solve for $i_{\text{s}}$ to find the voltage, $V(t)$. Hence (2.1f) is only necessary if we want to find $i_{\text{s}}$ having found $j$. We also take the leading order in $\gamma_{\text{dl}}$, so that the time derivatives in (2.1g) disappear.

In summary, we simplify the system (2.1) to the following equations for $c(x,t)$, $\varepsilon(x,t)$, $j(x,t)$, $\Phi(x,t)$ and $V(t)$:

[TABLE]

As shown in Table 1, the diffusional C-rate, $\mathcal{C}_{\text{d}}$, is equal to $0.6\mathcal{C}$, where $\mathcal{C}$ is the C-rate. Most practical applications have a C-rate below 0.25C, so the diffusional C-rate is usually small. Hence we perform an asymptotic expansion in the limit $\mathcal{C}_{\text{d}}\to 0$ and assume that we can expand all variables in (2.1c) in powers of $\mathcal{C}_{\text{d}}$:

[TABLE]

where $f=c,\varepsilon,\Phi,j$ and $V$. Hence (2.1c) becomes to leading order

[TABLE]

At first order, equating coefficients of $\mathcal{C}_{\text{d}}$ in (2.1c) gives

[TABLE]

Leading-order quasi-static solution

We now seek the solution to the lowest order problem. Integrating (2.1ea) with boundary conditions (2.1ef), then integrating again, gives $c^{(0)}=c^{(0)}(t)$. We then integrate (2.1ec), use boundary conditions (2.1ef), and integrate again, to find that $\Phi^{(0)}=\Phi^{(0)}(t)$. Hence $j^{(0)}_{\text{n}}$ and $j^{(0)}_{\text{p}}$ as defined by (2.1ed) and (2.1ee) are functions of time only; the boundary conditions (2.1eg) give

[TABLE]

Finally, by (2.1eb), $\varepsilon_{\text{n}}^{(0)}$, $\varepsilon_{\text{sep}}^{(0)}$ and $\varepsilon_{\text{p}}^{(0)}$ are functions of time only (in particular, $\varepsilon_{\text{sep}}^{(0)}\equiv\varepsilon_{\text{sep}}^{\text{max}}$). Hence to leading order, the whole problem is quasi-static. To determine $c^{(0)}$, we need to consider the first-order problem (2.1fa) for $c^{(1)}$. Integrating (2.1fa) from $x=0$ to $x=1$ and using (2.1g) and the boundary conditions (2.1fj) gives a solvability condition that determines $c^{(0)}$. We can combine this with (2.1eb), (2.1ed) and (2.1ee) to obtain a nonlinear differential-algebraic equation system governing $c^{(0)}$, $\varepsilon_{\text{n}}^{(0)}$, $\varepsilon_{\text{p}}^{(0)}$, $\Phi^{(0)}$ and $V^{(0)}$,

[TABLE]

with initial conditions (2.1eh). Integrate (2.1ha-c) and rearrange (2.1hd,e) to find the final leading-order solution,

[TABLE]

3.2 First-order quasi-static solution

We now solve the first-order system, (2.1f), to find the $\mathcal{O}(\mathcal{C}_{\text{d}})$ correction to the voltage. We solve (2.1f) as follows:

(i) find $c^{(1)}$ using (2.1fa), up to an arbitrary constant, $k(t)$;

(ii) find $k$ using a solvability condition on $c^{(2)}$, the $\mathcal{O}(\mathcal{C}_{\text{d}}^{2})$ correction to $c$;

(iii) find $\Phi^{(1)}$ using (2.1fc) up to an arbitrary constant, $A_{\text{n}}(t)$;

(iv) find $A_{\text{n}}$ using (2.1fd);

(v) find $V^{(1)}$ using (2.1fe).

Firstly, with known $c^{(0)}$ and $\varepsilon^{(0)}$, we can integrate (2.1fa) with respect to $x$ twice and use (2.1fj) to find an explicit equation for $c^{(1)}$ (given in B).

Having found $c^{(1)}$, we integrate (2.1fc), using (2.1g) and continuity of $\Phi^{(1)}$, to find

[TABLE]

where $A_{\text{n}}$ is an arbitrary constant.

We can now integrate (2.1fd) from $x=0$ to $x=\ell_{\text{n}}$ and integrate (2.1fe) from $x=1-\ell_{\text{p}}$ to $x=1$, using (2.1fk) each time, to find the correction term $V^{(1)}$:

[TABLE]

where we have introduced the averages

[TABLE]

3.3 Composite solution

The quasi-static solution developed in the ‘Leading-order quasi-static solution’ and ‘First-order quasi-static solution’ sections is valid when the current varies slowly, but fails to capture transient behaviour when the current changes more rapidly, such as a jump. To capture such transients, we could rescale time with $\tau=(t-t^{*})/\mathcal{C}_{\text{d}}$, where $t^{*}$ is the time of the jump in the current, define $C(\tau)=c(t)$ (and likewise for other variables) and expand in powers of $\mathcal{C}_{\text{d}}$. We give the details of such an approach in C.

Such a transient solution is valid at short times after a jump time $t^{*}$, but breaks down at times long after the jump time. To obtain a solution that is valid both at short times after a jump in current and at long times, without having to repeatedly ‘reset’ the transient solution, we use a ‘composite’ solution, which we now develop here.

We consider the lowest order and first order correction for the concentration by taking $\tilde{c}=c^{(0)}+\mathcal{C}_{\text{d}}\,c^{(1)}$. We then consider the PDE

[TABLE]

where $c^{(0)}$ and $\varepsilon^{(0)}$ are given by the quasi-static problem (2.1i) and $j^{(0)}$ is given by (2.1g). We note that for long times, $\mathcal{C}_{\text{d}}\,\partial c^{(1)}/\partial t$ is a higher-order term and we retrieve the quasi-static problem (2.1fa), while for short times, re-scaling $\tau=(t-t^{*})/\mathcal{C}_{\text{d}}$, $c^{(0)}$ is constant and we have the transient problem (2.1ca) for $c^{(1)}$. Hence (2.1m) is valid uniformly at both short times and long times.

The composite solution then consists of solving (2.1m) for $\tilde{c}$, then computing

[TABLE]

and finally finding $V^{(1)}$ through (2.1j) and (2.1k) with $c^{(1)}$ given by (2.1n).

4 Results

In the Solutions section,we derived four systems that are approximately equivalent to the full dimensionless system (2.1) with varying degrees of accuracy:

Numerical – part I 2. 2.

Leading-order quasi-static (LOQS) – (2.1i) 3. 3.

First-order quasi-static (FOQS) – (2.1a), (2.1j) and (2.1k) 4. 4.

Composite – (2.1m), (2.1j) and (2.1k)

The code used to solve the models and generate the results below is available publicly on GitHub [15]. Note that to obtain either the first-order quasi-static solution or the composite solution, we must first solve the leading-order quasi-static problem.

4.1 Reduced-order solutions

We now compare results from the four models. We treat the full numerical model as ‘ground truth’, and investigate the speed and accuracy of the three other models compared to the numerical model.

The most important output from the model is the voltage, since this is the variable that we can compare to experimental data (treating current as a known input). In Figure 1, we compare the voltage during a complete constant-current discharge at a range of C-rates. The discharge is deemed to be finished either when the concentration reaches zero anywhere in the cell, or when the voltage reaches a cut-off voltage of $10.5$V.

We observe that all three reduced-order solutions agree well with the numerical solution at very low C-rates (Figure 1a). As we increase the C-rate (Figures1b-d), only the first-order solutions (FOQS and composite) agree with the numerical solution; further, a discrepancy appears between the FOQS solution and the numerical solution at early times. Finally, for very high C-rates (Figures 1e-f) the composite solution still agrees very well with the numerical solution, but the FOQS solution does not, and terminates early, for reasons that we explain below.

To explain the behaviour observed in the voltages, we investigate internal variables, such as the concentration at various states of charge (Figure 2). At a very low C-rate of 0.1C (Figure 2a), the concentration remains almost uniform throughout the discharge; hence the LOQS solution, which does not take into account any spatial variations, provides a good fit to the numerical solution. At a higher C-rate of 0.5C (Figure 2b), the concentration in the numerical solution is no longer spatially homogeneous; this is non-uniformity is captured well by the FOQS and composite solutions, but not by the LOQS solution. However, even with the FOQS and composite solutions, there is a discrepancy in the concentration profiles in the positive electrode (Figures 2b,c, right-hand side of the spatial domain). This is because the solutions from the asymptotic methods assume a uniform interfacial current density, but in the numerical solution the interfacial current density is non-uniform.

Finally, at high C-rates (2C, Figure 2c), there is a diffusion transient at the start of the discharge; this is only captured by the composite solution, and not the FOQS solution. This initial diffusion transient also explains the discrepancy between the FOQS and numerical solutions at early times in Figure 1d. In addition to this, we can now see that the FOQS solution terminates early in Figure 1f because the concentration quickly reaches zero.

In Figure 3, we show the relative errors of the voltage obtained from reduced-order models compared to the voltage obtained from the numerical model.

Then, in Table 2, we compare the time taken to solve the various models. We see that the composite solution, leading-order quasi-static and first-order quasi-static solutions are roughly one, two and three orders of magnitude faster than the full numerical solution respectively. Coupled with the errors shown in Figure 3, the speeds shown in Table 2 suggest that in order to solve the model accurately and as quickly as possible, we should use the LOQS model for C-rates below 0.1C, the FOQS model for C-rates of 0.1-1C, and the composite model for C-rates above 1C.

The time taken for the leading-order quasi-static model is independent of grid size, while the time taken for the other models scales linearly with grid size. Note that we can expect to obtain a faster numerical solution by using a different spatial discretisation scheme than Finite Volumes, such as Chebyshev orthogonal collocation [16], and the relative speed-up of the composite solution by using the same discretisation would be similar.

4.2 Voltage breakdown

As well as obtaining a faster solution to the model, the composite solution allows us to identify the individual overpotentials that contribute to the total drop in voltage from full charge. We write the total dimensional voltage as the sum of the initial voltage,

[TABLE]

and individual voltage drops [17],

[TABLE]

where $V_{\text{U},i}$, $i=\text{n},\text{p}$ are the open-circuit voltages; $V_{\text{k},i}$, ${i}=\text{n},\text{p}$ are the kinetic overpotentials, accounting for losses due to the reactions at the electrode-electrolyte interfaces; $V_{\text{c}}$ is the concentration overpotential, accounting for losses due to concentration gradients; and $V_{\text{o}}$ is the Ohmic overpotential in the electrolyte, accounting for losses due to the electric resistance of the electrolyte. Equation 2.1h would usually include a term to account for Ohmic losses in the solid electrodes, but in our reduced-order models this term is zero since $\iota_{\text{s}}$ is large (c.f. equation 2.1a).

Combining (2.1ie), (2.1j) and (2.1k), we identify

[TABLE]

Together with the quasi-static formulas (2.1ia) for $c^{(0)}$ and (2.1a) for $c^{(1)}$, equations (2.1h) and (2.1ia) give an exact formula for the voltage that is valid for most operating C-rates (below 0.5C). For higher C-rates, we must solve (2.1m) and use (2.1n), instead of (2.1a), to find $c^{(1)}$.

In Figure 4, we show the relative contribution from each of the terms in (2.1ia). At low C-rates (Figure 4a), the drop in voltage is almost entirely due to the change in the OCV of the two electrodes as the concentration changes, $V_{\text{U,n}}$; this is why the LOQS solution is accurate enough in this regime. As we increase the C-rate (Figure 1b), the first non-OCV effects to become important are the kinetic overpotentials, $V_{\text{k,n}}$, with the kinetic overpotential in the positive electrode being greater. Finally, at very high C-rates (Figure 1c), we see that the sharp initial drop in voltage is due to both kinetic overpotentials and the ohmic overpotential in the electrolyte in roughly equal measures.

The breakdown of voltages given by (2.1h) could be used to help guide any optimisation of a battery design. For example, the behaviour observed in Figure 4b suggests that an effective way to improve the power capacity of the battery (by reducing the voltage drop) might be to target reduction of the kinetic overpotential in the positive electrode. Alternatively, this allows us to explore the trade-off between energy capacity and power capacity of the battery, and hence optimise the behaviour of the battery for specific applications. For example, increasing the total width, $L$, of an electrode pair would increase the energy capacity of the battery. However, increasing $L$ would also increase the diffusional C-rate (Table 1), and so increase the contributions to the voltage drop in (2.1ia), decreasing the power capacity of the battery. It becomes evident that we need to consider first-order effects in order to fully understand the trade-off between energy capacity and power capacity; the LOQS model would naively suggest that increasing $L$ is always the optimal strategy. A third option could be to fix the total width $L$ (i.e. fix the energy capacity) and optimise the relative widths of the negative electrode, separator and positive electrode in order to minimise the voltage drop given by (2.1h), and hence maximise the power capacity of the battery. In all three examples suggested, a quantitative optimisation could be done very efficiently using the formulas given in this paper.

4.3 Parameter fitting

Another potential application for our reduced-order models is in parameter fitting. This is an important challenge when comparing models to experimental data, as many parameters in the model cannot be identified individually a priori. To demonstrate how our model can be useful for parameter fitting, we fit each model to six constant-current discharges of a 17 Ah BBOXX Solar Home battery at intervals of 0.5 A from 3 A to 0.5 A. Each constant-current discharge is followed by a two-hour rest period during which the current is zero.

In order to use as few fitting parameters as possible, we take standard values for all parameters (given in Tables 1 and 2) except for the following:

1. (i)

the maximum electrode porosities, $\varepsilon_{\text{n}}^{\text{max}}$ and $\varepsilon_{\text{p}}^{\text{max}}$, which we assume to be equal; 2. (ii)

the separator porosity, $\varepsilon_{\text{sep}}^{\text{max}}$; 3. (iii)

the exchange current-density for the negative electrode, $j_{\text{n}}^{\text{ref}}$ (we assume that ${j_{\text{p}}^{\text{ref}}=j_{\text{n}}^{\text{ref}}/10}$); 4. (iv)

a ‘resistance in the wires’, $R_{\text{circuit}}$, which accounts for ohmic losses outside the battery within the BBOXX Solar Home (we take away $I_{\text{circuit}}R_{\text{circuit}}$ from the voltage output by the model before fitting to data). 5. (v)

an initial SOC, $q^{0}$, for five of the curves (we assume that the first discharge, at 3A, starts from full SOC);

This gives a total of nine fitting parameters for six curves.

To fit our model to the data, we minimise the sum-of-squares of the voltage prediction error. We do this in Python with both a derivative-based (leastsquares [18]) or derivative-free (DFO-GN [19]) optimisation algorithm.

Since the highest discharge current is quite low (3A corresponds to a C-rate of 0.18C), both the first-order quasi-static model and composite model agree well with the full numerical model. These three models give a very good fit to the data except near the start of the discharge, and during the short rest period afterwards, as shown in Figure 5. In Table 3, we compare the final sum-of-squares error, as well as the time taken to find the minimum when starting from the parameters from literature of Tables 1 and 2, for each optimisation algorithm/model combination. As expected, the quasi-static models are faster than the composite model, with the numerical model being by far the slowest. However, this is compounded by the fact that the derivative-based algorithm from SciPy cannot find the best fit for the composite and numerical models, and so we must resort to derivative-free models. Thus optimising the first-order quasi-static model is thirty times faster than optimising the full numerical model.

As shown in Table 4, the optimal parameters found by SciPy with the FOQS model are almost the same as those found by DFO-GN with the numerical model, with the exception of the exchange current-density, $j_{\text{n}}^{\text{ref}}$. Hence one could first fit all the parameters using SciPy and the FOQS model, and then only need to fit $j_{\text{n}}^{\text{ref}}$ with DFO-GN and the numerical model.

The discrepancy at early times may be due to effects from charging: we assume that the discharge starts from equilibrium, but this might not be the case following a charge. The final relaxation timescale is much longer than the two timescales in our model – the diffusion timescale ($\sim 15$ minutes) and the capacitance timescale ($\sim 5$ seconds) – which suggests that there is an extra physical effect that we have not considered. In particular, non-uniformities in the $y$- and $z$-dimensions may be important for both early time and relaxation behaviour.

5 Conclusions

The asymptotic methods developed in this paper allow us to simulate a discharge of a lead-acid battery with the low complexity and high speed of equivalent-circuit models, while retaining the accuracy and physical insights of electrochemical models. Further, these methods give important physical insight into the structure of the problem that is not obvious from numerical solutions of the full model. In particular, we observe that at low C-rates, the voltage drop from open-circuit potential is mainly due to kinetic effects; ohmic and concentration overpotentials are relatively very small. It is also important to note that our models are tunable: for discharges at low C-rate, we should choose the leading-order quasi-static model, while for high C-rates we can use the composite model, which still provides a significant speed-up compared to the full model.

There are many exciting applications for the models resulting from our asymptotic methods. Firstly, these models can be used in Battery Management Systems to replace equivalent-circuit models without introducing the complexity of full electrochemical models. We expect to find significant advantages from doing this, as the parameters will be much more robust to different states of charge than less physically based resistances and capacitances. Secondly, as demonstrated in this paper, the asymptotic models can be used to estimate parameter values, and compare results to experimental data, much more quickly than using the full electrochemical models; we have found that models are often not compared to experimental data due to the prohibitively large cost of the full electrochemical models. Thirdly, since the voltage is given in an explicit form, we can easily perform a parameter sensitivity analysis. Fourthly, by analysing the individual voltage drops from open-circuit potential, we can optimise battery design for specific applications, for example trading off between power and capacity.

An important extension for this work is to apply asymptotic methods to more complex models, for example including side reactions occurring during overcharge, or other long-term degradation mechanisms for lead-acid batteries such as corrosion and irreversible sulfation. We expect the ideas developed here to continue to apply in these more complex settings.

Finally, the ideas explored in this paper can be applied to other battery chemistries. In particular, our leading-order quasi-static model is very similar to the Single-Particle Model for lithium-ion batteries [20, 21, 22, 23, 24], and our composite model is very similar to the Single-Particle Model with electrolyte  [25, 26, 27, 28, 29, 30]. However, different chemistries and geometries may lead to different parameter sizes, and hence different distinguished limits.

Acknowledgements

This publication is based on work supported by the EPSRC Centre 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

$(0)$

leading-order

$(1)$

first-order

eff

effective

surf

surface

$\square$

convective

Accents

$\hat{}$

dimensional

$\bar{}$

averaged

Appendix A Parameters

The concentration-dependent functions are given in Table 1.

Appendix B Concentration in the first-order quasi-static solution

[TABLE]

We note that the piece-wise quadratic form of (2.1a) justifies the assumption of Knauff [13]. To find $k(t)$, we expand (2.1a) to second-order in $\mathcal{C}_{\text{d}}$:

[TABLE]

Integrating (2.1d) from $x=0$ to $x=1$ and using the fact that ${\partial c^{(1)}/\partial x=\partial c^{(2)}/\partial x=0}$ at $x=0,1$, together with (2.1fb) and (2.1fk), we find that

[TABLE]

and so, by (2.1fl), $\int_{0}^{1}\!\varepsilon^{(0)}c^{(1)}\,\mathrm{d}x\equiv 0$. Hence

[TABLE]

Appendix C Transient solution

Denoting the time of the jump by $t^{*}$, we rescale time with $\tau=(t-t^{*})/\mathcal{C}_{\text{d}}$, and define $C(x,\tau)=c(x,t)$, $E(x,\tau)=\varepsilon(x,\tau)$, $J(x,\tau)=j(x,t)$, $P(x,\tau)=\Phi(x,t)$ and $\mathcal{V}(\tau)=V(t)$. Then (2.1c) becomes

[TABLE]

We now expand the variables in powers of $\mathcal{C}_{\text{d}}$, as done in (2.1d). Then (2.1a) becomes, to leading order in $\mathcal{C}_{\text{d}}$,

[TABLE]

and to first order in $\mathcal{C}_{\text{d}}$,

[TABLE]

The boundary conditions follow from (2.1ef,g) and (2.1fe,f):

[TABLE]

Equations (2.1b) decouple: we first solve (2.1bb) for $E^{(0)}$, then (2.1ba) with (2.1da) for $C^{(0)}$, then (2.1bc) for $P^{(0)}$ up to an additive constant $A$, and finally use (2.1bd) and (2.1db) to find $A$ and (2.1be) and (2.1dc) to find $\mathcal{V}^{(0)}$. Similarly, having found the leading-order solution, equations (2.1c) decouple and can be solved in the same order as (2.1b). Hence we solve (2.1b) and (2.1c) more easily than the full system (2.1a).

However, we note that (2.1b), (2.1c) and the boundary conditions (2.1d) imply that $\int_{0}^{1}\!C^{(0)}\,\mathrm{d}x$ is constant, but $\rvert C^{(1)}\rvert$ grows with time. Hence the asymptotic expansion breaks down for long times since, when $\tau\sim\mathcal{O}(1/\mathcal{C}_{\text{d}})$, \rvert\text{\mathcal{C}_{\text{d}}}\,C^{(1)}\rvert\sim\rvert C^{(0)}\rvert.

Appendix D Parameters from fit to data

The full list of parameters obtained from each combination of fitting algorithm (SciPy or DFOGN) and model (LOQS, FOQS, Composite or Numerical) is given in Table 1.