# An asymptotic derivation of a single particle model with electrolyte

Scott G. Marquis<EMAIL_ADDRESS>222Mathematical Institute,
University of Oxford, OX2 6GG, United Kingdom Valentin Sulzer 22footnotemark:
2 Robert Timms 22footnotemark: 2 333The Faraday Institution Colin P. Please
22footnotemark: 2 33footnotemark: 3 S. Jon Chapman 22footnotemark: 2
33footnotemark: 3

###### Abstract

The standard 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. Recently,
there has been a move to extend the SPM to include electrolyte effects to
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 so that these models are not 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
electrolyte behavior. Firstly, this allows us to quantify the error in the
reduced model in terms of ratios of key parameters, from which the range of
applicable operating conditions can be determined. Secondly, by comparing our
model with ad-hoc models from the literature, we show that existing models
neglect 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 over existing models whilst
maintaining the same degree of computational complexity.

##  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.

Figure 1: Schematic of a lithium-ion battery. Active material particles are
shaded cross-stitch and diagonally for the negative and positive active
materials, respectively.

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,

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

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

$\displaystyle\phi_{\text{s,n}}^{*},\ \phi_{\text{e,n}}^{*},\
c_{\text{e,n}}^{*},\ i_{\text{e,n}}^{*},\ N_{\text{e,n}}^{*}:\quad$
$\displaystyle x^{*}\in\Omega_{\text{n}}^{*},$
$\displaystyle\phi_{\text{e,s}}^{*},\ c_{\text{e,s}}^{*},\
i_{\text{e,s}}^{*},\ N_{\text{e,s}}^{*}:\quad$ $\displaystyle
x^{*}\in\Omega_{\text{s}}^{*},$ $\displaystyle\phi_{\text{s,p}}^{*},\
\phi_{\text{e,p}}^{*},\ c_{\text{e,p}}^{*},\
i_{\text{e,p}}^{*},N_{\text{e,p}}^{*}:\quad$ $\displaystyle
x^{*}\in\Omega_{\text{p}}^{*},$ $\displaystyle c_{\text{s,n}}^{*}:$
$\displaystyle r^{*}\in[0,R_{\text{n}}^{*}],\quad
x^{*}\in\Omega_{\text{n}}^{*},$ $\displaystyle c_{\text{s,p}}^{*}:$
$\displaystyle r^{*}\in[0,R_{\text{p}}^{*}],\quad
x^{*}\in\Omega_{\text{p}}^{*}.$

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:

Governing equations _Charge conservation:_
$\displaystyle\partialderivative{i_{\text{e,k}}^{*}}{x^{*}}=\begin{cases}a_{\text{k}}^{*}j_{\text{k}}^{*},\quad&\text{k}=\text{n,
p},\\\ 0,\quad&\text{k}=\text{s},\end{cases}$
$\displaystyle\text{k}\in\\{\text{n, s, p}\\},$ (1b) $\displaystyle
i_{\text{e,k}}^{*}=\epsilon_{\text{k}}^{b}\kappa_{\text{e}}^{*}(c_{\text{e,k}}^{*})\left(-\partialderivative{\phi_{\text{e,k}}^{*}}{x^{*}}+2(1-t^{+})\frac{R^{*}T^{*}}{F^{*}}\partialderivative{x^{*}}\left(\log(c_{\text{e,k}}^{*})\right)\right),$
$\displaystyle\text{k}\in\\{\text{n, s, p}\\},$ (1c) $\displaystyle
I^{*}-i_{\text{e,k}}^{*}=-\sigma_{\text{k}}^{*}\partialderivative{\phi_{\text{s,k}}^{*}}{x^{*}},$
$\displaystyle\text{k}\in\\{\text{n, p}\\}.$ (1d) _Molar conservation:_
$\displaystyle\epsilon_{\text{k}}\partialderivative{c_{\text{e,k}}^{*}}{t^{*}}=\partialderivative{N_{\text{e,k}}^{*}}{x^{*}}+\frac{1}{F^{*}}\partialderivative{i_{\text{e,k}}^{*}}{x^{*}},$
$\displaystyle\text{k}\in\\{\text{n, s, p}\\},$ (1e) $\displaystyle
N_{\text{e,k}}^{*}=\epsilon_{\text{k}}^{b}D_{e}^{*}(c_{\text{e,k}}^{*})\partialderivative{c_{\text{e,k}}^{*}}{x^{*}}+\frac{t^{+}}{F^{*}}i_{\text{e,k}}^{*},$
$\displaystyle\text{k}\in\\{\text{n, s, p}\\},$ (1f)
$\displaystyle\partialderivative{c_{\text{s,k}}^{*}}{t^{*}}=\frac{1}{(r^{*})^{2}}\partialderivative{r^{*}}\left((r^{*})^{2}D_{\text{s,k}}^{*}\partialderivative{c_{\text{s,k}}^{*}}{r^{*}}\right),$
$\displaystyle\text{k}\in\\{\text{n, p}\\}.$ (1g) _Electrochemical reactions:_
$\displaystyle
j_{\text{k}}^{*}=j_{\text{$0$,k}}^{*}\sinh\left(\frac{F^{*}\eta^{*}_{\text{k}}}{2R^{*}T^{*}}\right),$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (1h) $\displaystyle
j_{\text{$0$,k}}^{*}=m_{\text{k}}^{*}(c_{\text{s,k}}^{*})^{1/2}(c_{\text{s,k,\text{max}}}^{*}-c_{\text{s,k}}^{*})^{1/2}(c_{\text{e,k}}^{*})^{1/2}$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (1i)
$\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\text{k}\in\\{\text{n, p}\\}.$ (1j) Boundary conditions
_Current:_ $\displaystyle
i^{*}_{\text{e,n}}\big{|}_{x^{*}=0}=i^{*}_{\text{e,p}}\big{|}_{x^{*}=L^{*}}=0,$
(1k)
$\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^{*},$
(1l)
$\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^{*}.$
(1m) _Concentration in electrolyte:_ $\displaystyle
N_{\text{e,n}}^{*}\big{|}_{x^{*}=0}=0,\quad
N_{\text{e,p}}^{*}\big{|}_{x^{*}=L^{*}}=0,$ (1n) $\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}}},$
(1o) $\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}}}.$
(1p) _Concentration in the electrode active material:_
$\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\text{k}\in\\{\text{n, p}\\}.$ (1q) Reference potential
$\displaystyle\phi^{*}_{\text{s,n}}\big{|}_{x^{*}=0}=0.$ (1r) Initial
conditions $\displaystyle
c_{\text{s,k}}^{*}(x^{*},r^{*},0)=c_{\text{s,k,$0$}}^{*},$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (1s) $\displaystyle
c_{\text{e,k}}^{*}(x^{*},0)=c^{*}_{\text{e,typ}},$
$\displaystyle\text{k}\in\\{\text{n, s, p}\\}.$ (1t)

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 LiPF6 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:

Global: $\displaystyle x^{*}=L^{*}x,\quad t^{*}=\tau^{*}_{\text{d}}t,\quad
D^{*}_{\text{e}}=D^{*}_{\text{e,typ}}D_{\text{e}},$ (2)
$\displaystyle\kappa^{*}_{\text{e}}=\kappa_{\text{e,typ}}^{*}\kappa_{\text{e}}\quad
I^{*}=I^{*}_{\text{typ}}I$
$\displaystyle\phi^{*}_{\text{s,n}}=\frac{R^{*}T^{*}}{F^{*}}\phi_{\text{s,n}},\quad\phi^{*}_{\text{s,p}}=\left(U^{*}_{\text{p,ref}}-U^{*}_{\text{n,ref}}\right)+\frac{R^{*}T^{*}}{F^{*}}\phi_{\text{p,k}},$
$\displaystyle\text{For }\text{k}\in\\{\text{n, p}\\}:\ $ $\displaystyle
r^{*}_{\text{k}}=R^{*}_{\text{k}}r_{\text{k}},\quad
c_{\text{s,k}}^{*}=c^{*}_{\text{s,k,max}}c_{\text{s,k}},$ $\displaystyle
j^{*}_{\text{k}}=\frac{I^{*}_{\text{typ}}}{a^{*}_{\text{k}}L^{*}}j_{\text{k}},\quad
j^{*}_{\text{$0,$k}}=\frac{I^{*}_{\text{typ}}}{a^{*}_{\text{k}}L^{*}}j_{\text{$0,$k}},\quad
m_{\text{k}}^{*}=m_{\text{k,typ}}^{*}m_{\text{k}},$
$\displaystyle\eta^{*}_{\text{k}}=\frac{R^{*}T^{*}}{F^{*}}\eta_{\text{k}},\quad
U^{*}_{\text{k}}=U^{*}_{\text{k,ref}}+\frac{R^{*}T^{*}}{F^{*}}U_{\text{k}}.\quad$
$\displaystyle\text{For }\text{k}\in\\{\text{n, p}\\}:\ $ $\displaystyle
i^{*}_{\text{s,k}}=I^{*}_{\text{typ}}i_{\text{s,k}}.$ $\displaystyle\text{For
}\text{k}\in\\{\text{n, s, p}\\}:\ $ $\displaystyle
c^{*}_{\text{e,k}}=c^{*}_{\text{e,typ}}c_{\text{e,k}},\quad{N}^{*}_{\text{e,k}}=\frac{D_{\text{e,typ}}^{*}c_{\text{e,typ}}^{*}}{L^{*}}{N}_{\text{e,k}},$
$\displaystyle\phi^{*}_{\text{e,k}}=-U^{*}_{\text{n,ref}}+\frac{R^{*}T^{*}}{F^{*}}\phi_{\text{e,k}},\quad{i}^{*}_{\text{e,k}}=I^{*}_{\text{typ}}{i}_{\text{e,k}}.$

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:

Governing equations _Charge conservation_ :
$\displaystyle\partialderivative{i_{\text{e,k}}}{x}$
$\displaystyle=\begin{cases}j_{\text{k}},\quad&\text{k}=\text{n, p},\\\
0,\qquad&\text{k}=\text{s},\end{cases}$ $\displaystyle\text{k}\in\\{\text{n,
s, p}\\},$ (3b) $\displaystyle i_{\text{e,k}}$
$\displaystyle=\epsilon_{\text{k}}^{\text{b}}\hat{\kappa}_{\text{e}}\kappa_{\text{e}}(c_{\text{e,k}})\left(-\partialderivative{\phi_{\text{e,k}}}{x}+2(1-t^{+})\partialderivative{x}\left(\log(c_{\text{e,k}})\right)\right),$
$\displaystyle\text{k}\in\\{\text{n, s, p}\\},$ (3c) $\displaystyle
I-i_{\text{e,k}}$
$\displaystyle=-\sigma_{\text{k}}\partialderivative{\phi_{\text{s,k}}}{x},$
$\displaystyle\text{k}\in\\{\text{n, p}\\}.$ (3d) _Molar conservation_ :
$\displaystyle\mathcal{C}_{\text{e}}\epsilon_{\text{k}}\gamma_{\text{e}}\partialderivative{c_{\text{e,k}}}{t}$
$\displaystyle=-\gamma_{\text{e}}\partialderivative{N_{\text{e,k}}}{x}+\mathcal{C}_{\text{e}}\partialderivative{i_{\text{e,k}}}{x},\quad$
$\displaystyle\text{k}\in\\{\text{n, s, p}\\},$ (3e) $\displaystyle
N_{\text{e,k}}$
$\displaystyle=-\epsilon_{\text{k}}^{\text{b}}D_{\text{e}}(c_{\text{e,k}})\partialderivative{c_{\text{e,k}}}{x}+\frac{\mathcal{C}_{\text{e}}t^{+}}{\gamma_{\text{e}}}i_{\text{e,k}},$
$\displaystyle\text{k}\in\\{\text{n, s, p}\\},$ (3f)
$\displaystyle\mathcal{C}_{\text{k}}\partialderivative{c_{\text{s,k}}}{t}$
$\displaystyle=\frac{1}{r_{\text{k}}^{2}}\partialderivative{r_{\text{k}}}\left(r_{\text{k}}^{2}\partialderivative{c_{\text{s,k}}}{r_{\text{k}}}\right),\quad$
$\displaystyle\text{k}\in\\{\text{n, p}\\}.$ (3g) _Electrochemical reactions_
: $\displaystyle j_{\text{k}}$
$\displaystyle=j_{\text{$0,$k}}\sinh\left(\frac{\eta_{\text{k}}}{2}\right),$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (3h) $\displaystyle
j_{\text{$0,$k}}$
$\displaystyle=\frac{\gamma_{\text{k}}}{\mathcal{C}_{\text{r,k}}}c_{\text{s,k}}^{1/2}(1-c_{\text{k}})^{1/2}c_{\text{e,k}}^{1/2}\big{|}_{r_{\text{k}}=1},$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (3i)
$\displaystyle\eta_{\text{k}}$
$\displaystyle=\phi_{\text{s,k}}-\phi_{\text{e,k}}-U_{\text{k}}(c_{\text{s,k}}\big{|}_{r_{\text{k}}=1}),$
$\displaystyle\text{k}\in\\{\text{n, p}\\}.$ (3j) Boundary conditions
_Current_ : $\displaystyle i_{\text{e,n}}\big{|}_{x=0}=0,\quad
i_{\text{e,p}}\big{|}_{x=1}=0,$ (4b)
$\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,$
(4c)
$\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.$
(4d) _Concentration in electrolyte_ : $\displaystyle
N_{\text{e,n}}\big{|}_{x=0}=0,\quad N_{\text{e,p}}\big{|}_{x=1}=0,$ (4e)
$\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}}},$
(4f) $\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}}}.$
(4g) _Concentration in the electrode active material_ :
$\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}\\}.$ (4h) _Reference potenital_ :
$\displaystyle\phi_{\text{s,n}}\big{|}_{x=0}=0.$ (4i) Initial conditions
$\displaystyle c_{\text{s,k}}(x,r,0)=c_{\text{s,k,0}},$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (4j) $\displaystyle
c_{\text{e,k}}(x,0)=1,$ $\displaystyle\text{k}\in\\{\text{n, s, p}\\}.$ (4k)

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

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

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

$\displaystyle V=$
$\displaystyle\phantom{+}\phi_{\text{s,n}}\big{|}_{x=x_{\text{n}}}-\phi_{\text{s,n}}\big{|}_{x=0}\quad$
(Negative electrode) (6)
$\displaystyle+\phi_{\text{e,n}}\big{|}_{x=x_{\text{n}}}-\phi_{\text{s,n}}\big{|}_{x=x_{\text{n}}}\quad$
(Negative electrochemical reaction)
$\displaystyle+\phi_{\text{e,p}}\big{|}_{x=x_{\text{p}}}-\phi_{\text{e,n}}\big{|}_{x=x_{\text{n}}}\quad$
(Electrolyte)
$\displaystyle+\phi_{\text{s,p}}\big{|}_{x=x_{\text{p}}}-\phi_{\text{e,p}}\big{|}_{x=x_{\text{p}}}\quad$
(Positive electrochemical reaction)
$\displaystyle+\phi_{\text{s,p}}\big{|}_{x=1}-\phi_{\text{s,p}}\big{|}_{x=x_{\text{p}}}\quad$
(Positive electrode).

We define the pointwise OCV to be

$U_{\text{eq}}\big{|}_{x_{\text{n}},x_{\text{p}}}=U_{\text{p}}(c_{\text{s,p}}\big{|}_{r=1})\big{|}_{x=x_{\text{p}}}-U_{\text{n}}(c_{\text{s,n}}\big{|}_{r=1})\big{|}_{x=x_{\text{n}}},$
(7)

and the pointwise solid-phase Ohmic losses to be

$\Delta\Phi_{\text{Solid}}\big{|}_{x_{\text{n}},x_{\text{p}}}=\left(\phi_{\text{s,p}}\big{|}_{x=1}-\phi_{\text{s,p}}\big{|}_{x=x_{\text{p}}}\right)+\left(\phi_{\text{s,n}}\big{|}_{x=x_{\text{n}}}-\phi_{\text{s,n}}\big{|}_{x=0}\right).$
(8)

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

$V=U_{\text{eq}}\big{|}_{x_{\text{n}},x_{\text{p}}}+\eta_{\text{p}}\big{|}_{x=x_{\text{p}}}-\eta_{\text{n}}\big{|}_{x=x_{\text{n}}}+\phi_{\text{e,p}}\big{|}_{x=x_{\text{p}}}-\phi_{\text{e,n}}\big{|}_{x=x_{\text{n}}}+\Delta\Phi_{\text{Solid}}\big{|}_{x_{\text{n}},x_{\text{p}}}.$
(9)

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

$V=\overline{U}_{\text{eq}}+\overline{\eta}_{\text{p}}-\overline{\eta}_{\text{n}}+\overline{\phi}_{\text{e,p}}-\overline{\phi}_{\text{e,n}}+\overline{\Delta\Phi}_{\text{Solid}}.$
(10)

where we use an overbar to represent the operation

$\overline{f}:=\frac{1}{L_{\text{p}}}\int_{1}^{1-L_{\text{p}}}\left(\frac{1}{L_{\text{n}}}\int_{0}^{L_{\text{n}}}f\,\text{d}x_{\text{n}}\right)\,\text{d}x_{\text{p}}.$
(11)

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

$\int_{0}^{L_{\text{n}}}c_{\text{e,n}}\,dx+\int_{L_{\text{n}}}^{1-L_{\text{p}}}c_{\text{e,s}}\,dx+\int_{1-L_{\text{p}}}^{1}c_{\text{e,p}}\,dx=1.$
(12)

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

$\displaystyle\sigma_{\text{k}}=\frac{\sigma_{\text{k}}^{\prime}}{\mathcal{C}_{\text{e}}},\quad\sigma_{\text{k}}^{\prime}=\mathcal{O}(1),\quad\text{k}\in\\{\text{n,
p}\\},$
$\displaystyle\hat{\kappa}_{\text{e}}=\frac{\hat{\kappa}_{\text{e}}^{\prime}}{\mathcal{C}_{\text{e}}},\quad\hat{\kappa}_{\text{e}}^{\prime}=\mathcal{O}(1).$

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

$\displaystyle c_{\text{s,k}}\sim
c_{\text{s,k}}^{0}+\mathcal{C}_{\text{e}}c_{\text{s,k}}^{1}+\mathcal{C}_{\text{e}}^{2}c_{\text{s,k}}^{2}+\dots,$

etc.

####  3.2.1 Leading-order model

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

$\partialderivative{N_{\text{e,k}}^{0}}{x}=0,\quad
N_{\text{e,k}}^{0}=-\epsilon_{\text{k}}^{b}D_{\text{e}}(c_{\text{e,k}}^{0})\partialderivative{c_{\text{e,k}}^{0}}{x},\quad\text{k}\in\\{\text{n,
s, p}\\}.$ (13)

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

$N_{\text{e,k}}^{0}=0,\quad c_{\text{e,k}}^{0}=1,\quad\text{k}\in\\{\text{n,
s, p}\\}.$ (14)

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

$\partialderivative{\phi_{\text{e,k}}^{0}}{x}=0,\quad\partialderivative{\phi_{\text{s,k}}^{0}}{x}=0.$
(15)

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

$\displaystyle i_{\text{e,n}}^{0}=\frac{xI}{L_{\text{n}}},\quad
i_{\text{e,s}}^{0}=I,\quad i_{\text{e,p}}^{0}=\frac{(1-x)I}{L_{\text{p}}},$
(16) $\displaystyle j_{\text{n}}^{0}=\frac{I}{L_{\text{n}}},\quad
j_{\text{p}}^{0}=-\frac{I}{L_{\text{p}}}.$ (17)

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

$\displaystyle\eta_{\text{n}}^{0}=2\sinh^{-1}\left(\frac{I}{j_{\text{$0$,n}}^{0}L_{\text{n}}}\right),\quad\eta_{\text{p}}^{0}=-2\sinh^{-1}\left(\frac{I}{j_{\text{0,p}}^{0}L_{\text{p}}}\right).$
(18)

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

$\phi_{\text{e,p}}^{0}-\phi_{\text{e,n}}^{0}=0,$ (19)

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

$\displaystyle\mathcal{C}_{\text{k}}\partialderivative{c_{\text{s,k}}^{0}}{t}$
$\displaystyle=\frac{1}{r_{\text{k}}^{2}}\partialderivative{r_{\text{k}}}\left(r_{\text{k}}^{2}\partialderivative{c_{\text{s,k}}^{0}}{r_{\text{k}}}\right),$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (20a)
$\displaystyle\partialderivative{c_{\text{s,k}}^{0}}{r_{\text{k}}}\bigg{|}_{r_{\text{k}}=0}$
$\displaystyle=0,\quad-\frac{a_{\text{k}}\gamma_{\text{k}}}{\mathcal{C}_{\text{k}}}\partialderivative{c_{\text{s,k}}^{0}}{r_{\text{k}}}\bigg{|}_{r_{\text{k}}=1}=\begin{cases}\frac{I}{L_{\text{n}}},\quad&\text{k}=\text{n},\\\
-\frac{I}{L_{\text{p}}},\quad&\text{k}=\text{p},\end{cases}$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (20b) $\displaystyle
c_{\text{s,k}}^{0}(r_{\text{k}},0)$ $\displaystyle=c_{\text{s,k,0}},$
$\displaystyle\text{k}\in\\{\text{n, p}\\}.$ (20c) Since $c_{\text{s,k}}^{0}$
is independent of $x$, the leading-order electrode-averaged OCV is simply
$\overline{U}_{\text{eq}}^{0}=U_{\text{p}}(c_{\text{s,p}}^{0}\big{|}_{r_{\text{p}}=1})-U_{\text{n}}(c_{\text{s,n}}^{0}\big{|}_{r_{\text{n}}=1}).$
(20d) Additionally, the leading-order electrode-averaged reaction
overpotentials are just $\overline{\eta}^{0}_{\text{k}}=\eta_{\text{k}}^{0}$.
Therefore, the leading-order voltage is given by
$V^{0}=\underbrace{\Bigg{.}U_{\text{p}}(c_{\text{s,p}}^{0})\big{|}_{r_{\text{p}}=1}-U_{\text{n}}(c_{\text{s,n}}^{0})\big{|}_{r_{\text{n}}=1}}_{\text{OCV}}\underbrace{\Bigg{.}-2\sinh^{-1}\left(\frac{I}{j_{\text{0,p}}^{0}L_{\text{p}}}\right)-2\sinh^{-1}\left(\frac{I}{j_{\text{0,n}}^{0}L_{\text{n}}}\right)}_{\text{Reaction
overpotentials}},$ (20e) where the leading-order component of the exchange-
current density, as given by (3i), is
$j_{\text{$0$,k}}^{0}=\frac{\gamma_{\text{k}}}{\mathcal{C}_{\text{r,k}}}(c^{0}_{\text{s,k}})^{1/2}(1-c_{\text{s,k}}^{0})^{1/2}.$
(20f)

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

$\displaystyle\partialderivative{N_{\text{e,k}}^{1}}{x}=\frac{1}{\gamma_{\text{e}}}\partialderivative{i_{\text{e,k}}^{0}}{x},$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (21) $\displaystyle
N_{\text{e,k}}^{1}=-\epsilon_{\text{k}}^{b}D_{\text{e}}(c^{0}_{\text{e,k}})\partialderivative{c^{1}_{\text{e,k}}}{x}+\frac{t^{+}}{\gamma_{\text{e}}}i_{\text{e,k}}^{0},$
$\displaystyle\text{k}\in\\{\text{n, p}\\}.$ (22)

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

$N_{\text{e,n}}=\frac{Ix}{\gamma_{\text{e}}L_{\text{n}}},\quad
N_{\text{e,s}}=\frac{I}{\gamma_{\text{e}}},\quad
N_{\text{e,p}}=\frac{I(1-x)}{\gamma_{\text{e}}L_{\text{p}}}.$ (23)

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

$\displaystyle
c_{\text{e,n}}^{1}=\frac{(1-t^{+})I}{\gamma_{\text{e}}6D_{\text{e}}(1)}\left(2\left(\frac{L_{\text{p}}^{2}}{\epsilon_{\text{p}}^{b}}-\frac{L_{\text{n}}^{2}}{\epsilon_{\text{n}}^{b}}\right)+\frac{3L_{s}}{\epsilon_{\text{s}}^{b}}(1+L_{\text{p}}-L_{\text{n}})+\frac{3}{\epsilon_{\text{n}}^{b}L_{\text{n}}}(L_{\text{n}}^{2}-x^{2})\right),$
(24a) $\displaystyle
c_{\text{e,p}}^{1}=\frac{(1-t^{+})I}{6\gamma_{\text{e}}D_{\text{e}}(1)}\left(2\left(\frac{L_{\text{p}}^{2}}{\epsilon_{\text{p}}^{b}}-\frac{L_{\text{n}}^{2}}{\epsilon_{\text{n}}^{b}}\right)+\frac{3}{\epsilon_{\text{s}}^{b}}(L_{\text{n}}^{2}-L_{\text{p}}^{2}+1-2x)\right),$
(24b) $\displaystyle
c_{\text{e,p}}^{1}=\frac{(1-t^{+})I}{6\gamma_{\text{e}}D_{\text{e}}(1)}\left(2\left(\frac{L_{\text{p}}^{2}}{\epsilon_{\text{p}}^{b}}-\frac{L_{\text{n}}^{2}}{\epsilon_{\text{n}}^{b}}\right)+\frac{3L_{s}}{\epsilon_{\text{s}}^{b}}(L_{\text{p}}-L_{\text{n}}-1)+\frac{3}{L_{\text{p}}\epsilon_{\text{p}}^{b}}((x-1)^{2}-L_{\text{p}}^{2})\right).$
(24c)

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

$\displaystyle
i_{\text{e,k}}^{0}=\epsilon_{\text{k}}^{b}\hat{\kappa}_{\text{e}}^{\prime}\kappa_{\text{e}}(1)\left(-\partialderivative{\phi_{\text{e,k}}^{1}}{x}+2(1-t^{+})\partialderivative{c_{\text{e,k}}^{1}}{x}\right).$
(25)

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

$\displaystyle\phi_{\text{e,n}}^{1}=\tilde{\phi}_{\text{e}}+2(1-t^{+})c_{\text{e,n}}^{1}-\frac{I}{\hat{\kappa}_{\text{e}}^{\prime}\kappa_{\text{e}}(1)}\left(\frac{x^{2}-L_{\text{n}}^{2}}{2\epsilon_{\text{n}}^{b}L_{\text{n}}}+\frac{L_{\text{n}}}{\epsilon_{\text{s}}^{b}}\right),$
(26a)
$\displaystyle\phi_{\text{e,s}}^{1}=\tilde{\phi}_{\text{e}}+2(1-t^{+})c_{\text{e,s}}^{1}-\frac{Ix}{\hat{\kappa}_{\text{e}}^{\prime}\kappa_{\text{e}}(1)\epsilon_{\text{s}}^{b}},$
(26b)
$\displaystyle\phi_{\text{e,n}}^{1}=\tilde{\phi}_{\text{e}}+2(1-t^{+})c_{\text{e,p}}^{1}-\frac{I}{\hat{\kappa}_{\text{e}}^{\prime}\kappa_{\text{e}}(1)}\left(\frac{x(2-x)+L_{\text{p}}^{2}-1}{2\epsilon_{\text{p}}^{b}L_{\text{p}}}+\frac{1-L_{\text{p}}}{\epsilon_{\text{s}}^{b}}\right),$
(26c)

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

$I-i_{\text{e,k}}^{0}=-\sigma^{\prime}_{\text{k}}\partialderivative{\phi_{\text{s,k}}^{1}}{x},\quad\text{k}\in\\{\text{n,
p}\\},$ (27)

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

$\displaystyle\phi_{\text{s,n}}^{1}=\phi_{\text{s,n}}\big{|}_{x=0}+\frac{Ix}{2\sigma^{\prime}_{\text{n}}L_{\text{n}}}\left(2L_{\text{n}}-x\right),$
(28a)
$\displaystyle\phi_{\text{s,p}}^{1}=\phi_{\text{s,p}}\big{|}_{x=1}+\frac{I(x-1)(1-2L_{\text{p}}-x)}{2\sigma^{\prime}_{\text{p}}L_{\text{p}}}.$
(28b)

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

$\displaystyle\partialderivative{i_{\text{e,k}}^{1}}{x}=j_{\text{k}}^{1},$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (29a) $\displaystyle
i_{\text{e,n}}^{1}\big{|}_{x=0,L_{\text{n}}}=i_{\text{e,p}}^{1}\big{|}_{1-L_{\text{p}},1}=0.$
(29b)

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

$\int_{0}^{L_{\text{n}}}j_{\text{n}}^{1}\,\text{d}x=0,\ \text{and}\
\int_{1-l_{\text{p}}}^{1}j_{\text{p}}^{1}\,\text{d}x=0.$ (30)

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

$\displaystyle\mathcal{C}_{\text{k}}\partialderivative{\overline{c}_{\text{s,k}}^{1}}{t}=\frac{1}{r_{\text{k}}^{2}}\partialderivative{r_{\text{k}}}\left(r_{\text{k}}^{2}\partialderivative{\overline{c}_{\text{s,k}}^{1}}{r_{\text{k}}}\right),$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (31a)
$\displaystyle\partialderivative{\overline{c}_{\text{s,k}}^{1}}{r_{\text{k}}}\bigg{|}_{r_{\text{k}}=0}=\partialderivative{\overline{c}_{\text{s,k}}^{1}}{r_{\text{k}}}\bigg{|}_{r_{\text{k}}=1}=0,$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (31b)
$\displaystyle\overline{c}_{\text{s,k}}^{1}(r_{\text{k}},0)=0,$
$\displaystyle\text{k}\in\\{\text{n, p}\\}.$ (31c)

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

$\overline{c}_{\text{s,n}}^{1}=0,\quad\overline{c}_{\text{s,p}}^{1}=0.$ (32)

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

$\overline{U}_{\text{eq}}^{1}=U_{\text{p}}^{\prime}(c_{\text{s,p}}^{0})\overline{c}_{\text{s,p}}^{1}-U_{\text{n}}^{\prime}(c_{\text{s,n}}^{0})\overline{c}_{\text{s,n}}^{1}=0.$
(33)

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

$\displaystyle
j_{\text{k}}^{1}=j_{\text{$0$,k}}^{1}\sinh\left(\frac{\eta_{\text{k}}^{0}}{2}\right)+\frac{j_{\text{$0$,k}}^{0}\eta_{\text{k}}^{1}}{2}\cosh(\frac{\eta_{\text{k}}^{0}}{2}),$
(34a)
$\displaystyle\eta_{\text{k}}^{1}=\phi^{1}_{\text{s,k}}-\phi^{1}_{\text{e,k}}-U_{\text{k}}^{\prime}(c^{0}_{\text{s,k}})c_{\text{s,k}}^{1}\big{|}_{r_{\text{k}}=1},$
(34b) $\displaystyle
j_{\text{$0$,k}}^{1}=\frac{j_{\text{$0$,k}}^{0}}{2}\left(\frac{c_{\text{s,k}}^{1}}{c_{\text{s,k}}^{0}}-\frac{c_{\text{s,k}}^{1}}{1-c_{\text{s,k}}^{0}}+c_{\text{e,k}}^{1}\right)\bigg{|}_{r_{\text{k}}=1},$
(34c)

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

$\displaystyle\overline{\eta}_{\text{k}}^{1}=-\overline{c}_{\text{e,k}}^{1}\tanh\left(\frac{\eta_{\text{k}}^{0}}{2}\right)=\frac{\overline{c}_{\text{e,k}}^{1}I}{\sqrt{(j_{0,k}^{0}L_{k})^{2}+I^{2}}}.$
(35)

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

$\displaystyle V^{1}=$
$\displaystyle\left(\frac{\overline{c}_{\text{e,p}}I}{\sqrt{\big{.}\smash[b]{(j_{\text{0,p}}^{0}L_{\text{p}})^{2}+I^{2}}}}+\frac{\overline{c}_{\text{e,n}}I}{\sqrt{\big{.}(j_{\text{0,n}}^{0}L_{\text{n}})^{2}+I^{2}}}\right)+2(1-t^{+})\left(\overline{c}_{\text{e,p}}^{1}-\overline{c}_{\text{e,n}}^{1}\right)$
(36)
$\displaystyle-\frac{I}{\hat{\kappa}_{\text{e}}^{\prime}\kappa_{\text{e}}(1)}\left(\frac{L_{\text{n}}}{3\epsilon_{\text{n}}^{b}}+\frac{L_{\text{s}}}{\epsilon_{\text{s}}^{b}}+\frac{L_{\text{p}}}{3\epsilon_{\text{p}}^{b}}\right)-\frac{I}{3}\left(\frac{L_{\text{p}}}{\tilde{\sigma}_{\text{p}}}+\frac{L_{\text{n}}}{\tilde{\sigma}_{\text{n}}}\right).$

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

$\displaystyle\overline{j}_{\text{0,n}}=\frac{1}{L_{\text{n}}}\int_{0}^{L_{\text{n}}}\frac{\gamma_{\text{n}}}{\mathcal{C}_{\text{r,n}}}(c_{\text{s,n}}^{0})^{1/2}(1-c_{\text{s,n}}^{0})^{1/2}(1+\mathcal{C}_{\text{e}}c_{\text{e,n}}^{1})^{1/2}\,\text{d}x,$
(37a)
$\displaystyle\overline{j}_{\text{0,p}}=\frac{1}{L_{\text{p}}}\int_{1-L_{\text{p}}}^{1}\frac{\gamma_{\text{p}}}{\mathcal{C}_{\text{r,p}}}(c_{\text{s,p}}^{0})^{1/2}(1-c_{\text{s,p}}^{0})^{1/2}(1+\mathcal{C}_{\text{e}}c_{\text{e,p}}^{1})^{1/2}\,\text{d}x,$
(37b)

and use the fact that

$\displaystyle-2\sinh^{-1}\left(\frac{I}{\overline{j}_{\text{0,k}}L_{\text{k}}}\right)=-2\sinh^{-1}\left(\frac{I}{j_{\text{0,k}}^{0}L_{\text{k}}}\right)+\mathcal{C}_{\text{e}}\frac{\overline{c}_{\text{e,k}}^{1}I}{\lambda\sqrt{(j_{\text{$0$,k}}^{0}L_{\text{k}})^{2}+I^{2}}}+\mathcal{O}\left(\mathcal{C}_{\text{e}}^{2}\right).$

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

$V=\overline{U}_{\text{eq}}+\overline{\eta}_{\text{r}}+\overline{\eta}_{\text{c}}+\overline{\Delta\Phi}_{\text{Elec}}+\overline{\Delta\Phi}_{\text{Solid}},$
(38)

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

$\displaystyle\overline{U}_{\text{eq}}=U_{\text{p}}(c_{\text{s,p}}^{0}\big{|}_{r_{\text{p}}=1})-U_{\text{n}}(c_{\text{s,n}}^{0}\big{|}_{r_{\text{n}}=1}),$
(39a)
$\displaystyle\overline{\eta}_{\text{r}}=-2\sinh^{-1}\left(\frac{I}{\overline{j}_{\text{p}}L_{\text{n}}}\right)-2\sinh^{-1}\left(\frac{I}{\overline{j}_{\text{n}}L_{\text{n}}}\right),$
(39b)
$\displaystyle\overline{\eta}_{\text{c}}=2\mathcal{C}_{\text{e}}(1-t^{+})\left(\overline{c}_{\text{e,p}}^{1}-\overline{c}_{\text{e,n}}^{1}\right),$
(39c)
$\displaystyle\overline{\Delta\Phi}_{\text{Elec}}=-\frac{I}{\hat{\kappa}_{\text{e}}\kappa_{\text{e}}(1)}\left(\frac{L_{\text{n}}}{3\epsilon_{\text{n}}^{b}}+\frac{L_{\text{s}}}{\epsilon_{\text{s}}^{b}}+\frac{L_{\text{p}}}{3\epsilon_{\text{p}}^{b}}\right),$
(39d)
$\displaystyle\overline{\Delta\Phi}_{\text{Solid}}=-\mathcal{C}_{\text{e}}(\overline{\varphi}_{\text{s,p}}^{1}-\overline{\phi}_{\text{s,n}}^{1}).$
(39e)

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:

Governing equations
$\displaystyle\mathcal{C}_{\text{k}}\partialderivative{c_{\text{s,k}}^{0}}{t}$
$\displaystyle=-\frac{1}{r_{\text{k}}^{2}}\partialderivative{r_{\text{k}}}\left(r_{\text{k}}^{2}\partialderivative{c_{\text{s,k}}^{0}}{r_{\text{k}}}\right),\quad$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (40b)
$\displaystyle\mathcal{C}_{\text{e}}\epsilon_{\text{k}}\gamma_{\text{e}}\partialderivative{c_{\text{e,k}}^{1}}{t}$
$\displaystyle=-\gamma_{\text{e}}\partialderivative{N_{\text{e,k}}^{1}}{x}+\begin{cases}\frac{I}{L_{\text{n}}},\quad&\text{k}=\text{n},\\\
0,\quad&\text{k}=\text{s},\\\
-\frac{I}{L_{\text{p}}},\quad&\text{k}=\text{p},\end{cases}$
$\displaystyle\text{k}\in\\{\text{n, s, p}\\},$ (40c) $\displaystyle
N_{\text{e,k}}^{1}$
$\displaystyle=-\epsilon_{\text{k}}^{\text{b}}D_{\text{e}}(1)\partialderivative{c^{1}_{\text{e,k}}}{x}+\begin{cases}\frac{xt^{+}I}{\gamma_{\text{e}}L_{\text{n}}},\quad&\text{k}=\text{n},\\\
\frac{t^{+}I}{\gamma_{\text{e}}},\quad&\text{k}=\text{s},\\\
\frac{(1-x)t^{+}I}{\gamma_{\text{e}}L_{\text{p}}},\quad&\text{k}=\text{p},\end{cases}$
$\displaystyle\text{k}\in\\{\text{n, s, p}\\}.$ (40d) Boundary conditions
$\displaystyle\partialderivative{c_{\text{s,k}}^{0}}{r_{\text{k}}}\bigg{|}_{r_{\text{k}}=0}=0,\quad-\frac{a_{\text{k}}\gamma_{k}}{\mathcal{C}_{k}}\partialderivative{c_{\text{s,k}}^{0}}{r_{\text{k}}}\bigg{|}_{r_{\text{k}}=1}=\begin{cases}\frac{I}{L_{\text{n}}}\quad&\text{k}=\text{n},\\\
-\frac{I}{L_{\text{p}}}\quad&\text{k}=\text{p},\end{cases}$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (40e) $\displaystyle
N_{\text{e,n}}^{1}\big{|}_{x=0}=0,\quad N_{\text{e,p}}^{1}\big{|}_{x=1}=0,$
(40f) $\displaystyle
c_{\text{e,n}}^{1}|_{x=L_{\text{n}}}=c_{\text{e,s}}^{1}|_{x=L_{\text{n}}},\quad
N_{\text{e,n}}^{1}\big{|}_{x=L_{\text{n}}}=N_{\text{e,s}}^{1}\big{|}_{x=L_{\text{n}}},$
(40g) $\displaystyle
c_{\text{e,s}}^{1}|_{x=1-L_{\text{p}}}=c_{\text{e,p}}^{1}|_{x=1-L_{\text{p}}},\quad
N_{\text{e,s}}^{1}\big{|}_{x=1-L_{\text{p}}}=N_{\text{e,p}}^{1}\big{|}_{x=1-L_{\text{p}}}.$
(40h) Initial conditions $\displaystyle
c_{\text{s,k}}^{0}(r_{\text{k}},0)=c_{\text{k,$0$}},$
$\displaystyle\text{k}\in\\{\text{n, p}\\},$ (40i) $\displaystyle
c_{\text{e,k}}^{1}(x,0)=0,$ $\displaystyle\text{k}\in\\{\text{n, s, p}\\}.$
(40j) Terminal voltage $\displaystyle
V=\overline{U}_{\text{eq}}+\overline{\eta}_{r}+\overline{\eta}_{c}+\overline{\Delta\Phi}_{\text{Elec}}+\overline{\Delta\Phi}_{\text{Solid}},$
(40k) where
$\displaystyle\overline{U}_{\text{eq}}=U_{\text{p}}(c_{\text{s,p}}^{0})\big{|}_{r_{\text{p}}=1}-U_{\text{n}}(c_{\text{s,n}}^{0})\big{|}_{r_{\text{n}}=1},$
(40l)
$\displaystyle\overline{\eta}_{r}=-2\sinh^{-1}\left(\frac{I}{\overline{j}_{\text{$0$,p}}L_{\text{p}}}\right)-2\sinh^{-1}\left(\frac{I}{\overline{j}_{\text{$0$,n}}L_{\text{n}}}\right),$
(40m)
$\displaystyle\overline{\eta}_{c}=2\mathcal{C}_{\text{e}}(1-t^{+})\left(\overline{c}_{\text{e,p}}^{1}-\overline{c}_{\text{e,n}}^{1}\right),$
(40n)
$\displaystyle\overline{j}_{\text{$0$,n}}=\frac{1}{L_{\text{n}}}\int_{0}^{L_{\text{n}}}\frac{\gamma_{\text{n}}}{\mathcal{C}_{\text{r,n}}}(c_{\text{s,n}}^{0})^{1/2}(1-c_{\text{s,n}}^{0})^{1/2}(1+\mathcal{C}_{\text{e}}c_{\text{e,n}}^{1})^{1/2}\,\text{d}x,$
(40o)
$\displaystyle\overline{j}_{\text{$0$,p}}=\frac{1}{L_{\text{p}}}\int_{1-L_{\text{p}}}^{1}\frac{\gamma_{\text{p}}}{\mathcal{C}_{\text{r,p}}}(c_{\text{s,p}}^{0})^{1/2}(1-c_{\text{s,p}}^{0})^{1/2}(1+\mathcal{C}_{\text{e}}c_{\text{e,p}}^{1})^{1/2}\,\text{d}x,$
(40p)
$\displaystyle\overline{\Delta\Phi}_{\text{Elec}}=-\frac{I}{\hat{\kappa}_{\text{e}}\kappa_{\text{e}}(1)}\left(\frac{L_{\text{n}}}{3\epsilon_{\text{n}}^{b}}+\frac{L_{\text{s}}}{\epsilon_{\text{s}}^{b}}+\frac{L_{\text{p}}}{3\epsilon_{\text{p}}^{b}}\right),$
(40q)
$\displaystyle\overline{\Delta\Phi}_{\text{Solid}}=-\frac{I}{3}\left(\frac{L_{\text{p}}}{\sigma_{\text{p}}}+\frac{L_{\text{n}}}{\sigma_{\text{n}}}\right).$
(40r)

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.

Figure 2: Constant-current discharge comparison of DFN model (1), SPM (20),
and SPMe (40) over a range of typical C-rates. The root mean square (RMS)
voltage error for the reduced models with respect to the DFN model is provided
in Table 4.

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.

Figure 3: Overpotential errors between the SPMe and DFN.

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

$U_{\text{k}}^{\prime\prime}(c_{\text{s,k}}^{0})\sim\mathcal{O}(\mathcal{C}_{\text{e}}^{-1}).$
(41)

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.

Figure 4: RMS voltage error relative to the DFN model with $30$ $x$-grid
points in each of the domains: negative electrode, separator, and positive
electrode; as well as $30$ $r$-grid points in the particles. This reference
model with $30$ points in each domain has an average (across C-rates) run time
of $120\text{\,}\mathrm{s}$. Here, ‘$n$ pts’ refers to $n$ points in each of
the domains: negative electrode, separator, and positive electrode as well as
$n$ points in each of the particles.

###  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})$.

Figure 5: Comparison of DFN model and SPMe internal states during a
$1\text{\,}\mathrm{C}$ constant current discharge. The DFN model solution is
given by the solid lines and the SPMe solution by the closest black dashed
line. Note that some of the black dashed lines lie upon others so it appears
that there are fewer black dashed lines than solid lines.

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.

Figure 6: Comparison of internal state across multiple C-rates at a discharge
capacity of $15\text{\,}\mathrm{A}\mathrm{h}\mathrm{/}\mathrm{m}^{2}$. The DFN
model solution is given by the solid lines and the SPMe solution by the
closest black dashed line. Note that we plot the $x$-averaged particle
concentration here instead of the surface concentration plotted in Figure 5.

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:

$V=\overline{U}_{\text{eq}}+\overline{\eta}_{\text{r}}+\eta_{\text{c}}\big{|}_{x_{\text{n}}=0,x_{\text{p}}=1}+\Delta\Phi_{\text{Elec}}\big{|}_{x_{\text{n}}=0,x_{\text{p}}=1}.$
(42)

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)

$N_{\text{e,k}}^{1}=-\epsilon_{\text{k}}^{b}D_{\text{e}}(1+\mathcal{C}_{\text{e}}c_{\text{e,k}}^{1})\partialderivative{c_{\text{e,k}}^{1}}{x}+\begin{cases}\frac{xt^{+}I}{\gamma_{\text{e}}L_{\text{n}}},\quad&\text{k}=\text{n},\\\
\frac{t^{+}I}{\gamma_{\text{e}}},\quad&\text{k}=\text{s},\\\
\frac{(1-x)t^{+}I}{\gamma_{\text{e}}L_{\text{p}}},\quad&\text{k}=\text{p},\end{cases}\quad\text{for
}\text{k}\in\\{\text{n, s, p}\\}.$ (43)

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:

$\displaystyle\eta_{\text{c}}\big{|}_{x_{\text{n}}=0,x_{\text{p}}=1}=2(1-t^{+})\log\left(\frac{1+\mathcal{C}_{\text{e}}c_{\text{e,p}}^{1}\big{|}_{x=1}}{1+\mathcal{C}_{\text{e}}c_{\text{e,p}}^{1}\big{|}_{x=0}}\right),$
(44)
$\displaystyle\Delta\Phi_{\text{Elec}}\big{|}_{x_{\text{n}}=0,x_{\text{p}}=1}=-\frac{I}{2\hat{\kappa}_{\text{e}}\overline{\kappa}_{\text{e}}^{\text{eff}}}\left(L_{\text{n}}+2L_{\text{s}}+L_{\text{p}}\right),$
(45)

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

$\Delta\Phi_{\text{Elec}}\big{|}_{x_{\text{n}}=0,x_{\text{p}}=1}=-\frac{I}{2\hat{\kappa}_{\text{e}}}\left(\frac{L_{\text{n}}}{\epsilon_{\text{n}}^{b}\overline{\kappa}_{\text{e,n}}}+2\frac{L_{s}}{\epsilon_{\text{s}}^{b}\overline{\kappa}_{\text{e,s}}}+\frac{L_{\text{p}}}{\epsilon_{\text{p}}^{b}\overline{\kappa}_{\text{e,p}}}\right),$
(46)

where

$\displaystyle\overline{\kappa}_{\text{e,n}}=\int_{0}^{L_{\text{n}}}\kappa_{\text{e}}\left(1+\mathcal{C}_{\text{e}}c_{\text{e,n}}^{1}\right)\,\text{d}x,$
(47)
$\displaystyle\overline{\kappa}_{\text{e,s}}=\int_{L_{\text{n}}}^{1-L_{\text{p}}}\kappa_{\text{e}}\left(1+\mathcal{C}_{\text{e}}c_{\text{e,s}}^{1}\right)\,\text{d}x,$
$\displaystyle\overline{\kappa}_{e,p}=\int_{1-L_{\text{p}}}^{1}\kappa_{\text{e}}\left(1+\mathcal{C}_{\text{e}}c_{\text{e,p}}^{1}\right)\,\text{d}x.$

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.

Figure 7: Comparison of versions of the SPMe: Canonical (40), Perez et. al.
[15], Han et. al. [17], and Kemper et. al. [14]. The models are compared by
considering the RMS voltage error of the simplified model voltage prediction
vs the DFN model (1).

##  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:

Governing Equations
$\displaystyle\partialderivative{c_{\text{s,k}}^{*}}{t^{*}}$
$\displaystyle=\frac{1}{(r^{*})^{2}}\partialderivative{r^{*}}\left((r^{*})^{2}D_{\text{s,k}}^{*}(c_{\text{s,k}}^{*})\partialderivative{c_{\text{s,k}}^{*}}{r^{*}}\right),\quad\text{k}\in\\{\text{n,
p}\\},$ (48b)
$\displaystyle\epsilon_{\text{k}}\partialderivative{c_{\text{e,k}}^{*}}{t^{*}}$
$\displaystyle=-\partialderivative{N_{\text{e,k}}^{*}}{x^{*}}+\begin{cases}\frac{I^{*}}{F^{*}L_{\text{n}}^{*}},\quad&\text{k}=\text{n},\\\
0,\quad&\text{k}=\text{s},\\\
-\frac{I^{*}}{F^{*}L_{\text{p}}^{*}},\quad&\text{k}=\text{p},\end{cases}\quad\text{k}\in\\{\text{n,
s, p}\\},$ (48c) $\displaystyle N_{\text{e,k}}^{*}$
$\displaystyle=-\epsilon_{\text{k}}^{\text{b}}D_{\text{e}}^{*}(c_{\text{e,typ}}^{*})\partialderivative{c^{*}_{\text{e,k}}}{x}+\begin{cases}\frac{x^{*}t^{+}I^{*}}{F^{*}L_{\text{n}}^{*}},\quad&\text{k}=\text{n},\\\
\frac{t^{+}I^{*}}{F^{*}},\quad&\text{k}=\text{s},\\\
\frac{(L^{*}-x^{*})t^{+}I^{*}}{F^{*}L_{\text{p}}},\quad&\text{k}=\text{p},\end{cases}\quad\text{k}\in\\{\text{n,
s, p}\\}.$ (48d) Boundary Conditions
$\displaystyle\partialderivative{c_{\text{s,k}}^{*}}{r^{*}}\bigg{|}_{r^{*}=0}=0,\quad-
D_{\text{s,k}}^{*}(c_{\text{s,k}}^{*})\partialderivative{c_{\text{s,k}}^{*}}{r^{*}}\bigg{|}_{r^{*}=R_{\text{k}}^{*}}=\begin{cases}\frac{I^{*}}{F^{*}a_{\text{n}}^{*}L_{\text{n}}^{*}},\quad&\text{k}=\text{n},\\\
-\frac{I^{*}}{F^{*}a_{\text{p}}^{*}L_{\text{p}}^{*}},\quad&\text{k}=\text{p},\end{cases}\quad\text{k}\in\\{\text{n,
p}\\},$ (48e) $\displaystyle N_{\text{e,n}}^{*}\big{|}_{x^{*}=0}=0,\quad
N_{\text{e,p}}^{*}\big{|}_{x^{*}=L^{*}}=0,$ (48f) $\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}}^{*}},$
(48g) $\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}}^{*}}.$
(48h) Initial Conditions $\displaystyle
c_{\text{s,k}}^{*}(r^{*},0)=c_{\text{k,$0$}}^{*},\quad\text{k}\in\\{\text{n,
p}\\},$ (48i) $\displaystyle
c_{\text{e,k}}^{*}(x^{*},0)=c_{\text{e,typ}}^{*},\quad\text{k}\in\\{\text{n,
s, p}\\}.$ (48j) Terminal Voltage $\displaystyle
V^{*}=\overline{U}_{\text{eq}}^{*}+\overline{\eta}_{r}^{*}+\overline{\eta}_{c}^{*}+\overline{\Delta\Phi}_{\text{Elec}}^{*}+\overline{\Delta\Phi}_{\text{Solid}}^{*},$
(48k) where
$\displaystyle\overline{U}_{\text{eq}}^{*}=U_{\text{p}}^{*}\left(c_{\text{s,p}}^{*}\big{|}_{r^{*}=R_{\text{p}}^{*}}\right)-U_{\text{n}}^{*}\left(c_{\text{s,n}}^{*}\big{|}_{r^{*}=R_{\text{n}}^{*}}\right),$
(48l)
$\displaystyle\overline{\eta}_{r}^{*}=-\frac{2R^{*}T^{*}}{F^{*}}\sinh^{-1}\left(\frac{I^{*}}{a_{\text{p}}^{*}\overline{j}_{\text{$0$,p}}^{*}L_{\text{p}}^{*}}\right)-\frac{2R^{*}T^{*}}{F^{*}}\sinh^{-1}\left(\frac{I^{*}}{a_{\text{n}}^{*}\overline{j}_{\text{$0$,n}}^{*}L_{\text{n}}^{*}}\right),$
(48m)
$\displaystyle\overline{\eta}_{c}^{*}=\frac{2R^{*}T^{*}}{F^{*}c_{\text{e,typ}}^{*}}(1-t^{+})\left(\overline{c}_{\text{e,p}}^{*}-\overline{c}_{\text{e,n}}^{*}\right),$
(48n)
$\displaystyle\overline{j}_{\text{$0$,n}}^{*}=\frac{1}{L_{\text{n}}^{*}}\int_{0}^{L_{\text{n}}^{*}}m_{\text{n}}^{*}(c_{\text{s,n}}^{*})^{1/2}(c_{\text{s,n,max}}^{*}-c_{\text{s,n}}^{*})^{1/2}(c_{\text{e,n}}^{*})^{1/2}\,\text{d}x^{*},$
(48o)
$\displaystyle\overline{j}_{\text{$0$,p}}^{*}=\frac{1}{L_{\text{p}}^{*}}\int_{L^{*}-L_{\text{p}}^{*}}^{L^{*}}m_{\text{p}}^{*}(c_{\text{s,p}}^{*})^{1/2}(c_{\text{s,p,max}}^{*}-c_{\text{s,p}}^{*})^{1/2}(c_{\text{e,p}}^{*})^{1/2}\,\text{d}x^{*},$
(48p)
$\displaystyle\overline{\Delta\Phi}_{\text{Elec}}^{*}=-\frac{I^{*}}{\kappa_{\text{e}}^{*}(c_{\text{e,typ}}^{*})}\left(\frac{L_{\text{n}}^{*}}{3\epsilon_{\text{n}}^{b}}+\frac{L_{\text{s}}^{*}}{\epsilon_{\text{s}}^{b}}+\frac{L_{\text{p}}^{*}}{3\epsilon_{\text{p}}^{b}}\right),$
(48q)
$\displaystyle\overline{\Delta\Phi}_{\text{Solid}}^{*}=-\frac{I^{*}}{3}\left(\frac{L_{\text{p}}^{*}}{\sigma_{\text{p}}^{*}}+\frac{L_{\text{n}}^{*}}{\sigma_{\text{n}}^{*}}\right).$
(48r)

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

$\displaystyle\text{max}\left(\left(\frac{I_{\text{typ}}^{*}L^{*}}{D_{\text{e,typ}}^{*}F^{*}c_{\text{n,max}}^{*}}\right)^{2},\left(\frac{I_{\text{typ}}^{*}L^{*}}{D_{\text{e,typ}}^{*}}\right)^{2}\left(\frac{1}{F^{*}R^{*}T^{*}}\right)\big{|}(U_{\text{k}}^{*})^{\prime\prime}\big{|}\right)$
(49)

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

$\phi_{\text{e}}^{0}=-U_{\text{n}}(c_{\text{s,n}}^{0}\big{|}_{r_{\text{n}}=1})-\eta_{\text{n}}^{0}$
(50)

$\tilde{\phi}_{\text{e}}=-2(1+t^{+})\overline{c}_{\text{e,n}}^{1}+\frac{IL_{\text{n}}}{\gamma_{\text{e}}\kappa_{\text{e}}(1)}\left(\frac{1}{3\epsilon_{\text{n}}^{b}}-\frac{1}{\epsilon_{\text{s}}^{b}}\right)+\overline{\phi}_{\text{s,n}}^{1}-\overline{\eta}_{\text{n}}^{1}$
(51)

##  Appendix B Electrode-averaged quantities

$\displaystyle\overline{c}_{\text{e,n}}^{1}=\frac{(1-t^{+})I}{6\gamma_{\text{e}}D_{\text{e}}(1)}\left(2\left(\frac{L_{\text{p}}^{2}}{\epsilon_{\text{p}}^{b}}-\frac{L_{\text{n}}^{2}}{\epsilon_{\text{n}}^{b}}\right)+\frac{2L_{\text{n}}}{\epsilon_{\text{n}}^{b}}+\frac{3L_{\text{s}}}{\epsilon_{\text{s}}^{b}}(L_{\text{p}}-L_{\text{n}}+1)\right)$
(52)
$\displaystyle\overline{c}_{\text{e,p}}^{1}=\frac{(1-t^{+})I}{6\gamma_{\text{e}}D_{\text{e}}(1)}\left(2\left(\frac{L_{\text{p}}^{2}}{\epsilon_{\text{p}}^{b}}-\frac{L_{\text{n}}^{2}}{\epsilon_{\text{n}}^{b}}\right)-\frac{2L_{\text{p}}}{\epsilon_{\text{p}}^{b}}+\frac{3L_{\text{s}}}{\epsilon_{\text{s}}^{b}}(L_{\text{p}}-L_{\text{n}}-1)\right)$
(53)
$\displaystyle\overline{\phi}_{\text{e,n}}^{1}=\tilde{\phi}_{\text{e}}+2(1-t^{+})\overline{c}_{\text{e,n}}^{1}+\frac{IL_{\text{n}}}{\hat{\kappa}_{\text{e}}^{\prime}\kappa_{\text{e}}(1)}\left(\frac{1}{3\epsilon_{\text{n}}^{b}}-\frac{1}{\epsilon_{\text{s}}^{b}}\right)$
(54)
$\displaystyle\overline{\phi}_{\text{e,p}}^{1}=\tilde{\phi}_{\text{e}}+2(1-t^{+})\overline{c}_{\text{e,p}}^{1}+\frac{IL_{\text{p}}}{\hat{\kappa}_{\text{e}}^{\prime}\kappa_{\text{e}}(1)}\left(\frac{1}{\epsilon_{\text{s}}^{b}}-\frac{1}{3\epsilon_{\text{p}}^{b}}\right)-\frac{I}{\hat{\kappa}_{\text{e}}^{\prime}\kappa_{\text{e}}(1)\epsilon_{\text{s}}^{b}}$
(55)
$\displaystyle\overline{\phi}_{\text{s,n}}^{1}=-\frac{IL_{\text{n}}}{3\sigma_{\text{n}}^{\prime}}$
(56)
$\displaystyle\overline{\varphi}_{\text{s,p}}^{1}=\frac{IL_{\text{p}}}{3\sigma_{\text{p}}^{\prime}}$
(57)
$\displaystyle\overline{\phi}_{\text{s,p}}^{1}=\overline{\varphi}_{\text{s,p}}+V^{1}$
(58)

##  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:

$\displaystyle V=$ $\displaystyle
U_{\text{eq}}+\eta_{\text{r}}+\eta_{\text{c}}+\Delta\Phi_{\text{Elec}},$ (59)
$\displaystyle U_{\text{eq}}=$ $\displaystyle
U_{\text{p}}(c_{\text{s,p}}^{0}\big{|}_{r=1})-U_{\text{n}}(c_{\text{s,n}}^{0}\big{|}_{r=1}),$
(60)
$\displaystyle\begin{split}\eta_{\text{r}}=&-2\sinh^{-1}\left(\frac{I}{L_{\text{p}}j_{\text{0,p}}(1)}\right)-2\sinh^{-1}\left(\frac{I}{L_{\text{n}}j_{\text{0,n}}(0)}\right)\\\
&+\frac{2\mathcal{C}_{\text{e}}}{\gamma_{\text{e}}}\frac{\sigma_{\text{p}}}{\overline{\kappa}_{\text{e,p}}+\frac{\mathcal{C}_{\text{e}}}{\gamma_{\text{e}}}\sigma_{\text{p}}}\left(\sinh^{-1}\left(\frac{I}{L_{\text{p}}j_{\text{0,p}}(1)}\right)-\sinh^{-1}\left(\frac{I}{L_{\text{p}}j_{\text{0,p}}(1-L_{\text{p}})}\right)\right)\\\
&+\frac{2\mathcal{C}_{\text{e}}}{\gamma_{\text{e}}}\frac{\sigma_{\text{n}}}{\overline{\kappa}_{\text{e,n}}+\frac{\mathcal{C}_{\text{e}}}{\gamma_{\text{e}}}\sigma_{\text{n}}}\left(\sinh^{-1}\left(\frac{I}{L_{\text{n}}j_{\text{0,n}}(0)}\right)-\sinh^{-1}\left(\frac{I}{L_{\text{n}}j_{\text{0,n}}(L_{\text{n}})}\right)\right),\end{split}$
(61)
$\displaystyle\begin{split}\eta_{c}=&2(1-t^{+})\left(\frac{\overline{\kappa}_{\text{e,p}}}{\overline{\kappa}_{\text{e,p}}+\frac{\mathcal{C}_{\text{e}}}{\gamma_{\text{e}}}\sigma_{\text{p}}}\log\left(\frac{1+\mathcal{C}_{\text{e}}c_{\text{e,p}}^{1}\big{|}_{x=1}}{1+\mathcal{C}_{\text{e}}c_{\text{e,p}}^{1}\big{|}_{x=1-L_{\text{p}}}}\right)+\log\left(\frac{1+\mathcal{C}_{\text{e}}c_{\text{e,s}}^{1}\big{|}_{x=L_{\text{n}}}}{1+\mathcal{C}_{\text{e}}c_{\text{e,s}}^{1}\big{|}_{x=1-L_{\text{p}}}}\right)\right.\\\
&\hphantom{2(1-t^{+})(1-t)}\left.+\frac{\overline{\kappa}_{\text{e,n}}}{\overline{\kappa}_{\text{e,n}}+\frac{\mathcal{C}_{\text{e}}}{\gamma_{\text{e}}}\sigma_{\text{n}}}\log\left(\frac{1+\mathcal{C}_{\text{e}}c_{\text{e,n}}^{1}\big{|}_{x=L_{\text{n}}}}{1+\mathcal{C}_{\text{e}}c_{\text{e,n}}^{1}\big{|}_{x=0}}\right)\right),\end{split}$
(62) $\displaystyle\Delta\Phi_{\text{Elec}}=$
$\displaystyle-\frac{I}{\hat{\kappa}_{\text{e}}}\left(\frac{L_{\text{n}}}{\overline{\kappa}_{\text{e,n}}+\mathcal{C}_{\text{e}}\gamma_{\text{e}}\sigma_{n}}-\frac{L_{s}}{\overline{\kappa}_{\text{e,s}}}+\frac{L_{\text{p}}}{\overline{\kappa}_{\text{e,p}}+\mathcal{C}_{\text{e}}\gamma_{\text{e}}\sigma_{p}}\right),$
(63)

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

$j_{\text{$0$,k}}=\frac{\gamma_{\text{k}}}{\mathcal{C}_{\text{r,k}}}(c_{\text{s,k}}^{0})^{1/2}(1-c_{\text{s,k}}^{0})^{1/2}(1+\mathcal{C}_{\text{e}}c_{\text{e,k}}^{1})^{1/2}.$

##  Appendix D Tables

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

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

Symbol | Expression | Interpretation | Value $[$\text{\,}\mathrm{s}$]$
---|---|---|---
$\tau_{\text{d}}^{*}$ | ${F^{*}c^{*}_{\text{n,max}}L^{*}}/{I_{\text{typ}}^{*}}$ | Discharge timescale | $2.2598\times 10^{4}/\mathcal{C}$
$\tau_{\text{n}}^{*}$ | ${(R_{\text{n}}^{*})^{2}}/{D^{*}_{\text{s,p}}}$ | Diffusion timescale in the negative electrode particle | $2.5641\times 10^{3}$
$\tau^{*}_{\text{p}}$ | ${(R^{*}_{\text{p}})^{2}}/{D^{*}_{\text{s,p}}}$ | Diffusion timescale in the positive electrode particle | $1\times 10^{3}$
$\tau^{*}_{\text{e}}$ | ${(L^{*})^{2}}/{D^{*}_{\text{e,typ}}}$ | Diffusion timescale in the electrolyte | 94.803
$\tau_{\text{r,n}}^{*}$ | ${F^{*}}/({m^{*}_{\text{n}}a^{*}_{\text{n}}(c^{*}_{\text{e,typ}})^{1/2}})$ | Reaction timescale in the negative electrode | 847.534
$\tau_{\text{r,p}}^{*}$ | ${F^{*}}/({m^{*}_{\text{p}}a^{*}_{\text{p}}(c^{*}_{\text{e,typ}})^{1/2}})$ | Reaction timescale in the positive electrode | $3.3901\times 10^{4}$

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

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

Table 3: Typical dimensionless parameter values. Here $\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$-direction current density of $24\text{\,}\mathrm{A}\mathrm{m}^{-2}$. This is for a cell with an initial stoichiometry of $0.8$ in the negative electrode and $0.6$ in the positive electrode with a voltage cutoff of $3.2\text{\,}\mathrm{V}$. | $0.1\text{\,}\mathrm{C}$ | $0.5\text{\,}\mathrm{C}$ | $1\text{\,}\mathrm{C}$ | $2\text{\,}\mathrm{C}$ | $3\text{\,}\mathrm{C}$
---|---|---|---|---|---
SPM | $1.72\text{\,}\mathrm{m}\mathrm{V}$ | $9.62\text{\,}\mathrm{m}\mathrm{V}$ | $19.86\text{\,}\mathrm{m}\mathrm{V}$ | $40.67\text{\,}\mathrm{m}\mathrm{V}$ | $62.78\text{\,}\mathrm{m}\mathrm{V}$
SPMe | $0.17\text{\,}\mathrm{m}\mathrm{V}$ | $1.34\text{\,}\mathrm{m}\mathrm{V}$ | $3.04\text{\,}\mathrm{m}\mathrm{V}$ | $7.36\text{\,}\mathrm{m}\mathrm{V}$ | $13.34\text{\,}\mathrm{m}\mathrm{V}$
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. | 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 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$ pts’ refers to $n$ points in each domain. Parameter combination | Required size |  Interpretation
---|---|---
$\mathcal{C}_{\text{e}}=I_{\text{typ}}^{*}L^{*}/(D_{\text{e,typ}}^{*}F^{*}c_{\text{n,max}}^{*})$ | $\ll 1$ |  Lithium-ion migration timescale is small relative to typical discharge timescale
$R^{*}T^{*}\sigma_{\text{k}}^{*}/(F^{*}I_{\text{typ}}^{*}L^{*})$ | $\gg 1$ |  Thermal voltage is large relative to the typical potential drop in electrode k
$R^{*}T^{*}\kappa_{\text{e,typ}}^{*}/(F^{*}I_{\text{typ}}^{*}L^{*})$ | $\gg 1$ |  Thermal voltage is large relative to the typical potential drop in the electrolyte
$(R_{\text{k}}^{*})^{2}I_{\text{typ}}^{*}/(D_{\text{s,k}}^{*}F^{*}c_{\text{n,max}}^{*}L^{*})$ | $\ll 1/\mathcal{C}_{\text{e}}$ |  Solid diffusion occurs on a shorter or similar timescale to a discharge
$I_{\text{typ}}^{*}/(m_{\text{k}}^{*}a_{\text{k}}^{*}(c_{\text{e,typ}}^{*})^{1/2}c_{\text{n,max}}^{*}L^{*})$ | $\ll 1/\mathcal{C}_{\text{e}}$ |  Reactions occur on a shorter or similar timescale to a discharge
Table 6: The key conditions to be satisfied for the application of (48). In
addition, it is required that $\mathcal{C}_{\text{e}}\ll
L_{\text{k}}^{*}/L^{*}\ll 1/\mathcal{C}_{\text{e}}$,
$\mathcal{C}_{\text{e}}\ll c_{\text{p,max}}/c_{\text{n,max}}\ll
1/\mathcal{C}_{\text{e}}$, and $\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.

## References

  * [1] M. Doyle, T. F. Fuller, J. Newman, Modelling of galvanostatic charge and discharge of the lithium/polymer/insertion cell, Journal of the Electrochemical Society 140 (6) (1993) 1526–1533.
  * [2] T. F. Fuller, Simulation and Optimization of the Dual Lithium-Ion Insertion Cell, Journal of The Electrochemical Society 141 (1) (1994) 1.  doi:10.1149/1.2054684.
URL http://jes.ecsdl.org/cgi/doi/10.1149/1.2054684

  * [3] J. Newman, K. E. Thomas-Alyea, Electrochemical systems, John Wiley & Sons, 2012\.
  * [4] A. M. Bizeray, S. Zhao, S. R. Duncan, D. A. Howey, Lithium-ion battery thermal-electrochemical model-based state estimation using orthogonal collocation and a modified extended Kalman filter, Journal of Power Sources 296 (2015) 400–412.  arXiv:1506.08689, doi:10.1016/j.jpowsour.2015.07.019.
  * [5] L. Cai, R. E. White, Lithium ion cell modeling using orthogonal collocation on finite elements, Journal of Power Sources 217 (2012) 248–255.  doi:10.1016/j.jpowsour.2012.06.043.
URL http://linkinghub.elsevier.com/retrieve/pii/S0378775312010439

  * [6] Y. Zeng, P. Albertus, R. Klein, N. Chaturvedi, A. Kojic, M. Z. Bazant, J. Christensen, Efficient conservative numerical schemes for 1d nonlinear spherical diffusion equations with applications in battery modeling, Journal of The Electrochemical Society 160 (9) (2013) A1565–A1571.
  * [7] R. N. Methekar, V. Ramadesigan, J. C. Pirkle, V. R. Subramanian, A perturbation approach for consistent initialization of index-1 explicit differential-algebraic equations arising from battery model simulations, Computers and Chemical Engineering 35 (11) (2011) 2227–2234.  doi:10.1016/j.compchemeng.2011.01.003.
  * [8] V. R. Subramanian, V. Boovaragavan, V. Ramadesigan, M. Arabandi, Mathematical Model Reformulation for Lithium-Ion Battery Simulations: Galvanostatic Boundary Conditions, Journal of The Electrochemical Society 156 (4) (2009) A260.  doi:10.1149/1.3065083.
URL http://jes.ecsdl.org/cgi/doi/10.1149/1.3065083

  * [9] P. W. C. Northrop, V. Ramadesigan, S. De, V. R. Subramanian, Coordinate Transformation, Orthogonal Collocation, Model Reformulation and Simulation of Electrochemical-Thermal Behavior of Lithium-Ion Battery Stacks, Journal of The Electrochemical Society 158 (12) (2011) A1461.  doi:10.1149/2.058112jes.
URL http://jes.ecsdl.org/cgi/doi/10.1149/2.058112jes

  * [10] S. J. Moura, F. B. Argomedo, R. Klein, A. Mirtabatabaei, M. Krstic, Battery State Estimation for a Single Particle Model With Electrolyte Dynamics, IEEE Transactions on Control Systems Technology 25 (2) (2017) 453–468.  doi:10.1109/tcst.2016.2571663.
  * [11] A. M. Bizeray, J.-H. Kim, S. R. Duncan, D. A. Howey, Identifiability and parameter estimation of the single particle lithium-ion battery model, IEEE Transactions on Control Systems Technology (99) (2018) 1–16.
  * [12] M. Guo, G. Sikha, R. E. White, Single-particle model for a lithium-ion cell: Thermal behavior, Journal of The Electrochemical Society 158 (2) (2011) A122–A132.
  * [13] S. Dey, B. Ayalew, P. Pisu, Combined estimation of state-of-charge and state-of-health of li-ion battery cells using smo on electrochemical model, in: 2014 13th International Workshop on Variable Structure Systems (VSS), IEEE, 2014, pp. 1–6.
  * [14] P. Kemper, D. Kum, Extended single particle model of li-ion batteries towards high current applications, in: Vehicle Power and Propulsion Conference (VPPC), 2013 IEEE, IEEE, 2013, pp. 1–6.
  * [15] H. E. Perez, X. Hu, S. J. Moura, Optimal Charging of Batteries via a Single Particle Model with Electrolyte and Thermal Dynamics (10) (2016) 4000–4005.
  * [16] E. Prada, D. Di Domenico, Y. Creff, J. Bernard, V. Sauvant-Moynot, F. Huet, Simplified electrochemical and thermal model of lifepo4-graphite li-ion batteries for fast charge applications, Journal of The Electrochemical Society 159 (9) (2012) A1508–A1519.
  * [17] X. Han, M. Ouyang, L. Lu, J. Li, Simplification of physics-based electrochemical model for lithium-ion battery on electric vehicle. part ii: Pseudo-two-dimensional model simplification and state of charge estimation, Journal of Power Sources 278 (2015) 814–825.
  * [18] S. K. Rahimian, S. Rayman, R. E. White, Extension of physics-based single particle model for higher charge–discharge rates, Journal of Power Sources 224 (2013) 180–194.
  * [19] T. R. Tanim, C. D. Rahn, C.-Y. Wang, A temperature dependent, single particle, lithium-ion cell model including electrolyte diffusion, Journal of Dynamic Systems, Measurement, and Control 137 (1) (2015) 011005.
  * [20] E. Hinch, Perturbation methods. 1991.
  * [21] C. M. Bender, S. A. Orszag, Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory, Springer Science & Business Media, 2013.
  * [22] G. Richardson, G. Denuault, C. P. Please, Multiscale modelling and analysis of lithium-ion battery charge and discharge, Journal of Engineering Mathematics 72 (1) (2012) 41–72.  doi:10.1007/s10665-011-9461-9.
URL http://link.springer.com/10.1007/s10665-011-9461-9

  * [23] V. Sulzer, S. J. Chapman, C. P. Please, D. A. Howey, C. W. Monroe, Faster lead-acid battery simulations from porous-electrode theory: Ii. asymptotic analysis, arXiv preprint arXiv:1902.01774 (2019).
  * [24] I. R. Moyles, M. G. Hennessy, T. G. Myers, B. R. Wetton, Asymptotic reduction of a porous electrode model for lithium-ion batteries, arXiv preprint arXiv:1805.07093 (2018).
  * [25] J. Foster, I. Korotkin, M. Castle, G. Richardson, Generalized single particle model for high-rate operation of graded lithium-ion electrodes: systematic derivation and validation., preprint (2019).
  * [26] J. Newman, Dualfoil, http://www.cchem.berkeley.edu/jsngrp/fortran.html (2014).
  * [27] S. Moura, fastdfn, https://github.com/scott-moura/fastDFN (2016).
  * [28] A. M. Bizeray, J.-H. Kim, S. R. Duncan, D. A. Howey, Identifiability and parameter estimation of the single particle lithium-ion battery model (2017).  arXiv:1702.02471.
URL http://arxiv.org/abs/1702.02471

  * [29] A.C. Hindmarsh, P.N. Brown, K.E. Grant, S.L. Lee, R. Serban, D.E. Shumaker, C.S. Woodward, Suite of nonlinear and differential/algebraic equation solvers, ACM Transactions on Mathematical Software (TOMS) 31 (3) (2005) 363–396.
  * [30] pybamm-team, PyBaMM, https://github.com/pybamm-team/PyBaMM (2019).
  * [31] Reniers, J., Howey, D., Mulder, G. Review and performance comparison of mechanical-chemical degradation models for lithium-ion batteries. (2019)

## 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.