An asymptotic derivation of a single particle model with electrolyte

TL;DR

This paper systematically derives an improved single particle model with electrolyte effects for lithium-ion batteries, quantifying errors and highlighting neglected terms in previous ad-hoc models to enhance accuracy without increasing complexity.

Contribution

A rigorous asymptotic derivation of the SPM with electrolyte effects, identifying key neglected terms and providing error bounds for different operating conditions.

Findings

Quantifies the error of the reduced model based on parameter ratios.

Shows previous models neglect important terms affecting accuracy.

Maintains computational complexity while improving model accuracy.

Abstract

The standard continuum model of a lithium-ion battery, the Doyle-Fuller-Newman (DFN) model, is computationally expensive to solve. Typically simpler models, such as the single particle model (SPM), are used to provide insight for control purposes. Recently, there has been a move to extend the SPM to include electrolyte effects, which increase the accuracy and range of applicability. However, these extended models are derived in an ad-hoc manner, which leaves open the possibility that important terms may have been neglected, resulting in the model not being as accurate as possible. In this paper, we provide a systematic asymptotic derivation of both the SPM and a correction term that accounts for the behaviour in the electrolyte. Firstly, this allows us to quantify the error in the reduced model in terms of ratios of key parameters in the model, from which the range of applicable operating conditions can be determined. Secondly, in comparing our model with the ad-hoc models from the literature, we show that previous models have neglected a key set of terms. In particular, we make the crucial distinction between writing the terminal voltage in pointwise and electrode-averaged form, which allows us to gain additional accuracy whilst maintaining the same computational complexity as the existing models.

Figures

Tables

Table 1. Table 1 : Dimensional model parameters with values taken from [ 27 ] .

Parameter	Units	Description	${\Omega }_{\text{n}}^{\ast }$	${\Omega }_{\text{s}}^{\ast }$	${\Omega }_{\text{p}}^{\ast }$
${ϵ}_{\text{k}}$	-	Electrolyte volume fraction	$0.3$	1	$0.3$
$c_{\text{k,max}}^{\ast }$	$\mathrm{mol}/m^3$	Maximum lithium concentration	$2.4983\times {10}^4$	-	$5.1218\times {10}^4$
${\sigma }_{\text{k}}^{\ast }$	$S/m$	Solid conductivity	100	-	10
$D_{\text{s,k}}^{\ast }$	$m^2/s$	Electrode diffusivity	$3.9\times {10}^{-14}$	-	$1\times {10}^{-13}$
$R_{\text{k}}^{\ast }$	$\mu m$	Particle radius	10	-	10
$a_{\text{k}}^{\ast }$	$\mu m^{-1}$	Electrode surface area density	0.18	-	0.15
$m_{\text{k}}^{\ast }$	$\left(A/m^2\right){\left(m^3/\mathrm{mol}\right)}^{1.5}$	Reaction rate	$2\times {10}^{-5}$	-	$6\times {10}^{-7}$
$L_{\text{k}}^{\ast }$	$\mu m$	Thickness	100	25	100
$U_{\text{k,ref}}^{\ast }$	$V$	Reference OCP	0.18	-	3.94
$c_{\text{e, typ}}^{\ast }$	$\mathrm{mol}/m^3$	Typical lithium-ion concentration in electrolyte		$1\times {10}^3$
$D_{\text{e,typ}}^{\ast }$	$m^2/s$	Typical electrolyte diffusivity		$5.34\times {10}^{-10}$
${\kappa }_{\text{e,typ}}^{\ast }$	$S/m$	Typical electrolyte conductivity		$1.1$
$F^{\ast }$	$C/\mathrm{mol}$	Faraday’s constant		96485
$R^{\ast }$	$J/\left(\mathrm{mol}K\right)$	Universal gas constant		8.314472
$T^{\ast }$	$K$	Temperature		298.15
$b$	-	Bruggeman coefficient		1.5
$t^{+}$	-	Transference number		0.4
$I_{\text{typ}}^{\ast }$	$A/m^2$	Typical current density		24 ($1C$)

Table 2. Table 2 : Timescales associated with the physical processes occurring within the battery model. Here, 𝒞 𝒞 \mathcal{C} is the C-rate.

Symbol	Expression	Interpretation	Value $\left[s\right]$
${\tau }_{\text{d}}^{\ast }$	$F^{\ast }{c}_{\text{n,max}}^{\ast }{L}^{\ast }/I_{\text{typ}}^{\ast }$	Discharge timescale	$2.2598\times {10}^4/𝒞$
${\tau }_{\text{n}}^{\ast }$	${\left(R_{\text{n}}^{\ast }\right)}^2/D_{\text{s,p}}^{\ast }$	Diffusion timescale in the negative electrode particle	$2.5641\times {10}^3$
${\tau }_{\text{p}}^{\ast }$	${\left(R_{\text{p}}^{\ast }\right)}^2/D_{\text{s,p}}^{\ast }$	Diffusion timescale in the positive electrode particle	$1\times {10}^3$
${\tau }_{\text{e}}^{\ast }$	${\left(L^{\ast }\right)}^2/D_{\text{e,typ}}^{\ast }$	Diffusion timescale in the electrolyte	94.803
${\tau }_{\text{r,n}}^{\ast }$	$F^{\ast }/\left(m_{\text{n}}^{\ast }{a}_{\text{n}}^{\ast }{\left(c_{\text{e,typ}}^{\ast }\right)}^{1/2}\right)$	Reaction timescale in the negative electrode	847.534
${\tau }_{\text{r,p}}^{\ast }$	$F^{\ast }/\left(m_{\text{p}}^{\ast }{a}_{\text{p}}^{\ast }{\left(c_{\text{e,typ}}^{\ast }\right)}^{1/2}\right)$	Reaction timescale in the positive electrode	$3.3901\times {10}^4$

Table 3. Table 3 : Typical dimensionless parameter values. Here 𝒞 = I ∗ / ( 24   Am − 2 ) 𝒞 superscript 𝐼 times 24 superscript Am 2 \mathcal{C}=I^{*}/($24\text{\,}\mathrm{A}\mathrm{m}^{-2}$) is the C-Rate where we have taken a 1C rate to correspond to a typical x 𝑥 x -direction current density of 24   Am − 2 times 24 superscript Am 2 24\text{\,}\mathrm{A}\mathrm{m}^{-2} . This is for a cell with an initial stoichiometry of 0.8 0.8 0.8 in the negative electrode and 0.6 0.6 0.6 in the positive electrode with a voltage cutoff of 3.2   V times 3.2 volt 3.2\text{\,}\mathrm{V} .

Parameter	Expression	Interpretation	n	s	p
$L_{\text{k}}$	$L_{\text{k}}^{\ast }/L^{\ast }$	Ratio of region thickness to cell thickness	0.4444	0.1111	0.4444
${𝒞}_{\text{k}}$	${\tau }_{\text{k}}^{\ast }/{\tau }_{\text{d}}^{\ast }$	Ratio of solid diffusion and discharge timescales	0.1134 $𝒞$	-	0.0442 $𝒞$
${𝒞}_{\text{r,k}}$	${\tau }_{\text{r,k}}^{\ast }/{\tau }_{\text{d}}^{\ast }$	Ratio of reaction and discharge timescales	0.0375 $𝒞$	-	1.5 $𝒞$
${\sigma }_{\text{k}}$	$\left(R^{\ast }{T}^{\ast }/F^{\ast }\right)/\left(\left(I_{\text{typ}}^{\ast }{L}^{\ast }/{\sigma }_{\text{k}}^{\ast }\right)\right)$	Ratio of thermal voltage and typical Ohmic drop in the solid	$475.791/𝒞$	-	$47.5791/𝒞$
$a_{\text{k}}$	$a_{\text{k}}^{\ast }{R}_{\text{k}}^{\ast }$	Product of particle radius and surface area density	1.8	-	1.5
${\gamma }_{\text{k}}$	$c_{\text{k,max}}^{\ast }/c_{\text{n,max}}^{\ast }$	Ratio of maximum lithium concentrations in solid	1	-	2.0501
${𝒞}_{\text{e}}$	${\tau }_{\text{e}}^{\ast }/{\tau }_{\text{d}}^{\ast }$	Ratio of electrolyte transport and discharge timescales		$4.19\times {10}^{-3}𝒞$
${\gamma }_{\text{e}}$	$c_{\text{e,typ}}^{\ast }/c_{\text{n,max}}^{\ast }$	Ratio of maximum lithium concentration in the negative electrode solid and typical electrolyte concentration		$0.04$
${\hat{\kappa }}_{\text{e}}$	$\left(R^{\ast }{T}^{\ast }/F^{\ast }\right)/\left(\left(I_{\text{typ}}^{\ast }{L}^{\ast }/{\kappa }_{\text{e,typ}}^{\ast }\right)\right)$	Ratio of thermal voltage to the typical Ohmic drop in the electrolyte		$5.2337/𝒞$
$c_{\text{k,0}}$	$c_{\text{k,0}}^{\ast }/c_{\text{k,max}}^{\ast }$	Ratio of initial lithium concentration to maximum lithium concentration in solid	0.8	-	0.6

Table 4. Table 4 : RMS voltage error between the reduced models and the DFN model for the finite-volume implementation with 30 points in each domain, as discussed in Section  5 .

	$0.1C$	$0.5C$	$1C$	$2C$	$3C$
SPM	$1.72\mathrm{mV}$	$9.62\mathrm{mV}$	$19.86\mathrm{mV}$	$40.67\mathrm{mV}$	$62.78\mathrm{mV}$
SPMe	$0.17\mathrm{mV}$	$1.34\mathrm{mV}$	$3.04\mathrm{mV}$	$7.36\mathrm{mV}$	$13.34\mathrm{mV}$

Table 5. Table 5 : Average computation time across C-rates (0.1 C, 0.5 C, 1 C, 2 C, 3 C) for the SPM ( 20 ), SPMe ( 40 ), and DFN model ( 1 ). Here,‘ n 𝑛 n  pts’ refers to n 𝑛 n points in each domain.

	5 pts	10 pts	20 pts	30 pts
SPM	0.24 s	0.25 s	0.34 s	0.41 s
SPMe	0.50 s	0.58 s	0.56 s	0.59 s
DFN	3.34 s	5.10 s	18.51 s	61.36 s

Table 6. Table 6 : The key conditions to be satisfied for the application of ( 48 ). In addition, it is required that 𝒞 e ≪ L k ∗ / L ∗ ≪ 1 / 𝒞 e much-less-than subscript 𝒞 e superscript subscript 𝐿 k superscript 𝐿 much-less-than 1 subscript 𝒞 e \mathcal{C}_{\text{e}}\ll L_{\text{k}}^{*}/L^{*}\ll 1/\mathcal{C}_{\text{e}} , 𝒞 e ≪ c p,max / c n,max ≪ 1 / 𝒞 e much-less-than subscript 𝒞 e subscript 𝑐 p,max subscript 𝑐 n,max much-less-than 1 subscript 𝒞 e \mathcal{C}_{\text{e}}\ll c_{\text{p,max}}/c_{\text{n,max}}\ll 1/\mathcal{C}_{\text{e}} , and 𝒞 e ≪ c e,typ / c n,max ≪ 1 / 𝒞 e much-less-than subscript 𝒞 e subscript 𝑐 e,typ subscript 𝑐 n,max much-less-than 1 subscript 𝒞 e \mathcal{C}_{\text{e}}\ll c_{\text{e,typ}}/c_{\text{n,max}}\ll 1/\mathcal{C}_{\text{e}} , which are true in practical situations.

Parameter combination	Required size	Interpretation
${𝒞}_{\text{e}}=I_{\text{typ}}^{\ast }{L}^{\ast }/\left(D_{\text{e,typ}}^{\ast }{F}^{\ast }{c}_{\text{n,max}}^{\ast }\right)$	$\ll 1$	Lithium-ion migration timescale is small relative to typical discharge timescale
$R^{\ast }{T}^{\ast }{\sigma }_{\text{k}}^{\ast }/\left(F^{\ast }{I}_{\text{typ}}^{\ast }{L}^{\ast }\right)$	$\gg 1$	Thermal voltage is large relative to the typical potential drop in electrode k
$R^{\ast }{T}^{\ast }{\kappa }_{\text{e,typ}}^{\ast }/\left(F^{\ast }{I}_{\text{typ}}^{\ast }{L}^{\ast }\right)$	$\gg 1$	Thermal voltage is large relative to the typical potential drop in the electrolyte
${\left(R_{\text{k}}^{\ast }\right)}^2{I}_{\text{typ}}^{\ast }/\left(D_{\text{s,k}}^{\ast }{F}^{\ast }{c}_{\text{n,max}}^{\ast }{L}^{\ast }\right)$	$\ll 1/{𝒞}_{\text{e}}$	Solid diffusion occurs on a shorter or similar timescale to a discharge
$I_{\text{typ}}^{\ast }/\left(m_{\text{k}}^{\ast }{a}_{\text{k}}^{\ast }{\left(c_{\text{e,typ}}^{\ast }\right)}^{1/2}{c}_{\text{n,max}}^{\ast }{L}^{\ast }\right)$	$\ll 1/{𝒞}_{\text{e}}$	Reactions occur on a shorter or similar timescale to a discharge

Equations

$${\Omega }_{\text{n}}^{\ast }=[0,L_{\text{n}}^{\ast }],{\Omega }_{\text{s}}^{\ast }=[L_{\text{n}}^{\ast },L^{\ast }-L_{\text{p}}^{\ast }],{\Omega }_{\text{p}}^{\ast }=[L^{\ast }-L_{\text{p}}^{\ast },L^{\ast }],$$

$${\varphi }_{\text{s,n}}^{\ast },{\varphi }_{\text{e,n}}^{\ast },c_{\text{e,n}}^{\ast },i_{\text{e,n}}^{\ast },N_{\text{e,n}}^{\ast }:$$

$${\varphi }_{\text{e,s}}^{\ast },c_{\text{e,s}}^{\ast },i_{\text{e,s}}^{\ast },N_{\text{e,s}}^{\ast }:$$

$${\varphi }_{\text{s,p}}^{\ast },{\varphi }_{\text{e,p}}^{\ast },c_{\text{e,p}}^{\ast },i_{\text{e,p}}^{\ast },N_{\text{e,p}}^{\ast }:$$

$$c_{\text{s,n}}^{\ast }:$$

$$c_{\text{s,p}}^{\ast }:$$

$$\text{partialderivative}{i}_{\text{e,k}}^{\ast }{x}^{\ast }=\{\begin{array}{ll}{a}_{\text{k}}^{\ast }{j}_{\text{k}}^{\ast },& \text{k}=\text{n, p},\\ 0,& \text{k}=\text{s},\end{array}$$

$$i_{\text{e,k}}^{\ast }={ϵ}_{\text{k}}^b{\kappa }_{\text{e}}^{\ast }(c_{\text{e,k}}^{\ast })\left(-\text{partialderivative}{\varphi }_{\text{e,k}}^{\ast }{x}^{\ast }+2(1-t^{+})\frac{R^{\ast }{T}^{\ast }}{F^{\ast }}\text{partialderivative}{x}^{\ast }\left(\log (c_{\text{e,k}}^{\ast })\right)\right),$$

$$I^{\ast }-i_{\text{e,k}}^{\ast }=-{\sigma }_{\text{k}}^{\ast }\text{partialderivative}{\varphi }_{\text{s,k}}^{\ast }{x}^{\ast },$$

$${ϵ}_{\text{k}}\text{partialderivative}{c}_{\text{e,k}}^{\ast }{t}^{\ast }=\text{partialderivative}{N}_{\text{e,k}}^{\ast }{x}^{\ast }+\frac{1}{F^{\ast }}\text{partialderivative}{i}_{\text{e,k}}^{\ast }{x}^{\ast },$$

$$N_{\text{e,k}}^{\ast }={ϵ}_{\text{k}}^b{D}_e^{\ast }(c_{\text{e,k}}^{\ast })\text{partialderivative}{c}_{\text{e,k}}^{\ast }{x}^{\ast }+\frac{t^{+}}{F^{\ast }}{i}_{\text{e,k}}^{\ast },$$

$$\text{partialderivative}{c}_{\text{s,k}}^{\ast }{t}^{\ast }=\frac{1}{(r^{\ast }{)}^2}\text{partialderivative}{r}^{\ast }\left((r^{\ast }{)}^2{D}_{\text{s,k}}^{\ast }\text{partialderivative}{c}_{\text{s,k}}^{\ast }{r}^{\ast }\right),$$

$$j_{\text{k}}^{\ast }=j_{\text{0,k}}^{\ast }\sinh \left(\frac{F^{\ast }{\eta }_{\text{k}}^{\ast }}{2R^{\ast }{T}^{\ast }}\right),$$

$$j_{\text{0,k}}^{\ast }=m_{\text{k}}^{\ast }(c_{\text{s,k}}^{\ast }{)}^{1/2}(c_{\text{s,k,max}}^{\ast }-c_{\text{s,k}}^{\ast }{)}^{1/2}(c_{\text{e,k}}^{\ast }{)}^{1/2}$$

$$\displaystyle\eta^{*}_{\text{k}}=\phi_{\text{s,k}}-\phi_{\text{e,k}}^{*}-U_{\text{k}}^{*}(c_{\text{s,k}}^{*}\big{|}_{r^{*}=R_{\text{k}}^{*}}),$$

$$\displaystyle i^{*}_{\text{e,n}}\big{|}_{x^{*}=0}=i^{*}_{\text{e,p}}\big{|}_{x^{*}=L^{*}}=0,$$

$$\displaystyle\phi^{*}_{\text{e,n}}\big{|}_{x^{*}=L^{*}_{\text{n}}}=\phi_{\text{e,s}}^{*}\big{|}_{x^{*}=L^{*}_{\text{n}}},\quad i^{*}_{\text{e,n}}\big{|}_{x^{*}=L^{*}_{\text{n}}}=i^{*}_{\text{e,s}}\big{|}_{x^{*}=L^{*}_{\text{n}}}=I^{*},$$

$$\displaystyle\phi^{*}_{\text{e,s}}\big{|}_{x^{*}=L^{*}-L^{*}_{\text{p}}}=\phi^{*}_{\text{e,p}}\big{|}_{x^{*}=L^{*}-L^{*}_{\text{p}}},\quad i^{*}_{\text{e,s}}\big{|}_{x=L^{*}-L^{*}_{\text{p}}}=i^{*}_{\text{e,p}}\big{|}_{x^{*}=L^{*}-L^{*}_{\text{p}}}=I^{*}.$$

$$\displaystyle N_{\text{e,n}}^{*}\big{|}_{x^{*}=0}=0,\quad N_{\text{e,p}}^{*}\big{|}_{x^{*}=L^{*}}=0,$$

$$\displaystyle c^{*}_{\text{e,n}}|_{x^{*}=L^{*}_{\text{n}}}=c^{*}_{\text{e,s}}|_{x^{*}=L^{*}_{\text{n}}},\quad N_{\text{e,n}}^{*}\big{|}_{x^{*}=L^{*}_{\text{n}}}=N_{\text{e,s}}^{*}\big{|}_{x^{*}=L^{*}_{\text{n}}},$$

$$\displaystyle c^{*}_{\text{e,s}}|_{x^{*}=L^{*}-L^{*}_{\text{p}}}=c^{*}_{\text{e,p}}|_{x^{*}=L^{*}-L^{*}_{\text{p}}},\quad N_{\text{e,s}}^{*}\big{|}_{x^{*}=L^{*}-L^{*}_{\text{p}}}=N_{\text{e,p}}^{*}\big{|}_{x^{*}=L^{*}-L^{*}_{\text{p}}}.$$

$$\displaystyle\partialderivative{c_{\text{s,k}}^{*}}{r^{*}}\bigg{|}_{r^{*}=0}=0,\quad-D^{*}_{\text{s,k}}\partialderivative{c_{\text{s,k}}^{*}}{r^{*}}\bigg{|}_{r^{*}=R_{\text{k}}^{*}}=\frac{j_{\text{k}}^{*}}{F^{*}},\quad$$

$$\displaystyle\phi^{*}_{\text{s,n}}\big{|}_{x^{*}=0}=0.$$

$$c_{\text{s,k}}^{\ast }(x^{\ast },r^{\ast },0)=c_{\text{s,k,0}}^{\ast },$$

$$c_{\text{e,k}}^{\ast }(x^{\ast },0)=c_{\text{e,typ}}^{\ast },$$

$$x^{\ast }=L^{\ast }x,t^{\ast }={\tau }_{\text{d}}^{\ast }t,D_{\text{e}}^{\ast }=D_{\text{e,typ}}^{\ast }{D}_{\text{e}},$$

$${\kappa }_{\text{e}}^{\ast }={\kappa }_{\text{e,typ}}^{\ast }{\kappa }_{\text{e}}{I}^{\ast }=I_{\text{typ}}^{\ast }I$$

$${\varphi }_{\text{s,n}}^{\ast }=\frac{R^{\ast }{T}^{\ast }}{F^{\ast }}{\varphi }_{\text{s,n}},{\varphi }_{\text{s,p}}^{\ast }=\left(U_{\text{p,ref}}^{\ast }-U_{\text{n,ref}}^{\ast }\right)+\frac{R^{\ast }{T}^{\ast }}{F^{\ast }}{\varphi }_{\text{p,k}},$$

$$\displaystyle\text{For }\text{k}\in\{\text{n, p}\}:\$$

$$j_{\text{k}}^{\ast }=\frac{I_{\text{typ}}^{\ast }}{a_{\text{k}}^{\ast }{L}^{\ast }}{j}_{\text{k}},j_{\text{0,k}}^{\ast }=\frac{I_{\text{typ}}^{\ast }}{a_{\text{k}}^{\ast }{L}^{\ast }}{j}_{\text{0,k}},m_{\text{k}}^{\ast }=m_{\text{k,typ}}^{\ast }{m}_{\text{k}},$$

$${\eta }_{\text{k}}^{\ast }=\frac{R^{\ast }{T}^{\ast }}{F^{\ast }}{\eta }_{\text{k}},U_{\text{k}}^{\ast }=U_{\text{k,ref}}^{\ast }+\frac{R^{\ast }{T}^{\ast }}{F^{\ast }}{U}_{\text{k}}.$$

$$\displaystyle\text{For }\text{k}\in\{\text{n, p}\}:\$$

$$\displaystyle\text{For }\text{k}\in\{\text{n, s, p}\}:\$$

$${\varphi }_{\text{e,k}}^{\ast }=-U_{\text{n,ref}}^{\ast }+\frac{R^{\ast }{T}^{\ast }}{F^{\ast }}{\varphi }_{\text{e,k}},i_{\text{e,k}}^{\ast }=I_{\text{typ}}^{\ast }{i}_{\text{e,k}}.$$

$$\text{partialderivative}{i}_{\text{e,k}}x$$

$$i_{\text{e,k}}$$

$$I-i_{\text{e,k}}$$

$${\mathcal{C}}_{\text{e}}{ϵ}_{\text{k}}{\gamma }_{\text{e}}\text{partialderivative}{c}_{\text{e,k}}t$$

$$N_{\text{e,k}}$$

$${\mathcal{C}}_{\text{k}}\text{partialderivative}{c}_{\text{s,k}}t$$

$$j_{\text{k}}$$

$$j_{\text{0,k}}$$

$${\eta }_{\text{k}}$$

$$\displaystyle i_{\text{e,n}}\big{|}_{x=0}=0,\quad i_{\text{e,p}}\big{|}_{x=1}=0,$$

$$\displaystyle\phi_{\text{e,n}}\big{|}_{x=L_{\text{n}}}=\phi_{\text{e,s}}\big{|}_{x=L_{\text{n}}},\quad i_{\text{e,n}}\big{|}_{x=L_{\text{n}}}=i_{\text{e,s}}\big{|}_{x=L_{\text{n}}}=I,$$

$$\displaystyle\phi_{\text{e,s}}\big{|}_{x=1-L_{\text{p}}}=\phi_{\text{e,p}}\big{|}_{x=1-L_{\text{p}}},\quad i_{\text{e,s}}\big{|}_{x=1-L_{\text{p}}}=i_{\text{e,p}}\big{|}_{x=1-L_{\text{p}}}=I.$$

$$\displaystyle N_{\text{e,n}}\big{|}_{x=0}=0,\quad N_{\text{e,p}}\big{|}_{x=1}=0,$$

$$\displaystyle c_{\text{e,n}}\big{|}_{x=L_{\text{n}}}=c_{\text{e,s}}|_{x=L_{\text{n}}},\quad N_{\text{e,n}}\big{|}_{x=L_{\text{n}}}=N_{\text{e,s}}\big{|}_{x=L_{\text{n}}},$$

$$\displaystyle c_{\text{e,s}}|_{x=1-L_{\text{p}}}=c_{\text{e,p}}|_{x=1-L_{\text{p}}},\quad N_{\text{e,s}}\big{|}_{x=1-L_{\text{p}}}=N_{\text{e,p}}\big{|}_{x=1-L_{\text{p}}}.$$

$$\displaystyle\partialderivative{c_{\text{s,k}}}{r_{\text{k}}}\bigg{|}_{r_{\text{k}}=0}=0,\quad-\frac{a_{\text{k}}\gamma_{\text{k}}}{\mathcal{C}_{\text{k}}}\partialderivative{c_{\text{s,k}}}{r_{\text{s,k}}}\bigg{|}_{r_{\text{k}}=1}=j_{\text{k}},\quad\text{k}\in\{\text{n, p}\}.$$

$$\displaystyle\phi_{\text{s,n}}\big{|}_{x=0}=0.$$

$$c_{\text{s,k}}(x,r,0)=c_{\text{s,k,0}},$$

$$c_{\text{e,k}}(x,0)=1,$$

$$V=\phi_{\text{s,p}}\big{|}_{x=1}-\phi_{\text{s,n}}\big{|}_{x=0}.$$

Full text

1 Introduction

1.1 Background

Lithium-ion batteries are used extensively in consumer electronics and industry, and as a result, understanding the physical processes occurring in such systems is key. Mathematical models are an essential tool for the design and management of battery systems. The standard mathematical model of a lithium-ion battery is the Doyle-Fuller-Newman (DFN) model, which was developed by John Newman and collaborators [1, 2, 3]. This model is sometimes referred to as the pseudo-two-dimensional (P2D) model or simply the Newman model. The model consists of a set of highly-coupled nonlinear parabolic and elliptic partial differential equations (PDEs). In the literature, this system of equations has been solved using a variety of different numerical methods including finite-difference methods, control volumes, finite-element methods, and orthogonal collocation, among others [1, 4, 5, 6, 7, 8, 9]. However, even when employing sophisticated numerical techniques, the DFN model remains too computationally complex for some applications. For instance, battery management systems (BMS) adopt equivalent circuit models, which consist of only a handful of ordinary differential equations (ODEs) and are favoured over the DFN model due to a combination of factors such as speed, memory requirements, and numerical convergence. Additionally, in the study of coupled systems of battery cells (e.g. thermal effects in a battery pack), spending computational resources (both RAM and CPU) on a model as detailed as the DFN model is often unnecessary. For these applications and others, physics-based models that are simpler than the DFN model are desired. The simplest of these, the SPM, has been employed in several settings in recent years [10, 11, 12, 13, 31]. There has also been a number of papers that provide justification for the SPM and suggest correction terms that may increase the accuracy of the predicted voltage [14, 15, 16, 17, 18, 19]. However, these approaches generally rely on a number of ad-hoc assumptions. In this work, we provide a systematic mathematical derivation of the SPM and an additional correction for the electrolyte by applying asymptotic methods to the DFN model. Asymptotic methods are widely applied within many subdisciplines of mathematics, but they are still relatively underutilized in battery modelling [20, 21]. However, they have been successfully applied in the derivation of the DFN model (by means of asymptotic homogenization) and in the derivation of reduced-order lead-acid battery models [22, 23]. Asymptotic methods have also been applied to the reduction of lithium-ion battery models that neglect the effects of the particles [24].

In our approach we consider approximations that can be found by exploiting two physically-relevant limits: i) that the electrical conductivity is large in the electrodes and electrolyte (such that the typical potential drop in each material is small relative to the thermal voltage), and ii) that the timescale of lithium-ion migration in the electrolyte is small relative to the typical timescale for a discharge. The validity of applying both of these limits is determined directly from the input parameter values which allows for the errors in the reduced models to be estimated a-priori. By comparison, in [10] for example, six assumptions that can only be validated a-posteriori by comparison with the full DFN model are required (e.g. that the current profile assumes a specific form). A key result of this work is the derivation of a single additional partial differential equation (PDE) and algebraic correction to the terminal voltage that accounts for nonuniform effects in the electrolyte and greatly improves the accuracy of the predictions when compared to the SPM. Additionally, we identify a key step overlooked in the ad-hoc derivations and show, through direct numerical comparison, that performing this step allows our reduced model to outperform the other ad-hoc models in the literature whilst retaining the same level of computational complexity.

After completion of the work in this paper, we became aware of another in-progress paper that employs asymptotic methods to simplify the DFN model [25]. In [25], the asymptotic limit of large changes in the open-circuit voltage (OCV) relative to the thermal voltage is taken; this is a different (and in fact complementary) limit to that taken here. Their limit recovers a variant of the SPM at leading order because the reaction overpotentials are small. In contrast, in our limit the small gradients in the electrolyte concentration, electrolyte potential, and electrode potentials give rise to the homogeneous behavior of electrode particles and hence the SPM. Whilst both limits recover a variant of SPM at leading order, the correction terms are different, for example [25] requires a second-order ODE for the electrolyte potential to be solved numerically. This results in a model that is slightly more computationally expensive than the model presented here but is significantly less computationally expensive than the DFN model. Therefore, if the model parameters for a particular case are not appropriate for the limit considered here, we encourage the reader to consider the reduced model in [25].

We begin by providing a brief overview of the DFN model, after which we nondimensionalize the model by introducing typical scalings. At this stage, we re-write the terminal voltage in electrode-averaged form, which is essential in the derivation of our reduced models. We then identify key dimensionless parameters and perform a systematic asymptotic reduction in the distinguished limit in which the electrical conductivities are of a comparable size to the ratio of the typical discharge timescale to the lithium-ion migration timescale. We make a uniform asymptotic expansion and, at leading order, recover the SPM. By extending the asymptotic expansion to first order, we obtain an additional PDE for the concentration of lithium ions in the electrolyte and an additional algebraic correction to the terminal voltage. We summarize our reduced model and compare its computational complexity and accuracy with the DFN model when the finite-volume discretization method is applied to both models. We also discuss the reduced computational complexity of our reduced model when alternative discretization methods are employed. Following our discussion on computational complexity, we compare our model with a selection of ad-hoc models from the literature, showing that our model best recovers the predictions of the DFN model. Finally, we provide a summary of the dimensional version of our model in Section 7, alongside a set of conditions which are to be satisfied for our model to be applicable and a set of a-priori error estimates.

2 Doyle-Fuller-Newman (DFN) model

Lithium-ion batteries consist of two electrodes, a porous separator, an electrolyte, and two current collectors, as displayed in Figure 1. Each electrode consists of active material particles within which lithium can be stored, and a binder (not shown) which holds the electrode together and maintains an electrical connection between the active material particles and the current collectors.

Upon discharge, lithium intercalated in the negative electrode particles diffuses to the surface of the particles where an electrochemical reaction occurs. This electrochemical reaction produces a lithium ion free to move through the electrolyte and an electron free to move through the electrode. The electron travels through the electrode, into the current collector, through an external wire, and towards the positive electrode. Meanwhile, the lithium ion migrates through the electrolyte towards the positive electrode. At the surface of the positive electrode particles, the lithium ion and the electron combine through another electrochemical reaction to form a lithium atom intercalated in the positive electrode particle. To charge the battery, a voltage is applied across the cell and the whole process occurs in reverse.

Here, we summarize the DFN model, which is the standard model of a lithium-ion battery [1, 2, 3]. The model is derived either by volume averaging [3] or the method of multiple scales [22]. Throughout, we use a superscript ‘$*$’ to indicate dimensional quantities. As indicated in Figure 1, the thicknesses of the negative electrode, separator, and positive electrode are $L^{*}_{\text{n}}$, $L_{\text{s}}^{*}$, and $L_{\text{p}}^{*}$, respectively. We denote the distance between the negative and positive current collectors by $L^{*}=L_{\text{n}}^{*}+L_{\text{s}}^{*}+L_{\text{p}}^{*}$. The active material particles in the negative and positive electrodes are assumed to be spheres with radii $R_{\text{n}}^{*}$ and $R_{\text{p}}^{*}$, respectively. Additionally, we assume that the behavior within each particle is spherically symmetric. We use the spatial coordinate $x^{*}\in[0,L^{*}]$ to indicate the location through the thickness of the battery and the spatial coordinate $r^{*}\in[0,R_{\text{k}}^{*}]$, $\text{k}\in\{\text{n, p}\}$ to indicate the location within each particle of active material. We define the following regions of the battery,

[TABLE]

which correspond to the negative electrode, separator, and positive electrode regions, respectively. We denote electric potentials by $\phi^{*}$, current densities by $i^{*}$, lithium concentrations by $c^{*}$ (in the electrolyte $c^{*}$ refers to lithium-ion concentrations), and molar fluxes by $N^{*}$. To indicate the region within which each variable is defined, we include a subscript $\text{k}\in\{\text{n, s, p}\}$ which corresponds to $\Omega^{*}_{\text{n}}$, $\Omega^{*}_{\text{s}}$, and $\Omega^{*}_{\text{p}}$, respectively. To distinguish variables in the electrolyte from those in the solid phase of the electrode, we employ an additional subscript ‘e’ for electrolyte variables and an additional subscript ‘s’ for solid-phase variables. For clarity, the variables in the model and their corresponding regions are

[TABLE]

We note that $c_{\text{s,n}}^{*}$ and $c_{\text{s,p}}^{*}$ depend on the macroscopic spatial variable, $x^{*}$, the microscopic spatial variable, $r^{*}$, and time, $t^{*}$, whereas all other variables depend on $x^{*}$ and $t^{*}$ only. When stating the governing equations, we take the region in which an equation holds to be implicitly defined by the subscript, $\text{k}\in\{\text{n, s, p}\}$, of the variables. With this in mind, the DFN model is summarized as:

[TABLE]

The functional forms for $U_{\text{n}}^{*}(c_{\text{s,n}}^{*})$, $U_{\text{p}}^{*}(c_{\text{s,p}}^{*})$, and $D_{\text{e}}^{*}(c_{\text{e}}^{*})$, which are the open-circuit potentials (OCPs) in the negative and positive electrodes and the lithium-ion diffusivity in the electrolyte, respectively, are taken from Newman’s DUALFOIL code [26]. A full list of parameters and their values is provided in Table 1. The values in this table are taken from Scott Moura’s fastDFN code [27], which are in turn adapted from the parameter values used in Newman’s DUALFOIL code [26]. These functions and parameters correspond to a cell with a graphite negative electrode, a \ceLiPF6 in EC:DMC electrolyte, and a lithium-cobalt-oxide positive electrode. In (1d), (1k), and (1l) the actual current density $I^{*}$ is distinct from the typical current density $I^{*}_{\text{typ}}$ in Table 1. We make this distinction to easily accommodate non-constant currents within the dimensionless form of our reduced models.

3 Asymptotic reduction of DFN model

3.1 Dimensionless form of DFN model

In the following, we will use asymptotic methods to systematically reduce the DFN model to simpler forms. In order to do this, we first nondimensionalize the DFN model by making the following scalings:

[TABLE]

We then identify the key timescales in the model, which are presented in Table 2. We also identify a set of dimensionless parameters, which we define in terms of the dimensional variables together with their physical meaning and calculated values in Table 3.

Substituting the scalings (2) into the dimensional model (1), we obtain the dimensionless version of the DFN model which is summarized as:

[TABLE]

[TABLE]

Before proceeding with model reduction, we note some helpful relations. The terminal voltage is given by

[TABLE]

We can re-write this in a more useful form by considering a particular path that the current follows through the cell. We consider the current entering through the negative current collector at $x=0$ and travelling through the solid phase of the negative electrode to a point $x_{\text{n}}\in[0,L_{\text{n}}]$. At this point, an electrochemical reaction occurs so that the current is transferred into the electrolyte. The current then travels through the electrolyte until it reaches a point $x_{\text{p}}\in[1-L_{\text{p}},1]$ where another electrochemical reaction occurs transferring the current into the solid phase of the positive electrode. Finally, the current travels through the positive electrode until it reaches the positive current collector. The terminal voltage can be written in terms of the potential drops associated with each section of this path as

[TABLE]

We define the pointwise OCV to be

[TABLE]

and the pointwise solid-phase Ohmic losses to be

[TABLE]

Then, using the definition of the reaction overpotential given in (3j), equation (6) becomes

[TABLE]

Equation (9) provides the voltage in terms the internal pointwise OCV, Ohmic losses, and overpotentials for a particular current path through the cell. However, we can also express the voltage in terms of the average (across all paths) OCV, Ohmic losses, and overpotentials by averaging (9) over each electrode. This is done by integrating (9) with respect to $x_{\text{n}}\in[0,L_{\text{n}}]$ and dividing by the negative electrode thickness, $L_{\text{n}}$, and then integrating with respect to $x_{\text{p}}\in[1-L_{\text{p}},1]$ and dividing by the positive electrode thickness, $L_{\text{p}}$. The result is that we can express the terminal voltage as

[TABLE]

where we use an overbar to represent the operation

[TABLE]

It is this electrode-averaged form of the voltage expression that will play a central role in our extension of the SPM. Without this electrode-averaged form, the extended model would need to be more complex than the model we present in order to obtain the same degree of accuracy.

Another useful property can be found by integrating (3e) with respect to $x$ over $\Omega_{\text{k}}$ for each $\text{k}\in\{\text{n, s, p}\}$, applying the boundary conditions (4b), (4c), (4d), (4e), (4f), (4g), and applying the initial condition (4k), to obtain

[TABLE]

Equation (12) is a statement of the conservation of the total number of lithium ions in the electrolyte.

3.2 The limit $\mathcal{C}_{\text{e}}\rightarrow 0$

We consider the limit of high electrical conductivity in the electrodes and electrolyte (such that the typical potential drop in each material is small relative to the thermal voltage) and the timescale for the migration of lithium ions in the electrolyte is small relative to the typical timescale of a discharge. This corresponds to the limit $\mathcal{C}_{\text{e}}\rightarrow 0$, where $\mathcal{C}_{\text{e}}$ is the ratio of the typical timescale for lithium-ion migration to the typical discharge timescale, $\sigma_{\text{k}}\rightarrow\infty$, where $\sigma_{\text{k}}$ is ratio of the thermal voltage to the typical Ohmic drop in the solid, and $\hat{\kappa}_{\text{e}}\rightarrow\infty$, where $\hat{\kappa}_{\text{e}}$ is the ratio of the thermal voltage to the typical Ohmic drop in the electrolyte. We take the distinguished limit in which $\sigma_{\text{k}}\mathcal{C}_{\text{e}}$ and $\hat{\kappa}_{\text{e}}\mathcal{C}_{\text{e}}$ both tend to a constant as $\mathcal{C}_{\text{e}}\rightarrow 0$, $\sigma_{\text{k}}\rightarrow\infty$, and $\hat{\kappa}_{\text{e}}\rightarrow\infty$ by setting

[TABLE]

We then expand all variables in powers of $\mathcal{C}_{\text{e}}$ in the form

[TABLE]

etc.

3.2.1 Leading-order model

At leading order in $\mathcal{C}_{\text{e}}$, (3e) and (3f) are

[TABLE]

Upon application of the leading-order boundary conditions (4e), (4f), (4g), the leading-order initial condition (4k), and the leading-order component of the condition (12), equation (13) yields

[TABLE]

Thus there is no depletion of the electrolyte at leading order. Equations (3c) and (3d), at leading order in $\mathcal{C}_{\text{e}}$, are then

[TABLE]

Since $c_{\text{e,k}}^{0}$, $\phi_{\text{e,k}}^{0}$, and $\phi_{\text{s,k}}^{0}$ are all independent of $x$ and $c_{\text{s,k}}^{0}$ is initially independent of $x$, it follows from the leading-order components of (3h), (3i), and (3j), that $j_{\text{k}}^{0}$, j_{\text{0,k}}^{0}, and $\eta_{\text{k}}^{0}$ do not depend on $x$. Therefore, we integrate (3b) with respect to $x$ over $\Omega_{\text{k}}$ for $\text{k}=\text{n, p}$ and apply (4b), (4c), (4d), to obtain

[TABLE]

From the leading-order component of (3h), we then have

[TABLE]

Employing (15) in conjunction with the leading-order interior boundary conditions (4c) and (4d), we obtain

[TABLE]

so that to leading order there is no potential drop in the electrolyte (for both the pointwise and electrode-averaged cases). Further, (15) implies that there are no solid-phase Ohmic losses (for both the pointwise and electrode-averaged cases) at leading order.

We are now in position to summarize the leading-order model. The leading-order description of lithium in the electrode particles is given by taking the leading-order component of (3g) and inserting (17) into the leading-order component of (4h) to obtain

[TABLE]

We identify (20) as the dimensionless form of the SPM [15, 28]. The name refers to the requirement to only solve a diffusion equation in one particle in each electrode, rather than solving a diffusion equation a particle at every macroscopic point as in the DFN model. This model should not be interpreted as replacing the many particles in an electrode by a single particle. Instead, in this limit all the particles in an electrode behave in exactly the same way and it is therefore sufficient to solve for just one representative particle.

3.2.2 First-order correction

We now proceed to calculate the $\mathcal{O}(\mathcal{C}_{\text{e}})$ corrections terms, which will extend the range of applicability of the reduced-order model to higher C-rates. At $\mathcal{O}(\mathcal{C}_{\text{e}})$, (3e) and (3f) become

[TABLE]

Integrating (21) with respect to $x$ over $\Omega_{\text{k}}$ for each $\text{k}\in\{\text{n, s, p}\}$ and applying the $\mathcal{O}(\mathcal{C}_{\text{e}})$ components of the boundary conditions (4e), (4f), and (4g), we obtain

[TABLE]

We substitute (23) into (22) and then integrate with respect to $x$ over $\Omega_{\text{k}}$ for each $\text{k}\in\{\text{n, s, p}\}$ using the $\mathcal{O}(\mathcal{C}_{\text{e}})$ components of the continuity boundary conditions (4f), (4g), and lithium-ion conservation condition (12) to determine the constants of integration. From this, we get

[TABLE]

At $\mathcal{O}(\mathcal{C}_{\text{e}})$, (3c) becomes

[TABLE]

We substitute (16) into (25), integrate with respect to $x$ over $\Omega_{\text{k}}$ for each $\text{k}\in\{\text{n, s, p}\}$, and determine two of the three constants of integration by applying the $\mathcal{O}(\mathcal{C}_{\text{e}})$ components of the interior boundary conditions (4c) and (4g) to obtain

[TABLE]

where $\tilde{\phi}_{\text{e}}$ is a constant, the form of which is provided in Appendix A.

At $\mathcal{O}(\mathcal{C}_{\text{e}})$, (3d) is

[TABLE]

which upon using (16) and integrating with respect to $x$, gives

[TABLE]

At $\mathcal{O}(\mathcal{C}_{\text{e}})$, (3b), (4b), (4c), and (4d) give

[TABLE]

Here we approach a key step in our derivation. Integrating (29a) with respect to $x$ over $\Omega_{\text{k}}$ for $\text{k}\in\{\text{n, p}\}$ and applying (29b) gives the conditions

[TABLE]

Thus the electrode-averaged corrections to the reaction currents $j_{\text{n}}$ and $j_{\text{p}}$ are zero. This means that after averaging the $\mathcal{O}(\mathcal{C}_{\text{e}})$ components of (3g) and (4h) in $x$ over $\Omega_{\text{k}}$ for each $\text{k}\in\{\text{n, p}\}$ and using (30), we get

[TABLE]

Crucially, there is no average flux on the surface $r_{\text{k}}=1$ in (31). The solution to (31) is therefore simply

[TABLE]

Therefore, the $\mathcal{O}(\mathcal{C}_{\text{e}})$ component of the electrode-averaged OCV is

[TABLE]

The $\mathcal{O}(\mathcal{C}_{\text{e}})$ components of (3h), (3j), and (3i) are

[TABLE]

where $\eta_{\text{k}}^{0}$ and j_{\text{0,k}}^{0} are given by (18) and (20f), respectively. Averaging (34a) over electrode $\Omega_{\text{k}}$, gives

[TABLE]

We can also easily obtain $\overline{\phi}_{\text{e,k}}^{1}$ and $\overline{\phi}_{\text{s,k}}^{1}$ by electrode-averaging (26) and (28). The expressions for these can be found in Appendix B.

We now have all the components necessary to obtain the first-order correction to the voltage. Inserting (33), (35), and the electrode-averaged potential expressions in Appendix B into (10), we obtain

[TABLE]

This is a purely algebraic expression and is therefore obtained without the numerical solution of any additional PDEs. As a result, no additional states need to be stored at each time step, which is beneficial for microcontrollers where RAM is limited. Further, the additional computational cost compound to the SPM (20) is evaluating this algebraic expression and is therefore negligible. The first term of (36) corresponds to corrections to the reaction overpotentials due to concentration variations in the electrolyte, the second to the concentration overpotential, the third to electrolyte Ohmic losses, and the fourth to solid-phase Ohmic losses.

3.3 Combined voltage expression

To write an expression for the combined leading- and first-order voltage, $V=V^{0}+\mathcal{C}_{\text{e}}V^{1}$, we define the electrode-averaged exchange-current densities to be

[TABLE]

and use the fact that

[TABLE]

The combined voltage accurate to $\mathcal{O}\left(\mathcal{C}_{\text{e}}^{2}\right)$ is then

[TABLE]

with the combined leading- and first-order electrode-averaged components of the voltage expression given by

[TABLE]

We refer to this model that includes the first-order correction to the terminal voltage as the SPMe(S). That is, the single particle model with electrolyte with the (S) referring to the fact that this formulation considers the electrolyte to be in steady state.

Before proceeding, we would like to draw further attention to what we view as the key step in the derivation of the SPMe(S), namely the electrode-averaging step. Electrode averaging is essential as it provides us with a well-defined problem: given the electrode-averaged current, find the electrode-averaged potential differences. If instead we try to evaluate the pointwise voltage expression, we must determine both the current (since we only know the electrode-averaged current) and potential difference at a particular point. To get around this issue, it is typical for ad-hoc models in the literature, which use a pointwise voltage expression, to implicitly assume that the current at a particular point is equal to the electrode-averaged current. This is not in general true, since it implies that the concentrations in every particle of the DFN model would be the same across all operating conditions. Of course, this is not a bad first assumption to make in the limit we have been considering, and it is indeed valid for very low currents; this is the reason the SPM is a reasonable approximation. However, we have shown systematically in our asymptotic expansion that this assumption is only true at leading order and not at first order. Electrode averaging removes the requirement for this assumption, and ensures our expressions are also valid at first order, so we gain additional accuracy over the ad-hoc models for negligible additional computational cost.

4 Canonical SPMe

The SPMe(S) holds for cases where the electrolyte can be taken to be in quasi-steady state (e.g. when the current varies over a longer timescale than the electrolyte diffusion timescale). However, for many applications transient effects in the electrolyte are important, particularly after a step change in current. These transient effects occur on the timescale of migration of lithium ions in the electrolyte. To include these effects it is therefore natural to scale time with the timescale for migration of lithium ions, $t^{*}=\tau_{\text{e}}\tilde{t}$. On this short migration-timescale, at leading and first order the concentrations in the electrode particles remain constant and the terms describing exchange of lithium within the electrolyte equation are negligible. We do not present the corresponding systematic asymptotic analysis for this migration-timescale problem here but instead simply state the composite model which produces the correct result on both the diffusion and discharge timescales. The main difference between this model and the SPMe(S) is that we must now solve a PDE to obtain the first-order correction for the electrolyte concentration, $c_{\text{e,k}}^{1}$. We shall take this model to be the canonical SPMe and therefore simply refer to it as the SPMe. The SPMe is summarized as:

[TABLE]

Here the overbar terms are electrode-averaged quantities. The SPMe (40) consist of two independent linear PDE problems that describe the concentration of lithium in the negative and positive particles and an independent linear PDE problem that describes the concentration of lithium ions in the three regions of the electrolyte. The terminal voltage is obtained post-solution through a simple and easily interpreted algebraic expression. Since all three PDE problems are independent, the problem has a naturally parallel structure. The linearity of the PDEs is also advantageous for the application of numerical methods and the determination of simpler/analytic solutions in special cases (e.g. the SPMe(S)).

5 Model comparisons

5.1 Finite-volume implementation

In this section, we compare the DFN model (1), SPM (20), and the SPMe (40). We implement the DFN model by discretizing the spatial dimensions using the finite-volume method to convert the system of PDEs into a system of differential algebraic equations (DAEs) of index one. Before solving this system, a set of consistent initial conditions for the potentials are found numerically using Newton’s method. The time evolution of the system is then performed using the SUNDIALs DAE solvers interfaced though PyBaMM (Python Battery Mathematical Modelling) [29, 30]. PyBaMM is a battery modelling software implemented in Python which is designed to facilitate the comparison of battery models by providing a common interface to discretization methods and numerical solvers. Similarly, we use the finite-volume method to discretize the spatial dimensions of the SPM and SPMe and again use SUNDIALs for the time evolution. We use the same mesh to discretize the SPM, SPMe, and DFN model. In the $x$-direction, we use $30$ points in the negative electrode, $20$ points in the separator, and $30$ points in the positive electrode. In the $r$-direction, we use $15$ points. Numerical errors associated with the spatial discretization, are therefore of order $10^{-2}$. However, we have limited their influence upon the comparison by applying the same numerical method to each model. Since we aim to compare models and not numerical methods, we only concern ourselves with the relative reduction in computational complexity obtained by using the SPMe instead of the DFN model, whilst retaining the same numerical method. Here, we use the finite-volume method, but of course alternative methods could be applied to both models to further increase speed and/or reduce memory requirements of each.

Computational complexity consists of space complexity (memory) and time complexity. To demonstrate the reduction in space complexity, we consider our finite-volume implementation with $30$ points in the electrodes, $20$ points in the separator, and $15$ points in the particles. In this case, the DFN model requires $2\times(30\times 15)=900$ states for the concentration of lithium in the particles, $80$ states for the concentration of lithium ions in the electrolyte, $80$ states for the electrolyte potential, and $60$ states for the electrode potentials. This leads to a total of $1120$ internal states, which are to be stored at each time step. On the other hand, the SPMe requires only $30$ states for the concentration of lithium in the particles and $80$ states for the concentration of lithium ions in the electrolyte. Therefore, for this discretization the SPMe requires just over $10\%$ of the memory required by the DFN model. With regards to time complexity, the DFN model is limited in two respects. Firstly, it is limited by the large number of states which must be computed at every time step and secondly it is limited by the stiff nature of the system, which arises from the discretization of the mixed parabolic and elliptic PDEs resulting in a system of DAEs. As previously discussed, the SPMe addresses the first of these limitations, however, it also addresses the second limitation because the model consists of only three linear parabolic PDEs, which upon discretization lead to a well-conditioned system of ODEs amenable to taking larger time steps than the stiff DAE system of the DFN model. Furthermore, being a system of DAEs the DFN model suffers from convergence issues if inconsistent initial conditions are provided or if there is a large change in current, for example, in switching between charging and discharging. This is an inherent issue with DAE based algorithms and efforts have been made to limit the range of inconsistent initial conditions [7]. However, this robustness problem persists and convergence of the DFN model cannot always be ensured for non-constant currents.

5.2 Voltage comparison

We compare the SPM, SPMe, and DFN model by considering the case of a single constant-current discharge over a range of C-rates. The initial stoichiometries of the negative and positive electrodes are 0.8 and 0.6, respectively, and we terminate the discharge when the terminal voltage reaches $3.2\text{\,}\mathrm{V}$. For this cell a C-rate of $1\text{\,}\mathrm{C}$ corresponds to a discharge current density of $24\text{\,}\mathrm{A}\mathrm{/}\mathrm{m}^{2}$. As provided in Table 3, we have $\mathcal{C}_{\text{e}}=5.1\times 10^{-3}\mathcal{C}$ where $\mathcal{C}$ is the C-rate. The predicted terminal voltage of each model is presented in Figure 2.

At low C-rates, all three models match well with a RMS voltage error of just $1.72\text{\,}\mathrm{m}\mathrm{V}$ for the SPM at $0.1\text{\,}\mathrm{C}$, as expected. However, at higher C-rates we observe that the SPM prediction deviates from the DFN model solution, producing a RMS voltage error of $19.86\text{\,}\mathrm{m}\mathrm{V}$ at $1\text{\,}\mathrm{C}$ and $62.78\text{\,}\mathrm{m}\mathrm{V}$ at $3\text{\,}\mathrm{C}$. The SPMe greatly improves upon this, with a RMS voltage error of just $3.04\text{\,}\mathrm{m}\mathrm{V}$ at $1\text{\,}\mathrm{C}$ and $13.34\text{\,}\mathrm{m}\mathrm{V}$ at $3\text{\,}\mathrm{C}$. Unfortunately, there is a discrepancy in the voltage curves near then end of the discharge at $3\text{\,}\mathrm{C}$. To investigate the source of this discrepancy, we plot the error in each component of the voltage during a $3\text{\,}\mathrm{C}$ discharge in Figure 3.

We observe that for the majority of the discharge all components of the SPMe voltage agree well with the voltage predicted by the DFN model. However, near the end of the discharge, there is a large increase in the error of our solution, as observed in Figure 2. Around 60-70% of this error is due to a poor estimation of the electrode-averaged OCV. This error occurs when the OCV becomes highly nonlinear. If we extend our asymptotic expansion of the OCP, $U_{\text{k}}(c_{\text{k}})$ to second order, we obtain the term $\mathcal{C}_{\text{e}}^{2}U_{\text{k}}^{\prime\prime}(c_{\text{s,k}}^{0})(c_{\text{s,k}}^{1})^{2}/2$. Hence, when the OCV is highly nonlinear, $U_{\text{k}}^{\prime\prime}(c_{\text{s,k}}^{0})$ becomes large and the higher-order terms neglected by the SPMe become significant. To account for this behavior, we can consider the distinguished limit in which

[TABLE]

In this limit, we cannot avoid solving for the concentrations in all the particles in each electrode and therefore lose a large portion of the computational simplicity of the SPM and SPMe. We have developed and implemented a numerical scheme for this limit and found that it does indeed correct the voltage discrepancy. However, because of the increased computational complexity, we do not discuss it in detail within this paper.

5.3 Grid dependence and computation time

In Figure 4, we compare the solutions of the SPMe and DFN model with $5$, $10$, $20$ and $30$ points in each domain (negative electrode, separator, positive electrode, negative particle, positive particle) across a range of C-rates. Here, we measure the RMS voltage error with respect to the DFN model solution with $30$ points in each domain. It is important here to only focus upon the relative timings of the models and not absolute times. Of course, a pure C++ (instead of Python) implementation of the models with a more sophisticated numerical method would increase speed of both models, but that has not been our focus. It should also be noted that this comparison is for a constant-current discharge; the non-constant case has been reported to give rise to even longer computation times for the DFN model [4]. The key observation to make from Figure 4 is that for a relatively small increase in RMS voltage error, particularly at low C-rates, an order of magnitude decrease in computation time is achieved by using the SPMe instead of the DFN model. Achieving such large decreases in computation time is consistent across all C-rates and numbers of grid points used. Further, we can observe that the SPMe generally increases the accuracy of the SPM by an order of magnitude for a particular C-rate whilst maintaining a similar computation time. When a small number of grid points (e.g. 10 points) are used, and a current with a C-rate above $0.5\text{\,}\mathrm{C}$ is applied, the SPMe produces a solution that is not limited by the discretization error; we discuss the limiting asymptotic errors in Section 7. Therefore, to achieve the asymptotic accuracy of the SPMe, a coarser spatial discretization is often sufficient, which could also allow for further increases in speed and savings in memory.

5.4 Comparison of internal states

To further confirm the accuracy of the SPMe, we compare the internal states predicted by the SPMe and DFN model. These are presented for a $1\text{\,}\mathrm{C}$ constant current discharge in Figure 5. We observe good agreement between the two sets of model predictions. However, two key discrepancies can be observed: the first in the negative electrode stoichiometry and the second in the electrolyte concentration at late times. We note that the apparent discrepancy in the negative electrode potential is only due to the scale employed in Figure 5. This discrepancy is in fact $\mathcal{O}(\mathcal{C}_{\text{e}}^{2})$.

The discrepancy in the negative electrode stoichiometry is a result of the SPMe using only the leading-order equations in $\mathcal{C}_{\text{e}}$ within the electrode particles; this is equivalent to considering only the electrode-averaged concentration in the particles. This is the same approximation as employed by the SPM and the SPMe is therefore only accurate to $\mathcal{O}(\mathcal{C}_{\text{e}})$ for the concentration in the particles. It is possible to additionally solve for the first-order correction in the particles, but this requires solving a diffusion equation in each particle instead of a single diffusion equation in the electrode-averaged particle. This increases the computational complexity of the model and has therefore been omitted from this paper.

The discrepancy in the electrolyte concentration at late times is fundamentally connected to the discrepancy in the voltage curves observed in Figure 2. As already discussed, this is due to nonlinear nature of the OCV and the term $U_{\text{k}}^{\prime\prime}(c_{\text{s,k}}^{0})$ becoming large. We must then consider the distinguished limit in which (41) holds. In this limit, we consider a heterogeneous interfacial current density, and so the electrode-averaged current density source/sink terms in (40c) and (40d) are replaced by heterogeneous versions. We have implemented a numerical scheme for this limit and can confirm that this discrepancy is accounted for in this way. However, this limit requires one to solve for the concentration in each of the particles and is therefore much more computationally expensive than the models of concern in this paper.

In Figure 6, we compare the internal states of the SPMe and DFN model across a range of C-rates. Here, we display the lithium concentration in the x-averaged particle to demonstrate the ability of our model to capture the concentration profile inside particles as well as the surface concentration. As before, we observe good agreement of the internal states, with the exception of the electrolyte concentrations and the electrolyte potentials at large C-rates. Both of these discrepancies are a result of the the nonlinearity of the OCV, as discussed previously.

5.5 Alternative numerical methods

In the previous section, we only discussed the benefits of the SPMe compared to the DFN model when both models were implemented using the finite-volume method. However, there are other more sophisticated methods that have been applied to the DFN model, which can also be applied to the SPMe (e.g. [9]). The approach implemented in [9] involves volume averaging the PDEs in the particles and employing a polynomial approximation to obtain the surface concentration. As a result, in their formulation of the DFN model (Table IV of [9]), there are 4 PDEs, and 1 algebraic equation in each electrode along with 2 PDEs in the separator to be solved. Upon applying their orthogonal collocation discretization procedure, $5(N_{\text{n}}+1)+2(N_{\text{s}}+1)+5(N_{\text{p}}+1)$ DAEs are obtained (see just after equation (26) of [9] for more details), where $N_{\text{k}}$ is the number of cosine basis functions used in each domain. If we were to reformulate the governing equations in the particles in the SPMe using their polynomial approximation approach, we would have 1 PDE for the electrolyte concentration and 1 ODE for the $r$-averaged lithium concentration in the average particle in each electrode (note that the algebraic constraint for the surface concentration can now be solved analytically) and 1 PDE for the electrolyte concentration in the separator. Therefore, upon application of their discretization, we obtain $((N_{\text{n}}+1)+1)+(N_{\text{s}}+1)+((N_{\text{p}}+1)+1)$ ODEs. Using a single cosine in each domain results in $8$ ODEs for the SPMe compared to $24$ DAEs for the DFN model. Alternatively, using $N_{\text{n}}=7$, $N_{\text{s}}=3$, and $N_{\text{p}}=7$ produces a system of 88 DAEs for the DFN model but only $22$ ODEs for the SPMe. Therefore, applying this discretization to the SPMe instead of the DFN model can result in the memory requirements being reduced to below a third. This is a step towards employing a practical implementation of a physics-based model in a BMS where RAM is limited. Alternatively, a system of battery packs consisting of more than three times as many cells can be solved with similar RAM requirements using the SPMe instead of the DFN model. In addition, it is also typically faster to solve fewer ODEs (two of which are independent of the others) than a greater number of highly-coupled DAEs. Furthermore, the convergence issues related to DAE systems are not an issue for the ODE system produced by discretizing the SPMe. These types of benefits will be consistent across different discretization approaches.

To provide an example of the benefits of using the SPMe with the discretization in [9] instead of the DFN model, we consider the case where there is a memory limit such that only $24$ states can be stored at each time step. In this case, the DFN model is restricted to employing just one cosine basis function in each domain, whereas the SPMe can make use of $N_{\text{n}}=7$, $N_{\text{s}}=3$, $N_{\text{p}}=7$ basis functions. The RMS voltage error for a $1\text{\,}\mathrm{C}$ discharge for the DFN model using a single cosine function in each domain relative to a full finite difference implementation (with 50 points in each electrode and 35 points in the separator) is $17.84\text{\,}\mathrm{m}\mathrm{V}$, as reported in [9]. In [9], using $N_{\text{n}}=7$, $N_{\text{s}}=3$, $N_{\text{p}}=7$ basis functions, the RMS voltage error for the DFN model is just $0.91\text{\,}\mathrm{m}\mathrm{V}$. If we assume that this numerical error carries across to an SPMe implementation (this appears to be the case in Figure 4, for example for $0.1\text{\,}\mathrm{C}$ with $5\text{\,}\mathrm{p}\mathrm{t}\mathrm{s}$), then the main source of error will be modelling error. For a $1\text{\,}\mathrm{C}$ discharge, the SPMe produces an error of just $3.04\text{\,}\mathrm{m}\mathrm{V}$ in the finite-volume implementation, see Table 4 (this RMS voltage error also holds for 50 points in the electrodes and 35 points in the separator). Assuming that this error is entirely modelling error, then we expect the RMS voltage error of a $N_{\text{n}}=7$, $N_{\text{s}}=3$, $N_{\text{p}}=7$ discretization of the SPMe to be approximately $3.04\text{\,}\mathrm{m}\mathrm{V}$. Therefore, with the same memory budget employing the SPMe instead of the DFN model can improve the accuracy of the voltage prediction by an order of magnitude.

6 Critical assessment of variations on the SPMe in the literature

6.1 Overview of models

There are a number of alternative models in the literature that extend the SPM in an ad-hoc manner to account for electrolyte effects [14, 15, 16, 17, 18, 19]. In this section, we highlight the key differences between these models and the canonical SPMe (40) presented here. We have chosen to compare a subset of the models, which cover the variety of ad-hoc models available. To do this, we have converted the models to dimensionless form using the scalings in (2). In some papers the model is discretized during development. We view the choice of discretization to be a numerical method instead of a feature of the model itself. Therefore, we have converted each model into continuum form so as to only highlight the differences in the underlying models. We do not aim to study the benefits and drawbacks of different numerical methods.

A common theme in the models from the literature is to replace the electrode-averaged concentration overpotential and electrode-averaged electrolyte Ohmic losses with pointwise versions. It is also common to neglect the solid-phase Ohmic losses (this is a reasonable assumption since we have already observed these to be small). The general form of the terminal voltage expression used in the literature is then:

[TABLE]

This expression consists of a combination of both electrode-averaged and pointwise terms, and therefore an accuracy of $\mathcal{O}(\mathcal{C}_{\text{e}}^{2})$ cannot be ensured.

We begin by considering the model proposed by Perez et. al. in [15]. Firstly, the electrolyte flux (40d) is replaced by the expression (after converting to our notation)

[TABLE]

Due to the the presence of the nonlinear diffusion coefficient, the electrolyte problem is nonlinear, while that for the SPMe (40) is linear. Since the other terms are unchanged, this nonlinear form is also accurate up to $\mathcal{O}(\mathcal{C}_{\text{e}}^{2})$. Equation (43) includes some (but not all) higher-order terms, which in practice may increase accuracy but this cannot be ensured (there is also a chance it could reduce accuracy). Secondly, the electrode-averaged concentration overpotential (40n) and the electrode-averaged electrolyte Ohmic losses (40q) are replaced by their pointwise versions:

[TABLE]

respectively (note that to get (45) we have corrected the sign of the expression in [15]). Here, $\overline{\kappa}_{\text{e}}^{\text{eff}}$ is the dimensionless effective conductivity averaged across the entire cell, with the effective conductivity defined by $\kappa_{\text{e}}^{\text{eff}}(c_{\text{e,k}})=\epsilon_{\text{k}}^{b}\kappa_{\text{e}}(c_{\text{e,k}})$. The terminal voltage is then given by (42). Additionally, (45) requires that $\kappa_{e}^{\text{eff}}(c_{\text{e,k}})\approx\overline{\kappa}_{e}^{\text{eff}}$ throughout the cell. With this assumption, $\mathcal{O}(\mathcal{C}_{\text{e}}^{2})$ accuracy cannot be ensured for all values of $\epsilon_{\text{n}}$, $\epsilon_{\text{s}}$, and $\epsilon_{\text{p}}$. Finally, solid-phase Ohmic losses are neglected and Ohmic losses due to the presence of SEI are included; we shall neglect the SEI terms in our comparisons, noting that (40) can be easily extended to include them.

The model presented by Prada et. al. [16] also employs (44) for the concentration overpotential. However, the electrolyte Ohmic losses are taken to be

[TABLE]

where

[TABLE]

Whilst (46) does not rely upon the assumption that $\kappa_{\text{e}}^{\text{eff}}(c_{\text{e,k}})\approx\overline{\kappa}_{\text{e}}^{\text{eff}}$, its form is still a result of considering the pointwise electrolyte potential difference instead of the electrode-averaged difference. In addition to these differences, Prada et. al. [16] take the exchange-current densities $\overline{j}_{\text{0,n}}$ and $\overline{j}_{\text{0,p}}$ to be constant. In terms of accurately reproducing the results of the DFN model, this simplification has a clear disadvantage as the reaction overpotentials are strong functions of the lithium and lithium-ion concentrations.

The model developed by Han et. al [17] is the same as that presented by Prada et. al. [16] without the additional assumption of constant exchange-current densities. That is, Han et. al. [17] employ (44) and (46), which are the pointwise concentration overpotential, and electrolyte Ohmic losses, respectively. Han et. al. [17] note the tendency for their model to over-correct the voltage when accounting for electrolyte effects. We suspect the use of pointwise terms is the cause.

The model presented by Kemper et. al. [14] is somewhat different to the others we have discussed. Firstly, the model is presented as a set of ODEs instead of PDEs. These ODEs are derived by spatially discretizing the underlying PDEs. Whilst this particular discretization may be useful, we consider this to be a numerical method and not a feature of the model itself. Since we aim to compare the underlying simplified models directly, we have converted these ODEs back into PDEs. The resulting PDEs that describe the concentrations in the electrode particles and the electrolyte, are equivalent to those used in our model. However, the expression for the terminal voltage is very different and it is not easy to prescribe meaning to each of the individual components; we have however attempted this. We have converted the voltage from [14] into dimensionless form and provided details in Appendix C to be clear about the exact model we are comparing. It was unclear if the components of the voltage correspond to electrode-averaged or pointwise quantities so we have left this unstated.

6.2 Model comparison

We compare the variations of the SPMe in the literature and our canonical SPMe. For this section, we use a finite-volume implementation of each of the models implemented in MATLAB with ODE15s being used for the time integration. We consider a range of constant-current discharge rates and then consider the RMS voltage error of each model relative to the DFN model. For each model, we discretize using 30 points in each electrode, 20 in the separator, and 15 in each particle. Our results are presented in Figure 7, where we compare the models in [15, 17, 14]. Each version of the SPMe consists of three parabolic PDEs, one in the negative particle, one in the positive particle, and one in the electrolyte alongside an algebraic expression for the voltage. As a result, each SPMe must store $2\times 15+(30+20+30)=120$ states at each time step. Therefore, the memory requirements of each model are similar. Upon evaluation, each version of the model takes on average $0.07\text{\,}\mathrm{s}$ (at the time of writing the ODE15s implementation is faster the Python–SUNDIALs implementation). We observe that across all discharge rates, our canonical SPMe outperforms the other models from the literature. In particular, our canonical SPMe is consistently an order of magnitude more accurate than the models in Perez et. al. [15] and Kemper et. al [14]. Furthermore, at higher C-rates, the RMS voltage errors in the models from the literature approach being of the order of $0.1\text{\,}\mathrm{V}$ whereas the RMS voltage errors of our SPMe only reach the order of $0.01\text{\,}\mathrm{V}$. Additionally, as we would expect, our model converges to the DFN model solution at a faster rate than the other models. We attribute the main gains of our model to the consistent electrode-averaged OCPs, overpotentials, and Ohmic losses in our terminal voltage expression.

7 Dimensional model summary and conditions for application

We now present a summary of our dimensional SPMe alongside the conditions that should be met for the model to be valid. The purpose of this section is to serve as a reference, which is easily accessible to those not interested in the full details of the model derivation. For the purpose of this section, we will drop the superscript notation, which was used to indicated the asymptotic order of a variable. Also, to be consistent with the rest of the paper, we employ a superscript $*$ to denote dimensional quantities. By reapplying our scalings in (2) to (40) and combining the leading- and first-order equations for the electrolyte concentration, we obtain:

[TABLE]

We have provided the conditions that ensure the validity of (48) in Table 6.

If the conditions in Table 6 are met then the model error at a particular time is of size

[TABLE]

where the term \big{|}(U_{\text{k}}^{*})^{\prime\prime}\big{|} is the absolute value of the second derivative of the OCP in electrode k. This term is responsible for the decrease in accuracy when the OCP is significantly nonlinear, as discussed in Section 5. In situations where transient effects are not important (for example when the current varies on timescales longer than the electrolyte diffusion timescale), the same degree of accuracy can be achieved by employing the dimensional equivalent to the SPMe(S); this is achieved by neglecting the time derivative term in the electrolyte equation (48c).

8 Summary and further work

We have systematically derived simplified mathematical models from the standard DFN model through the use of asymptotic methods. The leading-order model, the SPM (20) is commonly used in the control community. By considering higher-order effects, we have extended this model to develop a canonical SPMe that is more accurate and applicable over a larger range of operating conditions. The canonical SPMe has been shown to give good agreement with the DFN model whilst providing dramatic decreases in computational complexity in the form of both reduced memory requirements and reduced computation time, both of which are highly desirable features for BMS, parameter estimation, and optimization. We have also shown our version of the SPMe be more accurate than the other reduced models of similar computational complexity available in the literature. A key result of this paper is to identify the requirement of writing the output voltage expression in terms of the electrode-averaged OCV, overpotentials, and Ohmic losses. This step has been overlooked in previous literature. Finally, our systematic approach has allowed us to identify the reasons for discrepancies in the predictions of the SPMe and identify the minimal extensions required for them to be corrected.

There are a number of possible additional physical mechanisms that are of interest to incorporate into a battery model. These include mechanical effects, thermal effects, and degradation mechanisms. One approach would be to simply introduce these effects in an ad-hoc manner to existing simple models without consideration of their interactions within the context of a more complicated model such as the DFN model. However, as demonstrated here, it is important to derive reduced-order models in a systematic fashion; not doing so can lead to inconsistent terms and unnecessary loss of accuracy. The approach in this paper can be applied to models that include additional physical effects, ensuring that the resulting reduced-order models retain the most important terms to faithfully represent the underlying physics. Further, by conducting asymptotic analysis, the error introduced in developing simplified models can be properly quantified in terms of parameter groupings, allowing for the applicability of a model to be determined a-priori.

Appendix A Electrolyte constants of integration

[TABLE]

[TABLE]

Appendix B Electrode-averaged quantities

[TABLE]

Appendix C Dimensionless voltage from [14]

Although the form of the voltage in [14] is not given explicitly, we assume that it is given by $\phi_{\text{s}}(0^{+},t)-\phi_{\text{s}}(0^{-},t)$ (in their notation), which is the potential difference between the solid-phase potential evaluated at the positive and negative current collectors. We also note that in [14] current has been defined in the direction of flow of positive charge (opposite of the convention used here) so we account for the sign change here. The model in [14] is then taken to be:

[TABLE]

[TABLE]

where $\overline{\kappa}_{\text{e,k}}$ is given by (47) and

[TABLE]

Appendix D Tables

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 Siemens Corporate Technology and BBOXX. This work was also supported with funding provided by The Faraday Institution, grant number EP/S003053/1, FIRG003.