Systematic derivation and validation of a reduced thermal-electrochemical model for lithium-ion batteries using asymptotic methods

The widely used Doyler-Fuller-Newman (DFN) model for lithium-ion batteries is too computationally expensive for certain applications, which has motivated the appearance of a plethora of simpler models. These models are usually posed in an ad hoc manner, leading to inconsistencies with the DFN model and to multiple formulations of the same model, and the Single Particle Model (SPM) is a very good example of the latter. In this work, we discuss the concept of SPM-type models showing that, despite the multiple formulations found in the literature, these models always follow the same structure, and we extend this discussion to models accounting for thermal effects. Then, we present a Thermal Single Particle Model with electrolyte (TSPMe) derived in a systematic manner using asymptotic techniques. The validation of the TSPMe against a thermal DFN model shows very high accuracy with a computational cost almost forty times smaller. The comparison against experimental data shows that the model does a reasonable job predicting the behaviour of a real battery, but a very good parameter set is required to obtain accurate predictions.


Introduction
With the electrification of vehicles and the spread of portable electronic devices, lithium-ion batteries have become a key technology for energy storage. This has motivated a quest to design batteries that can store more energy and last for longer, and to manage them efficiently and safely. Mathematical models are an invaluable tool for battery design and control as they provide a cheap, safe and fast alternative to experiments.
There exist multiple approaches to deterministic models for lithium-ion batteries, and they can usually be classified into two broad categories: equivalent-circuit models and physics-based models. Equivalentcircuit models describe the battery by assuming it has the structure of a particular electrical circuit. The morphology of this circuit is a modelling assumption, and the values of the different parameters need to be determined from fitting the model to experimental data (see [1,2] for examples). Equivalent-circuit models are widely used in battery management systems (BMS) because they are computationally cheap to simulate [3]. However, they offer little insight on the battery internal states and therefore it is very complicated to extend these models to account for additional physics. A detailed description of equivalent-circuit models can be found in [4].
In this work, we focus on physics-based models. These models are derived from fundamental physical laws and have the advantage that they provide a lot of insight on the internal states of the battery. However, due to their complexity, they are computationally more expensive than equivalent-circuit models. The most widely used physics-based model is the Doyle-Fuller-Newman (DFN) model, originally posed in [5,6].
This model accounts for conservation of mass and charge in both electrode and electrolyte, and the four equations are coupled together through the reaction current density term that describes the electrodeelectrolyte intercalation reaction. One of the key aspects for the success of the DFN model is that it assumes a simplified geometry which captures the main features of a battery but is much more affordable to solve than trying to resolve all the phenomena at a microscale level [7]. The geometry of the DFN model assumes that most of the variables (electrolyte concentration and potentials both in electrode and electrolyte) can be resolved at a macroscale level, in which the material is assumed to be homogeneous. Then, the effects of the porous material on these equations are accounted for by taking effective parameters. The transport of intercalated lithium is described by assuming that at each point of the electrode there is a representative spherical particle in which lithium diffuses. Even though it is not explicitly derived in [5], this is a homogenised model and it can be formalised from the mathematical point of view (see [8] for details). A more detailed description of the DFN model can be found in the handbooks [4,9].
The DFN model is still quite complex as it involves a coupled system of differential-algebraic equations (DAEs). Hence, this model is not suitable for applications where speed is crucial, such as battery optimization and control, so simpler models have been posed in the literature. One of the best-known examples is the Single Particle Model (SPM). The key idea of this model is that we can use one single particle in each electrode to represent the behaviour of all particles. SPMs may [10] or may not [11] include electrolyte dynamics. In both cases the model can be drastically simplified so all the remaining variables of interest can be easily computed from the concentrations in the electrode particles and the electrolyte.
The SPMs are widely used, but in their origin there are still two main challenges. First, most of these models are posed on an ad hoc basis and therefore the connection to the DFN (or similar) model is lost. This is far from ideal because, in the process of derivation, inconsistencies might be introduced and certain features lost. In addition, it makes it hard to add new physics in a consistent way. The second issue, closely related to the first, is that we can find multiple Single Particle Models in the literature following different formulations [10,11,12,13,14,15].
Over the last few years, there has been a surge in the asymptotic analysis of battery models [14,15,16,17,18,19,20] with the goal of obtaining reduced models in a systematic manner. The different authors take different assumptions to reduce the models. One big difference is that [17,18,19,20] assume fast diffusion in the electrodes and thus do not reduce to SPM-type models, while [14,15,16] retain this feature. Given that our interest on this work is on SPM-type models, we will focus our attention on the latter two.
The key idea behind asymptotic methods is to analyse the dimensionless groupings that appear in the model, determine which can be assumed to be small and perform an expansion in the limit where these parameters tend to zero to obtain a simpler model. More details on these methods can be found in the handbook [21]. The main advantage of these methods is that they offer a systematic approach to model reduction, so they can be applied to very complex models, and they ensure consistency between the full and reduced models, allowing for the assumptions to be validated and the error estimated a priori (i.e. before solving the model).
So far we referred to electrochemical models only, but a similar description holds for coupled thermalelectrochemical models. These are usually posed as a coupling between the DFN model, to describe the electrochemistry, and some sort of thermal model [4,22,23]. Due to their complexity, simpler models have been posed ad hoc, usually based on SPM-type electrochemical models [13]. The asymptotic analysis of thermal-electrochemical models is even more recent as it is a natural evolution of the asymptotic reduction of pure electrochemical models. To our knowledge, the only references in the literature on asymptotic reductions of thermal-electrochemical models are [20,24]. In [20] a multilayer cell is considered but with the electrode fast diffusion electrochemical model presented in [18]; while in [24] a model for single-layer pouch cells based on DFN and the SPMe in [14].
The goal of this paper is to derive a SPM-type thermal-electrochemical model for a multi-layer cell, using an asymptotic method based on minimal assumptions. The resultant model, which we refer to as Thermal Single Particle Model with electrolyte (TSPMe), is validated against the full thermal DFN model and experimental data to show its validity and proof that it is suitable for practical applications.  Each block is a model (or part of a model) that converts the input variable into the output variable (or variables). The colour of the blocks indicates both the type of model and its complexity: red for a system of differential-algebraic equations (DAEs, high complexity), yellow for a system of partial differential equations (PDEs, medium complexity) and green for explicit expressions (low complexity).

Discussion of Single Particle Models
Before presenting the formulation for the Thermal Single Particle Model with electrolyte, we discuss the nature of Single Particle Models. Some of the most popular reduced models for lithium-ion batteries are Single Particle Models (SPM), so in the literature we can find several different models referred to as such [10,11,12,13,14,15]. These are usually based on slightly different sets of assumptions and therefore their formulations are different as well. Because the different models are referred to by the same name, it makes the term "Single Particle Model" very confusing. With the hope of sheding some light into this issue, in this section we discuss what are the fundamental features that constitute a Single Particle Model, understanding it as a type of model rather than a specific formulation. We refer to the models of this kind as SPM-type models, which applies to both models with and without electrolyte dynamics. However, given that the model without electrolyte dynamics (SPM) is a particular case to the model with electrolyte dynamics (SPMe), we consider the latter in this analysis. By considering what are the essential features of an SPM-type model we can better understand and classify existing models in the literature.

Electrochemical models
We start focusing on purely electrochemical models. By electrochemical models we refer to models that only account for mass and charge transport (in both the electrodes and the electrolyte) but do not include further physical effects such as thermal, mechanical or degradation. The classic example of electrochemical model is the Doyle-Fuller-Newman (DFN) model [5], which we use as a starting point in our model reduction (see Appendix B for details).
To better understand these models, we can represent them as a block diagram, like the ones shown in Figure 1. In these diagrams, we have an input (or inputs) that are known and want to calculate an output which is unknown. To convert an input into an output we use a model, represented in the diagram by a block. Note that this model could be anything that converts the input into an output, but here we consider it from the mathematical point of view, in which the model is composed by a set of differential and/or algebraic equations that need to be solved, without entering into the method used to solve them.
The most typical setup for electrochemical models is to use the applied current as input and calculate the terminal voltage as an output. Other setups, such as inputting the voltage are common as well, but make the reduced models significantly more expensive as we discuss later. In addition, sometimes there are other variables that might be of interest and hence be part of the model output.
Before discussing the advantages of SPM-type models, let's consider the DFN model [5]. The model is composed by a system of differential-algebraic equations (DAEs). Then, in the block diagram it is represented as a single block so, to calculate the output, we need to solve the full model in a single step (see Figure 1a). Solving numerically a system of DAEs is a complex task, hence the model is very expensive to compute.
The key aspect of SPM-type models is that they break this block into simpler blocks (see Figure 1b), and we can group these blocks into two main steps. The first step involves calculating from the input current the concentrations in the representative electrode particles, and in the electrolyte (if not neglected). Then, all the other variables of interest can be calculated from these concentrations using explicit expressions, or even closed-form expressions in some cases.
The first step in the SPM-type models (yellow blocks in Figure 1b) involves solving up to three decoupled partial differential equations: two PDEs for the concentration of lithium in the particles (one for each electrode) and one PDE for the ion concentration in the electrolyte. If the electrolyte dynamics are not considered, then we only need to solve the two PDEs for the particles. Because the PDEs are decoupled, we can solve them independently which results into a relatively simple problem. These PDEs can be solved numerically using the method of lines combined with standard spatial methods (e.g. finite volume methods, finite element methods...) and further techniques can be used to speed up the calculations of the solution to the resulting ODE system (e.g residue grouping [25], balanced truncation [26]...).
The second step (green block in Figure 1b) consists on evaluating the explicit expression for the terminal voltage, and for other variables of interest (such as currents or potentials) if required. which takes the electrolyte and electrode concentrations as inputs. Because this step involves evaluating an explicit expression rather than solving an equation, it is computationally cheap. In addition, we can limit this step to evaluating only the variables we are interested in.
The main advantage of this approach is that it fully decouples the problem, so each concentration is the solution of an independent PDE and therefore a lot cheaper to solve than the DFN model. Note, however, that this decoupling only occurs in the case where current is the input of the model. If, instead, we want to simulate a potentiostatic discharge, we need to impose an algebraic constraint on the voltage and, thus, we need to solve a system of DAEs. This system, even though it is more complex than the standard decoupled SPM-type models, is still a lot simpler than the DFN model.
In the literature, we find several instances of SPM-type models, posed either in an ad hoc manner [10,11,12,13] or a systematic manner [14,15]. Even though in most cases both approaches yield very similar results, asymptotic methods are systematic and, therefore, more reliable. Moreover, the systematic approach allows us to easily extend the models to account for extra physical phenomena.
These systematically derived models are obtained by applying asymptotic techniques, in which the size of the dimensionless groupings appearing in the model is used to decide which terms of the equations can be neglected. There are different sets of assumptions that lead to the derivation of SPM-type models (see [14,15] for two different examples) but they yield very similar results [15,27]. However, all these different sets of assumptions are based on the same key idea that allows for SPM-type models: the battery operates in a regime where deviations from the equilibrium potential are small, which means that the battery operates at moderate to low C-rates.
In this scenario, the analysis is as follows. If the deviations from equilibrium potential are small, then all the particles within an electrode are approximately at the same potential, which is around the equilibrium (or open-circuit) potential, and therefore they behave very similarly. In consequence, we can define a representative particle (or average particle) in each electrode on which to solve the lithium diffusion problem. Because all the particles behave similarly, the reaction current density for the lithium intercalation must be identical for all particles within an electrode and thus, using a conservation of charge argument, we can calculate it directly from the input current. This is a crucial result because the reaction current density was the condition that coupled the four conservation equations in the DFN model (conservation of mass and charge in both electrode and electrolyte). If the reaction current density is known a priori, each of the four problems can now be treated separately and, moreover, the potential equations both in the electrodes and the electrolyte can be integrated analytically. Therefore, the only equations to solve  Each block is a model (or part of a model) that converts the input variable into the output variable (or variables). The colour of the blocks indicates both the type of model and its complexity: red for a system of differential-algebraic equations (DAEs, high complexity), yellow for a system of partial differential equations (PDEs, medium complexity) and green for explicit expressions (low complexity).
numerically are the ones governing the concentrations and everything else can be calculated from explicit expressions, as shown in Figure 1b. This reasoning can be formalised into a rigorous mathematical analysis using asymptotic techniques, as shown in Appendix B.2. If additional assumptions are taken, there might be further simplifications that take place in the model which may yield different expressions for the potentials (green block in Figure 1b). These further simplifications are discussed in Section 3.1.

Thermal-electrochemical models
After considering electrochemical SPM-type models, we focus on the coupled thermal-electrochemical models. These models, apart from the physics discussed in the previous section, describe how the temperature evolves in the battery. The two models are strongly coupled because the parameters of the electrochemical model depend on temperature and the heat generated in the battery depends on the electrochemical processes. This coupling adds an extra layer of complexity to thermal-electrochemical models compared to pure electrochemical models.
We can visualise this extra layer of complexity in the form of a block diagram, as shown in Figure 2. Compared to the electrochemical model in Figure 1, the thermal-electrochemical model has additional blocks to process the thermal inputs and outputs. In particular, the new inputs are current and ambient temperature, and the new outputs are voltage and cell temperature. Depending on the model used, the cell temperature can be the average temperature, the surface temperature or a temperature distribution across the cell.
The classic example of a thermal-electrochemical model would be the DFN coupled with a thermal model as shown in [4,22,23]. This model is even more computationally expensive than the standard DFN model, and another issue arises when considering a battery composed of multiple layers (we refer to each layer as a "cell"). As discussed in [8], in this situation the electrochemical behaviour needs to be described at a cell level, but the thermal behaviour needs to be described at a battery level (i.e. across multiple cells). This brings up a multi-scale coupling in the model that requires a hierarchy of submodels, similar to the one between electrodes and particles in the DFN model.
The thermal model is posed at the battery level because heat can flow from one cell to the neighbouring ones. However, ions cannot flow from one cell to the neighbouring ones, so the electrochemical models should battery (∼ 1 cm) cell (∼ 100 µm) particle (∼ 1 µm) thermal model cell model particle model be posed independently for each cell. Even in the case for cylindrical and prismatic cells, composed of big electrodes rolled up around a mandrel, the path for ions to go from one layer to the neighbouring ones is so long compared to the electrode thickness that these effects can be neglected. The two problems are strongly coupled, because the electrochemical properties of each cell can depend on the temperature while the heat generated by the electrochemistry impacts the evolution of temperature. The coupling between these two problems at different scales imposes a hierarchy similar to the classic DFN model. There, we solve for the electrode potential and the electrolyte concentration and potential across the cell, and at each point of the electrode domain there is a representative particle that describes the intercalated lithium concentration at the microscale. Here, we solve for temperature across the battery and at each point there is a representative cell that describes the electrochemical behaviour. Moreover, if the chosen electrochemical model is the DFN, then the hierarchy is at three levels: battery, cell and particle. These relations are shown in Figure 3.
This thermal-electrochemical coupling at multiple scales results into a system of DAEs again, even when using an SPM-type electrochemical model, which was one of the main causes of the complexity of the DFN model. Therefore, in this case too, we look for methods to simplify the model. The starting point is the DFN model coupled to a thermal model (see [4,22,23] for details). The simplifications on the electrochemical part of the model follow the procedure explained in Section 2.1, so here we focus on the thermal part only. In order to get rid of the DAEs again, we need to work in a regime in which the different instances of the electrochemical model (i.e. the different cells) can be considered to be very similar. It is the same idea as taking that all the particles work similarly within each electrode in the SPM model. There are two main assumptions that lead to this decoupling: the assumption that the temperature variation is small enough to not cause significant variations in the electrochemical parameters, and the assumption that the temperature is approximately homogeneous within the battery (see [24]). The validity of each assumption depends on the battery geometry and the cooling conditions. In this work, we focus on the latter assumption as it is reasonable for the batteries and the conditions we are interested in: LG M50 cylindrical batteries cycled in a thermal chamber without any active cooling. The full derivation of the reduced thermal model is shown in Appendix B.1.
Then, in terms of the block diagram, the TSPMe works in the following way. A system of coupled PDEs and ODEs (the latter for the temperature) is solved from the input variables, which in this case are the input current and ambient temperature. Then, all the other quantities of interest can be calculated from them. The main difference compared to the purely electrochemical models is that, if the parameters depend on temperature, the system of PDEs and ODEs is now coupled so concentrations and temperature need to be solved simultaneously (see Figure 2).

Thermal Single Particle Model with electrolyte (TSPMe)
Having discussed the background and interpretation of SPM-type models, we can now present the Thermal Single Particle Model with electrolyte (TSPMe). The details on the full derivation from the DFN model are given in Appendix B.
The TSPMe is a reduced thermal-electrochemical model for lithium-ion batteries. It consists of an SPMe for the electrochemistry coupled with a lumped thermal model for the cell. It requires solving three diffusion equations for the concentrations (two for a representative particle in each electrode and one for the electrolyte) and a first-order ODE for the average cell temperature. This model can be further simplified for some particular cases, as presented in Section 3.1.
The geometry of the model is shown in Figure 4. The two electrode particles are defined in the domain 0 < r < R k (for k ∈ {n, p}), where n denotes the negative electrode/particle and p denotes the positive electrode/particle. The electrolyte domain spans across the two porous electrodes and the separator (0 < x < L), where the negative electrode is defined in 0 < x < L n , the separator in L n < x < L − L p and the positive electrode in L − L p < x < L. For the full cell, given that the heat equation is averaged over the whole cell, we do not need to define a domain.
We need to solve for the concentrations (and we are interested in certain quantities defined) both in the electrodes and the electrolyte, so we identify them with the subscripts n and p to denote the negative and positive electrodes/particles, respectively, and the subscript e to denote the electrolyte. Then, in each particle the governing equations for intercalated lithium are where c k is the intercalated lithium concentration, D k is the diffusion coefficient of intercalated lithium (which may depend on c k ), R k is the particle radius, J k is the reaction volumetric current density, a k is the surface area density, F is the Faraday constant, c k,init is the initial concentration, i app is the applied current density, and L k is the electrode thickness. Recall that the subscript k ∈ {n, p} denotes the negative and positive electrode, respectively. The electrolyte problem is given by where c e is the lithium-ion concentration in the electrolyte, D e is the bulk diffusion coefficient of lithium ions (which may depend on c e ), B(x) is the geometry factor, t + is the transference number, L is the total cell thickness, and c e,init is the initial concentration. In (3) the geometric factor B(x) and the porosity (or electrolyte volume fraction) ε(x) are defined as where the values in each subdomain B k and ε k are constant. The thermal problem is given by with where T is the average temperature of the cell (assumed to be homogeneous in space), θ is the lumped volumetric heat capacity of the cell, h is the heat exchange coefficient, a cool is the cooling surface area of the cell per unit of volume, T amb is the ambient temperature, Q s is the heat source term due to the electrode (Joule heating), Q e is the heat source term due to the electrolyte (both Joule heating and due to concentration gradients), Q irr is the irreversible heat source term of the intercalation reaction and Q rev is the reversible heat source term of the intercalation reaction. In addition, σ k is the electronic conductivity of the electrodes, σ e is the ionic conductivity of the electrolyte, R is the universal gas constant, j k is the reaction surface current density, and Π k is the Peltier term. After solving (1), (3) and (5) numerically, we can use the values of c k , c e and T to calculate any other variable of interest from explicit expressions.
For example, the potentials and current in the electrodes and the electrolyte can be calculated as The exchange current densities j n and j p are defined as Here, U k is the open circuit potential, m k is the intercalation reaction rate, and c max k is the maximum concentration in the electrode.
The output voltage is given by where We can interpret these terms in the following way: U eq is the equilibrium potential to which the terminal voltage converges when no current is applied, while all the other terms are deviations from equilibrium due to different effects. η r is due to the reaction overpotentials, η e is due to the concentration gradients in the electrolyte, ∆Φ e is due to Ohmic losses in the electrolyte, and ∆Φ s is due to Ohmic losses in the (solid) electrode.
The TSPMe is valid in a particular regime both in terms of electrochemical and thermal submodels. The electrochemical submodel is valid when the deviations from the equilibrium potential are slow, which is equivalent to moderate to low C-rates. In terms of the dimensionless parameters, this corresponds to the case where where i 0 and Φ 0 are the typical current and electrode potential, respectively, and σ e,typ is the typical value of σ e . Further details on these dimensionless parameters can be found in Appendix B. The thermal submodel is valid when the heat exchange with environment is much slower than the heat transfer inside of the battery. In terms of dimensionless parameters, this corresponds to the case where where L batt is the typical length scale of the battery, κ is the lumped thermal conductivity of the battery and t 0 is the typical discharge time.

Further simplifications
Notice that the TSPMe is already an SPM-type model, as defined in Section 2, which can be implemented at a very low computational cost. However, some of the expressions, especially those related with the electrolyte, are a bit convoluted and this can make them not suitable for certain applications and hard to grasp their physical meaning. In this section we present some further simplifications that can be done by taking additional assumptions, and the full details of the derivation can be found in Appendix C.
Here we present the dimensional form of three main different simplifications: quasi-steady-state electrolyte concentration, constant electrolyte conductivity and fast lithium diffusion. For the first two we consider as well a particular case that allows even further reduction of the model. Notice that the three simplifications are independent from each other, so we can choose to use a subset of them or all together.

Quasi-steady state electrolyte concentration
One simplification for the electrolyte concentration is to take the quasi-steady-state problem, which is a valid assumption when the current varies over a larger time scale than the electrolyte diffusion and electrolyte ion generation time scales. Mathematically, this corresponds to the case where The physical meaning of this limit is that electrolyte transient effects are negligible because they happen at a much shorter time scale than variations in the applied current. This is a reasonable assumption for constant-current discharge (or charge) experiments. With these assumptions, there is no time dependence in the electrolyte concentration equations so their integration is a one-off step. Additionally, if the ion transport properties in the electrolyte (D e and t + ) are assumed to be constant, we can find the following analytical expression for the concentration where the volume of pores per unit of electrode plate area is given by and L s = L − L n − L p is the thickness of the separator. An additional step to the previous simplification occurs when the variation in electrolyte concentration is small. This corresponds to the limit and it yields the result c e = c e,init . Therefore, we obtain a Single Particle Model as we do not need to solve for c e . Notice that this case, even though very similar to the one studied in [14], is not identical as we have taken the additional assumption λ 1. This assumption also helps simplifying the expression for the electrolyte potential because we can take σ e (c e ) as a constant. We provide these expressions in the next section.

Constant electrolyte conductivity
One assumption to simplify the electrolyte potential is to take that the electrolyte conductivity is a constant. Then, the integrals involving σ e in (5c), (6) and (8) can be computed analytically.
Hence, the potentials in both electrode and electrolyte can be calculated using Then, the Ohmic losses in the electrolyte can be written as and the electrolyte heat source term can be written as If the electrolyte conductivity is assumed to be high, then an additional simplification can be done by eliminating all the terms involving σ −1 e . Mathematically, this corresponds to the limit Σ e = RT amb F Li 0 σ e,typ 1.

Fast electrode diffusion
Finally, the last simplification that can be taken is fast electrode diffusion, which corresponds to the limit In this case, the particle problem (1) reduces to an ODE of the form and we find several instances of this model in the literature, such as [17,18,19,20].

Results
After presenting the TSPMe model, we now validate it to show its relevance and applicability. To show the accuracy and performance of the reduced model, we start by comparing the TSPMe with the thermal DFN (TDFN). Then, to show the applicability of the TSPMe in real applications we compare it against experimental data on the LG M50 cell. This is a 21700 cylindrical cell, with a nominal capacity of 5 Ah, and with an NMC811 positive electrode and a graphite and SiO x negative electrode (see [28] for further details). In both cases, we perform comparisons at different C-rates (C/2, 1C and 2C) and ambient temperatures (0 • C, 10 • C and 25 • C).

Comparison between TSPMe and DFN
We first look at the comparison between the reduced model (the TSPMe) and a thermal DFN model (using a lumped thermal model, see [16]) under different discharge conditions to assess how well does the TSPMe approximate the TDFN. Here we use a lumped thermal model coupled to the DFN for simplicity, as resolving the battery geometry would require a way heavier computational approach that is out of the scope of this work (see [29] for an example).
The models have been implemented in PyBaMM (Python Battery Mathematical Modelling) package, a Python-based open-source package to simulate battery models [30]. The code and data are available online (see "Data and code availability"). Both models are solved using a finite volumes method for the spatial discretisation using the same number of points in both methods for a given domain. For the results here, we used 30 points in the electrode particles and 20 points for each electrode and separator. This results in a system of 122 ODEs for the TSPMe model and a system of 1262 ODEs and 100 algebraic equations for the DFN model. Observe that, for the same mesh, the DFN model requires solving a system over 10 times larger than the TSPMe, and a lot more complex due to the algebraic constraints. In order to solve the models we used the SciPy ODE solver [31] for the TSPMe, while for the DFN we needed a DAE solver, so we used the CasADI solver [32].
For the simulations we used the parameter values shown in Table 1. The electrochemical parameters are for the LG M50 cell and have been taken from [28], while the thermal parameters (θ, k and h) have been taken from [33]. For the thermal parameters, given that they are lumped parameters for the whole battery, we took from [33] the values for each part and calculated the lumped value using the thicknesses from [28]. For the heat exchange coefficient, we took the intermediate value from the range discussed in [33]. Note that some of the parameters in Table 1 do not appear in the TSPMe model, but their values have been used to calculate the dimensionless parameters in Table A.6, which are needed in the asymptotic analysis presented in Appendix B. These parameters are identified with an asterisk in Table 1.
To capture the microstructure effects on the macroscopic transport, we use the Bruggeman correlation, so we set B k = ε 1.5 k . We also do not include the reversible heat generation effects Q rev because the entropic term Π k = T ∂U k ∂T needs to be consistent with the open-circuit potential U k and this would require a careful and exhaustive experimental analysis that is out of the scope of this work.
The comparison between both models is shown in Figures 5-7, with the error data in Table 2 and the computational times in Table 3. We observe that the TSPMe does a very good job representing the voltage and the cell temperature compared to the TDFN model, and that the discrepancy increases with the C-rate. This could be expected as, from the asymptotic analysis in Appendix B, we know that the some of the assumptions on the parameter sizes break down if the applied current is too large. The data in Table 2 shows a slight trend that the agreement between models improves as the temperature decreases, which could be caused because lower temperatures yield lower values of the thermal potential, which is one of the assumptions in the asymptotic analysis.
Meanwhile, the computational time to solve each model shows that the TSPMe is between approximately 10 times (2C) and 40 times (C/2) faster than the TDFN model. We also notice that, for the TSPMe, the computational time does not vary much with the C-rate, but for the TDFN model the computational time is inversely proportional to the C-rate. This is because we are analysing a full discharge and therefore the discharge time is inversely proportional to the C-rate as well. On the other hand, there is no significant effect of the experiment temperature on the computational time. The Heat exchange coefficient 20 Table 1: Dimensional parameters for the LG M50 cell. The electrochemical parameters are taken from [28] and the thermal parameters are taken from [33]. The parameters marked with an asterisk indicate that they are not directly used in the TSPMe but have been used to evaluate the dimensionless grouping for its derivation (see Appendix B).
solving each model 20 times on a laptop with an Intel Core i7-7660U (2.50Ghz) processor and 16 GB RAM, and the values shown in Table 3 are the mean and standard deviation for different C-rates and temperatures. This huge difference in the computational speed is because the PDEs arising in the TSPMe are discretised into a system of ODEs which is well-behaved, while the discretised TDFN model gives a DAE system which is ill-posed and, therefore, solving this system requires very specific numerical schemes that are computationally expensive [34]. In addition, further reduction methods, such as residue grouping [25] or balanced truncation [26], could be applied to the discretised system of ODEs arising from the TSPMe to reduced the computational time even more.   Table 2: Error between TSPMe and DFN models for different temperatures and C-rates. The first values are the root-meansquared-error (RMSE) and the second values (in brackets) are the peak error.
From this comparison, we conclude that the TSPMe provides a very good approximation to the thermal DFN model at a much lower computational cost. In the next section we assess if this model can be used to obtain meaningful predictions of experimental data.

Comparison with experimental data
Now we compare the TSPMe model with experimental data to assess if the model can be used in real applications. The experimental data is for a commercial cell, the LG M50. For the rate tests, cells were placed inside an Espec thermal chamber and, for accuracy, thermocouples were used to record the temperature inside the chamber. Cell temperature was monitored on the cell mid surface and four cells were used per C-rate and temperature. One of these cells had two additional thermocouples: one on the positive tab and one on the negative tab. A 10 A Digatron battery cycler was used for the tests. For the temperature measurements, K-type thermocouples from RS components were used, which is accurate up to a standard deviation of ±0.75%. Three temperature settings 0 • C, 10 • C and 25 • C were used and for each cycle. The cycle consisted of a constant current (CC) charge step of C/3 to a cut-off voltage of 4.2 V, followed by a constant voltage (CV) step at 4.2 V until the charge current dropped to 0.25 A, then a two-hour long rest period, before a constant current discharge step to 2.5 V and a final two-hour long rest before repeating 0.5C 1C 2C 25 • C 0.53 ± 0.02 0.53 ± 0.02 0.59 ± 0.01 10 • C 0.51 ± 0.03 0.53 ± 0.02 0.60 ± 0.02 0 • C 0.55 ± 0.14 0.52 ± 0.01 0.59 ± 0.02   it but for a different discharge C-rate. Different discharge C-rates were explored: 0.1C, 0.5C, 1C and 2C, running two cycles for each discharge rate. To reduce the data size, a variable time step combined with a voltage difference condition was used. For C/10 and C/2 discharge rates we used 60 s or 10 mV sampling and for 1C and 2C discharge rates we used 1 s or 10 mV sampling. Data was recorded once one of the two conditions was satisfied. Note that one of the cells cycled at 1C (cell 791) has not been included in the validation data because it gave faulty temperature measurements. However, the data for that cell is available in the online repository. We would like to highlight here that the purpose of this work is not to parameterise the thermal model for a specific cell, but to derive and validate a reduced thermal-electrochemical model. Therefore, in the simulations presented in this section we have used existing parameter sets and calibrated and validated them against the experimental data. Still, there are some aspects to discuss about the parameters before we continue. We have used as a starting point the parameters shown in Table 1, but some of them needed to be adjusted. The parameters related to the electrochemistry are taken from [28]. These parameters are assumed to not depend on temperature, except for the reaction rates, which follow an Arrhenius relation. Because the diffusion coefficient in the particles is taken to be constant (and thus it is an effective coefficient, see discussion in [28]), we adjust the diffusion coefficient in the negative electrode for each experiment condition. That is, for each combination of C-rate and temperature, we change the value of D n until we get a reasonable qualitative agreement with experimental data, leaving D p unchanged. In addition, for each temperature we have adjusted the initial concentration in the positive electrode for each temperature, so the initial equilibrium potential matches the rest potential before starting the discharge.
The parameters for the thermal model have been taken from [33] and then tuned to match the data (in particular they were tuned for 1C at 25 • C) obtaining the values of θ = 2.32 · 10 6 J K −1 m −3 and h = 16 W m −2 K −1 . Recall that we have neglected the reversible heat generation term Q rev . However, for the moderate C-rates in the experimental data its contribution should be relatively small. Finally, the value used for both the cell initial temperature and the ambient temperature (which as a modelling assumption we take them to be the same) is taken to be the average value of the final cell temperature of each experiment. The parameters that were adjusted for the comparison with experimental data are shown in Table 4.
The comparison between TSPMe and experimental data for voltage and temperature is shown in Figures 8-10, with the error metrics shown in Table 5. Looking at Figure 8, we observe a very good agreement of the voltage at 25 • C at all C-rates. This could be expected, given that the parameters in [28] were measured at room temperature. On the other hand, we observe that temperature shows very good agreement with data at 2C but not so good at C/2. This could be caused because at C/2 the reversible heat source, which we have neglected in the simulations, might have a relevant impact on the temperature. However, further work is required to validate this hypothesis. In all cases, we notice a very good agreement between TSPMe and temperature data during relaxation, which shows that the tuned value of the heat exchange coefficient   is a good estimation. Looking at the comparison at lower temperatures, we notice that they show worse agreement, especially for voltage. This is reasonable given that, again, the parameters were measured at room temperature. And even the reaction rates, which have a temperature dependence, were measured in the range between 25 • C and 60 • C and therefore the values at 10 • C and 0 • C are an extrapolation. The temperature from TSPMe shows good agreement with data at all experimental temperatures, and provides better agreement at higher C-rates.

Conclusions
In this work, we have derived a reduced thermal-electrochemical model that can be used to simulate the behaviour of lithium-ion batteries in a fast and accurate way (compared to the thermal Doyle-Fuller-Newman model). The reduced model, which we refer to as Thermal Single Particle Model with electrolyte (TSPMe), is a SPM-type model (Single Particle Model) but, as opposed to many of the SPM-type models found in Figure 9: Comparison between TSPMe and experimental data at 10 • C. The colour lines represent the experimental data for the different cells studied (with the cell number matching that of the dataset) and the black line represents the TSPMe. The black cross represents the initial equilibrium potential of the model. the literature, the TSPMe has been derived systematically using asymptotic techniques. This systematic method allows us to determine in advance the range of validity of the model and the accuracy of the reduced model.
Before presenting the TSPMe, we discussed the concept of the Single Particle Model (SPM), both for pure electrochemical models and coupled thermal-electrochemical models (Section 2). From that analysis we concluded that the fundamental feature of SPM-type models is that they split the model into two steps: the first step involves solving a system of partial differential equations (PDEs) to determine the intercalated lithium and electrolyte ion concentrations, while the second step involves finding any other variable of interest from closed form expressions. This is a crucial distinction because the computationally expensive part when solving the model is solving the PDEs. Moreover, when comparing the different SPM-type models in the literature we find that the main differences between them are in the second step, thus the system of PDEs is the same across most of the models. The same idea naturally extends when a thermal-electrochemical model is considered.
The reduction of the thermal DFN model to the TSPMe was based on only two assumptions: the deviations from the open-circuit potential are small and the heat transfer of the battery is limited by the heat exchange with the environment. Based on these assumptions, we could perform the reduction as any other result came as a consequence of the assumptions and the model. In terms of performance, we found that the TSPMe approximates very well the thermal DFN model, especially at low C-rates, but it is significantly faster. At high C-rates (2C) it is over ten times faster while at low C-rates (C/2) the TSPMe is about forty times faster. Moreover, despite all the limitations of the model (due to the assumptions taken to derive it) and the parameters used in the simulations (which were estimated at room temperature), the TSPMe showed good agreement with experimental data on a commercial cell. This good accuracy and low computational time make the TSPMe model (and other models of the same family) a very good candidate for battery design and control.
There are two main lines to explore as possible extensions of this work. On one side, a full thermal parameterisation, including thermal parameters and thermal dependence of electrochemical parameters in the relevant temperature range (0 • C to 25 • C), is required for a better comparison between the TSPMe and experimental data. This will allow us to rule out discrepancies due to the parameters rather than the  Table 5: Error between TSPMe and experimental data for different temperatures and C-rates. The first values are the rootmean-squared-error (RMSE) and the second values (in brackets) are the coefficient of determination R 2 . Note that the latter is dimensionless and the closer it is to one, the better the fit is.
model, and thus obtain a more critical validation of the model. On the other side, the model should also be validated against models resolving the full jellyroll structure of the cylindrical battery (e.g. [29]) to assess the discrepancies between the TSPMe, which assumes a lumped thermal model, and the fully resolved model.
positive particles negative particles negative electrode separator positive electrode Figure A.11: Geometry for the thermal DFN model. The battery is the domain y ∈ Ω batt which can have an arbitrary geometry (here, for example, it is depicted as a cylindrical domain). It is composed of multiple cells which are considered to be one-dimensional in the spatial coordinate x. In turn, porous electrode is modelled as an array of spherical particles, which are described by a one-dimensional spherically symmetric model in the spatial coordinate r.
The electrolyte equations, which are defined in the domain 0 ≤ x ≤ L, are The intercalation reaction between the electrode and the electrolyte is given by The boundary conditions are the following. At the current collector ends we impose The heat equation at the macroscale is given by To simplify the notation, when integrating an electrode variable over the domain 0 ≤ x ≤ L we imply splitting the integral, taking for k = n the domain 0 ≤ x ≤ L n and for k = p the domain L − L p ≤ x ≤ L.
We allow some of the parameters in the model to depend on certain variables. In particular, we take D k to depend on the lithium concentration c k , and D e and σ e to depend on the ion concentration c e . Here, we take t + to be a constant, which is a common modelling assumption. However, the results in this work can be generalised to account for a variable t + . The rest of the parameters are taken to be constant, unless stated otherwise.
We define the following scalings of the problem (A.4) and we choose the time scale t 0 to be the discharge time scale The parameters i 0 and Φ 0 are the typical current and electrode potential, respectively, and the subscript typ denotes the typical value of that parameter.
Then, we can write the dimensionless model as follows. The model for concentration in the particles reads and the potential in each electrode is given by Now, thex domain is defined to be 0 ≤x ≤ n if k = n, and 1 − p ≤x ≤ 1 if k = p. The electrolyte equations, which are defined in the domain 0 ≤x ≤ 1, are The intercalation reaction between the electrode and the electrolyte is given bŷ The boundary conditions at current collector ends arê The dimensionless parameters of the model are From [28,33] we find typical values of the dimensionless parameters, which are shown in Table A

Appendix B. Derivation of the base TSPMe
We now consider the asymptotic limits for the analysis. We take the limits and we assume all the other parameters to be O(1). The reasoning behind these limits is the following. The limits K 1 and Bi 1, as shown in Appendix B.1, give that the temperature is homogeneous in space, which highly simplifies the temperature contributions in the DFN model. This corresponds to the scenario in which the bottleneck in heat transfer is the heat dissipation to the environment. On the other hand, the limit λ 1 is the only assumption needed to break the DFN model into an SPMe as, combined with the implicit assumption that Σ k and Σ e are not small (i.e. O(1) or larger), implies that any deviations from the equilibrium potential are small. Note that at large C-rates this assumption is no longer true because Σ k and Σ e become much smaller than one.
In order to simplify the analysis, we now drop hats from the dimensionless variables as for the rest of the appendix we work with the dimensionless model.

Appendix B.1. Derivation of the reduced thermal model
We start considering the reduction of the thermal model (A.8), as the results of this analysis are helpful in the reduction of the electrochemical model in Appendix B.2. We use the limits λ 1, K 1 and Bi 1. For simplicity, we define δ = λ −1 as our small parameter, so we scale K = δ −1K and Bi = δBi. Then, we can introduce these scalings into (A.8) which, dropping hats, gives We now expand the temperature and the heat source term in powers of δ as We find that, at leading order, the governing equation is so we conclude that T 0 = T 0 (t). We now look at the O(δ) equations to determine T 0 , which are We can now average (B.5a) over the whole domain Ω batt . Applying the divergence theorem and using (B.5b) we find which represent, respectively, the battery averaged heat source term and the surface area per unit volume of the battery. We finally need to write down the leading order source term Q 0 . If we expand (B.2b) we find that there might be an O δ −1 term given by however, as we will see in Appendix B.2, the analysis shows that ∂Φ k0 ∂x = 0 so this term vanishes. Then, we have that the leading order heat source term is After simplifying the electrochemical model in the next section we can provide more detailed expressions for each of the terms that compose Q 0 .

Appendix B.2. Derivation of the reduced electrochemical model
We now reduce the electrochemical model to the base SPMe model just considering the limits λ 1 and using the fact that at leading order the temperature is homogeneous in space, as derived in Appendix B.1. We use again δ = λ −1 1 as our small parameter and rewrite the dimensionless model accordingly. We expand all the variables and derived quantities in powers of δ, using the notation and so on. We expand first (A.6f) which at leading order gives from which we conclude that the potential in the electrodes is spatially homogeneous so Φ k0 = Φ k0 (t). Then, from the leading order term in (A.6m) we conclude Therefore, if the open circuit potential U k is invertible, we have and thus the concentration at the boundary of the particles does not depend on x, i.e. it is the same for all the particles. Then, if the concentration is initially homogeneous in x, it must remain homogeneous at all times so we conclude that c k0 = c k0 (t, r). Using (A.6c) we conclude as well that J k0 = J k0 (t).
The assumption that the concentration is initially homogeneous in x is reasonable, as that should be the case if the battery is left to rest for long enough.
We now take the leading order expansion of (A.6e) with the boundary conditions (A.7). Integrating in each electrode separately we find This is a key result because J k were the terms coupling the four PDEs together. Now that they can be determined a priori the system of PDEs decouples and thus it is a much simpler problem to deal with. Now we can write write the leading order equations for c k and c e . For the electrode particles we have while for the electrolyte we have where, to obtain the boundary conditions (B.15b), we have combined the boundary conditions for molar flux and current in (A.7).
Having the concentration equations that compose the SPMe, we now focus on the potentials. We want to calculate the leading order term of Φ e and the first order term of Φ k so the potentials, in dimensional form, are accurate up to O RT amb F . We focus first on the leading order electrolyte potential Φ e0 . We start calculating i e0 , and expanding (A.6h) jointly with the boundary conditions (A.7) and integrating we obtain Now we can use the leading order expansion of (A.6j) to determine the leading order potential. We have Integrating and using the reference of potential boundary condition in (A.7a) we obtain Finally, we calculate Φ k1 taking the O(δ) equation of (A.6e) combined with the O(1) term in (A.6f), which give where A k and B k are integration constants that need to be determined. These constants are different for each electrode so we determine them separately. Applying the boundary conditions (A.7a) and (A.7b) we find We use (A.6l) and (A.6m) to determine the values of B k . Combining them we have which can be transformed into We could determine B k directly from the expression above, but that would require calculating c k1 . Instead, and given that B k are homogeneous in space, we average (B.23) over each electrode as this approach does not require calculating c k1 . We first show that the averaged c k1 over each electrode is zero. Here we show it for the negative electrode only, but the analysis for the positive electrode is analogous. The concentration c n1 follows the problem where we have used the fact that 1 n n 0 J n1 dx = 0, which can be shown in the same way that we determined J k0 in (B.13). From (B.25) we find thatc n1 ≡ 0 and following a similar argument we can show thatc p1 ≡ 0.
We now average (B.23) over each electrode and isolate B n and B p , which gives across the battery domain Ω batt , the leading order heat source term Q 0 is homogeneous in y and therefore the volume averaging is trivial. Substituting the expressions for the potential found in Appendix B.2 we can rewrite each contribution to the heat source term as In particular, we can reduce the second integral in Q e0 using integration by parts The dimensional version of these equations is presented in Section 3.

Appendix C. Further simplifications
In this section we provide the derivation for the further simplifications to the model presented in Section 3.1. The three main simplifications considered are: quasi-steady-state electrolyte concentration, constant electrolyte conductivity and fast lithium diffusion.
Appendix C.1. Quasi-steady state electrolyte concentration One simplification for the electrolyte concentration is to take the quasi-steady state problem, which is a valid assumption when the current varies over a larger time scale than the electrolyte diffusion and electrolyte ion generation time scales. This corresponds to the asymptotic limit C e = O(δ) and γ e = O(δ). In order to obtain analytical solutions to this problem, we also assume that the diffusion coefficient in the electrolyte D e does not depend on the concentration.
Introducing the scalings C e = δC e and γ e = δγ e into (A.6g), (A.6i) and (A.7) we find that the problem at leading order is given by Even though this problem seems well-posed as it has two boundary conditions, from conservation of charge in the electrolyte we find that they are equivalent. Therefore, we need an extra piece of information, which comes from imposing that the total ion concentration must be conserved. Hence, using the value for the initial concentration, we find 1 0 ε(x)c e0 dx = ε n n + ε s s + ε p p , (C.2) where s = 1 − p − n is the dimensionless separator thickness.
To integrate (C.1) we split the integral into three domain (negative electrode, separator and positive electrode) so B(x) and ε(x) are constants and we use continuity of concentration and flux to match the solutions in each domain, which are given by where ] + − denotes the difference between the limit from the right and the limit from the left at a given point. Then, we find  , if 1 − p ≤ x ≤ 1.
(C.4) The dimensional form is presented in (12). An additional simplification to the previous model is to take γ e = O(1) instead, which gives a small variation to the electrolyte concentration. Then, we can take c e0 = 1.

Appendix C.2. Constant electrolyte conductivity
Another assumption to simplify the electrolyte potential is to take the electrolyte conductivity σ e (c e0 ) = σ e to be a constant. Then, we can obtain analytical expressions for the integrals involving σ e .
The first term to consider is the integral arising in the definition of Φ e0 . We can use the fact that B(x) is piecewise constant and split the integral into integrals over each domain. Then we find 1 Σ e σ e We can now use these expressions to calculate the integrals in Φ n1 and Φ p1 , which give Finally we need to calculate the integral that appears in Q e0 , which is Then, the simplified TSPMe can be written as detailed in Appendix B.3 but with the new definitions of potential and heat source term (C.8a)