Asymptotic reduction, solution, and homogenisation of a thermo-electrochemical model for a lithium-ion battery

TL;DR

This paper develops and compares asymptotic and homogenised models for lithium-ion batteries, demonstrating their accuracy and potential for thermal management without thermal runaway risk.

Contribution

It introduces asymptotic and homogenisation techniques to simplify and analyze complex thermo-electrochemical battery models, enabling analytical solutions and improved thermal management.

Findings

Asymptotic solutions agree with numerical models up to 2C charge rate.

Homogenised models can be solved analytically for practical use.

Thermal runaway does not occur in the model, only potential shifts.

Abstract

We study two thermo-electrochemical models for lithium-ion batteries. The first is based on volume averaging the electrode microstructure whereas the second is based on the pseudo-two-dimensional (P2D) approach which treats the electrode as a collection of spherical particles. A scaling analysis is used to reduce the volume-averaged model and show that the electrochemical reactions are the dominant source of heat. Matched asymptotic expansions are used to compute solutions of the volume-averaged model for the cases of constant applied current, oscillating applied current, and constant cell potential. The asymptotic and numerical solutions of the volume-averaged model are in remarkable agreement with numerical solutions of the thermal P2D model for (dis)charge rates up to 2C, and reasonable agreement is found at 4C. Homogenisation is then used to derive a thermal model for a battery consisting of several connected lithium-ion cells. Despite accounting for the Arrhenius dependence of the reaction coefficients, we show that thermal runaway does not occur in the model. Instead, the cell potential is simply pushed closer to the open-circuit potential. We also show that in many cases, the homogenised battery model can be solved analytically, making it ideal for use in on-board thermal management systems.

Tables

Table 1. Table 1: Parameters applied to the battery

Parameter (Units)	Value (Reference)
$H$ (m)	65$\times {10}^{-3}$ [22]
$x_p$ (m)	70$\times {10}^{-6}$ [22]
$x_n$ (m)	95$\times {10}^{-6}$ [22]
$L$ (m)	129$\times {10}^{-6}$ [22]
${\lambda }_{\text{p}}^{\text{cc}}$ (m)	$1.0\times {10}^{-6}$ [22]
${\lambda }_{\text{n}}^{\text{cc}}$ (m)	$6.2\times {10}^{-6}$ [22]
$A_{\mathrm{cell}}$ (m2)	16.94$\times {10}^{-2}$ [22]
$F$ (C mol-1)	96487 [66]
$R$ (J mol-1 K-1)	8.314 [66]
$T_a$ (K)	298.15 [22]
$I_{\mathrm{app}}$ (A)	2.3 [69]
$i_0$ (A m-2)	13.6 [64]

Equations

$$i_{a,\text{i}}$$

$$\frac{\partial }{\partial x}\left({\varphi }_{a,\text{i}}{i}_{a,\text{i}}\right)$$

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

$$\frac{\partial }{\partial x}\left({\varphi }_{e,\text{i}}{i}_{e,\text{i}}\right)$$

$${\sigma }_e=F^2[{\mu }_L(T_{\text{i}})+{\mu }_A(T_{\text{i}})]c_{L,\text{i}},$$

$${\mu }_L(T_{\text{i}})=\frac{D_L(T_{\text{i}})}{R{T}_{\text{i}}},{\mu }_A(T_{\text{i}})=\frac{D_A(T_{\text{i}})}{R{T}_{\text{i}}},$$

$$\frac{\partial }{\partial x}\left({\varphi }_{a,\text{i}}{i}_{a,\text{i}}+{\varphi }_{e,\text{i}}{i}_{e,\text{i}}\right)=0.$$

$$\frac{\partial }{\partial t}\left({\varphi }_{a,\text{i}}{c}_{a,\text{i}}\right)=\frac{1}{F}\frac{\partial }{\partial x}\left({\varphi }_{a,\text{i}}{i}_{a,\text{i}}\right),$$

$$\frac{\partial c_{a,\text{i}}^{\ast }}{\partial t}=\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^2{D}_{a,\text{i}}\frac{\partial c_{a,\text{i}}^{\ast }}{\partial r}\right),$$

$$-D_{a,\text{i}}\frac{\partial c_{a,\text{i}}^{\ast }}{\partial r}=\frac{1}{F}\left(g_{\text{i}}+C_{\Gamma ,\text{i}}\frac{\partial }{\partial t}\left({\Phi }_{a,\text{i}}-{\Phi }_{e,\text{i}}\right)\right),r=R_{p,\text{i}},$$

$$\frac{\partial c_{a,\text{i}}^{\ast }}{\partial r}=0,r=0.$$

$$\frac{\partial c_{a,\text{i}}^{\ast }}{\partial t}=-\frac{3}{R_{p,\text{i}}F}\left(g_{\text{i}}+C_{\Gamma ,\text{i}}\frac{\partial }{\partial t}\left({\Phi }_{a,\text{i}}-{\Phi }_{e,\text{i}}\right)\right).$$

$$\frac{\partial }{\partial t}({\varphi }_{e,\text{i}}{c}_{L,\text{i}})$$

$$D_e(T_{\text{i}})=D_L(T_{\text{i}})-\theta (T_{\text{i}})\left[D_L(T_{\text{i}})-D_A(T_{\text{i}})\right].$$

$$\theta (T_{\text{i}})=\frac{{\mu }_L(T_{\text{i}})}{{\mu }_L(T_{\text{i}})+{\mu }_A(T_{\text{i}})}=\frac{D_L(T_{\text{i}})}{D_L(T_{\text{i}})+D_A(T_{\text{i}})},$$

$$\frac{\partial }{\partial t}({\varphi }_{e,\text{s}}{c}_{L,\text{s}})$$

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

$$\frac{\partial }{\partial x}\left({\varphi }_{e,\text{s}}{i}_{e,\text{s}}\right)$$

$$g_{\text{i}}=j_{\text{i}}\left(\exp \left[\frac{(1-{\beta }_{\text{i}})F}{R{T}_{\text{i}}}{\eta }_{\text{i}}\right]-\exp \left[\frac{-{\beta }_{\text{i}}F}{R{T}_{\text{i}}}{\eta }_{\text{i}}\right]\right),$$

$$j_{\text{i}}$$

$$K_{\text{eff},\text{i}}(T_{\text{i}})=[K_{a,\text{i}}(T_{\text{i}}){]}^{{\beta }_{\text{i}}}[K_{L,\text{i}}(T_{\text{i}}){]}^{1-{\beta }_{\text{i}}}.$$

$$K_{y,\text{i}}(T_{\text{i}})=K_{y,\text{i}}^0\exp \left[\frac{E_{y,\text{i}}}{R}\left(\frac{1}{T_a}-\frac{1}{T_{\text{i}}}\right)\right],$$

$$K_{\text{eff},\text{i}}=K_{\text{eff},\text{i}}^0\exp \left[\frac{E_{\text{eff},\text{i}}}{R}\left(\frac{1}{T_a}-\frac{1}{T_{\text{i}}}\right)\right],$$

$$U_{\text{i}}$$

$${\rho }_{\text{i}}{\mathsf{c}}_{\text{i}}\frac{\partial T_{\text{i}}}{\partial t}=\frac{\partial }{\partial x}\left(k_{\text{i}}\frac{\partial T_{\text{i}}}{\partial x}\right)+{\varphi }_{e,\text{i}}{q}_{e,\text{i}}+{\varphi }_{a,\text{i}}{q}_{a,\text{i}}+q_{r,\text{i}},$$

$${\rho }_{\text{i}}{\mathsf{c}}_{\text{i}}={\varphi }_{e,\text{i}}{\rho }_{\text{e}}{\mathsf{c}}_{\text{e}}+(1-{\varphi }_{e,\text{i}}){\rho }_{\text{solid},\text{i}}{\mathsf{c}}_{\text{solid},\text{i}},k_{\text{i}}={\varphi }_{e,\text{i}}{k}_e+(1-{\varphi }_{e,\text{i}})k_{\text{solid},\text{i}},$$

$$q_{a,\text{i}}$$

$$q_{e,\text{i}}$$

$$q_{r,\text{i}}$$

$$q_{r,\text{i}}^{\text{irr}}=a_{\text{i}}{g}_{\text{i}}{\eta }_{\text{i}},q_{r,\text{i}}^{\text{rev}}=a_{\text{i}}{g}_{\text{i}}{T}_{\text{i}}\frac{\partial U_{\text{i}}}{\partial T_{\text{i}}},$$

$$q_{r,\text{i}}=a_{\text{i}}{g}_{\text{i}}\left({\Phi }_{a,\text{i}}-{\Phi }_{e,\text{i}}+\frac{E_{L,\text{i}}-E_{a,\text{i}}}{F}\right).$$

$${\rho }_{\text{s}}{\mathsf{c}}_{\text{s}}\frac{\partial T_{\text{s}}}{\partial t}=\frac{\partial }{\partial x}\left(k_{\text{s}}\frac{\partial T_{\text{s}}}{\partial x}\right)+{\varphi }_{e,\text{s}}{q}_{e,\text{s}},$$

$$c_{L,\text{i}}-c_{L,\text{s}}$$

$${\varphi }_{e,\text{i}}\frac{\partial c_{L,\text{i}}}{\partial x}-{\varphi }_{e,\text{s}}\frac{\partial c_{L,\text{s}}}{\partial x}$$

$${\varphi }_{e,\text{i}}{i}_{e,\text{i}}-{\varphi }_{e,\text{s}}{i}_{e,\text{s}}$$

$${\Phi }_{e,\text{i}}-{\Phi }_{e,\text{s}}$$

$${\varphi }_{e,\text{i}}\frac{\partial {\Phi }_{e,\text{i}}}{\partial x}-{\varphi }_{e,\text{s}}\frac{\partial {\Phi }_{e,\text{s}}}{\partial x}$$

Full text

1 Introduction

Lithium-ion batteries (LIBs) are ubiquitous in modern society and are the primary energy source for portable electronic devices and electric cars. Research and development of LIBs has been driven by their high energy density, low self-discharge rate, and lack of memory effects [1]. Despite the attractive features of LIBs, there are major safety concerns due to repeated incidents involving overheating, combusting, or exploding LIB-powered devices [2]. Aside from safety concerns, heat generation also plays an important role in battery ageing [3] and degradation [4, 5], both of which are accelerated at high temperatures. Ensuring safe, reliable, and optimal battery usage requires detailed understanding of the link between the electrochemical processes that drive battery operation and the thermal processes that result in heat generation, distribution, and dissipation [6].

One of the primary causes of rapid heat generation in an LIB is from internal short circuiting, which can result from mechanical abuse, poor battery design, or the formation of lithium dendrites [7]. Large temperature rises can also occur at high discharge rates if a battery is insufficiently cooled, triggering the onset of thermal runaway due to the combination of exothermic chemical reactions and temperature-dependent material parameters. The decomposition of the electrolyte and active material at high temperatures can not only produce heat, but also highly flammable gases [8].

In order to combat the negative effects of heat generation, many battery packs are equipped with an on-board thermal management system (TMS) [9]. The TMS provides real-time monitoring of the temperature and makes adjustments to the battery performance and the cooling system to maintain safe and optimal operation. The TMS is guided by numerical simulations of thermal-electrochemical LIB models, which must be physically accurate yet simple enough to be simulated in real time. The need for rapid simulations is especially important for electric vehicles, as their battery packs consist of thousands of LIBs [10, 11] that need to be monitored individually and as a unit. Simultaneously achieving physical accuracy and low computational complexity remains a major challenge in the design of management systems [12].

Mathematical modelling plays an important role in LIB research because it provides cost-effective insights that can be difficult to obtain through experimental measurement. Modelling has led to improved understanding of many aspects of batteries, including intercalation kinetics [13], active-material utilisation [14] and distribution [15], mechanics [16, 17], phase separation [18, 19, 20, 21, 22, 23], electrode fabrication [24], and model parameter estimation [25, 26]. The development of porous electrode theory by Newman [27, 28] laid the foundation on which nearly every battery model is now built upon. This theory rests on the assumption that the complex microscale geometry of porous electrodes can be simplified through a volume averaging procedure. Doyle *et al. *[29] later argued that diffusion of intercalated lithium across the microscale can be important during battery operation and developed the pseudo-two-dimensional (P2D) model to couple microscale diffusion to macroscale electochemistry. The P2D framework, which treats the electrode as a regular array of spherical particles, is now one of the most common modelling approaches [30, 31].

A wide range of thermo-electrochemical models have been developed [32, 33, 34] in order to investigate short circuiting [35, 36, 37], thermal runaway [38, 39, 40], thermal abuse [41], dominant heat generation mechanisms [22], optimal battery design [42], and heat distribution in battery packs [43, 44]. Many studies are highly computational in nature and make use of large-scale simulations using commercial finite-element software [45], which can be time consuming and offer limited insights. The computational burden of solving a comprehensive thermal-electrochemical model based on the P2D framework can be reduced through cutting-edge numerical implementations or model reduction. As summarised by Jokar *et al. *[30], the electrochemical problem can be simplied by neglecting spatial gradients [46, 47] or the use of polynomial approximations [48, 49, 50], the Galerkin projection method [51, 52, 53], Kalman filtering [54, 55], and Padé approximations [56, 57]. The thermal problem is often simplified by treating the temperature as spatially uniform and integrating the source terms in the heat equation, resulting in a so-called lumped thermal model [55, 58, 59, 60]. Many of these reduction techniques are immediately applied without justification, leading to a lack of robustness for their applicability and range of use. A more systematic approach based on exploitation of the vastly different time scales of various physical processes can lead to a more organic simplification.

Asymptotic techniques provide a systematic means of carrying out model reduction and solution construction with clear ranges of validity, but have received relatively little attention from the battery community. Richardson *et al. *[61] uses homogenisation to upscale a P2D model. Sulzer *et al. *[62] apply asymptotic methods to simplify a model for a lead-acid battery while Marquis *et al. *[63] use similar techniques to obtain simplifications of the P2D model known as the single-particle model. Moyles *et al. *[64] use matched asymptotic expansions to obtain an equivalent circuit model from a volume-averaged model and show that it can reproduce experimental discharge curves for a variety of discharge rates. All of these works are based on isothermal situations, and therefore cannot account for the heat generation mechanisms and battery failures previously discussed. A systematic asymptotic reduction of a thermo-electrochemical model is therefore both novel and timely.

In this paper, we use asymptotic methods to reduce a thermo-electrochemical model and obtain solutions for common modes of battery operation (e.g., fixed current and fixed potential). This reduction is based on a volume-averaged model, similar in nature to the original model proposed by Newman, rather than the P2D model. We first consider the thermo-electrochemical problem within a single cell of an LIB. We account for Arrhenius reaction kinetics in order to explore whether thermal runaway can occur. We then validate our volume-averaged model and the asymptotic solutions against numerical simulations of a thermal P2D model. Due to rapid thermal diffusion, we find that a single-cell model does not predict the large temperature increases that are seen experimentally. Therefore, we apply homogenisation to obtain a model for an LIB that consists of many cells. The increase in thermal diffusion time leads to greater temperature increases, but we find that thermal runaway does not occur, even when Arrhenius reaction kinetics are accounted for. We also find that for realistic parameter values, the homogenised battery model can be solved analytically, making it ideal for use in a TMS.

The outline of the paper is as follows. In Sec. 2 we introduce a model for a single cell of an LIB and non-dimensionalise it. This model is reduced in Sec. 3 by neglecting non-singular terms. Asymptotic solutions are then constructed in Sec. 4 and compared against the P2D model. The cell model is upscaled to a battery model in Sec. 5. Discussions and conclusions follow in Sec. 6.

2 Cell-scale modelling

We consider the thermal-electrochemical processes that occur within the cell of an LIB composed of a positive ($P$) electrode, a separator ($S$), and a negative ($N$) electrode; see Fig. 1. The domain is partitioned so that the positive electrode exists on $0\leq x\leq x_{p}$, the separator on $x_{p}\leq x\leq x_{n}$, and the negative electrode on $x_{n}\leq x\leq L$. The electrode volume can be decomposed into three subdomains corresponding to void space occupied by an electrolyte and solid space occupied by active and inactive materials. The separator is also decomposed into a void and solid fraction but the solid only serves as an electrical insulator that permits ion flow. In each of the phases, conservation of mass, charge, and thermal energy is imposed. A one-dimensional model is employed because of the large aspect ratio between the height and length of a cell, consistent with a number of previous studies [40, 22, 65, 64]. In the equations that follow, the roman subscript $\text{i}=\text{n}$, p, s is used to denote the negative electrode, positive electrode, and separator. Thus, the quantity $\psi_{j}$ in component i of the LIB cell is written as $\psi_{j,\text{{i}}}$. If there is no roman subscript on a variable, then it is assumed that this quantity is uniform across the LIB cell.

2.1 Electrochemical model

The bulk equations of the electrochemical model have been derived from a volume-averaging procedure by Moyles *et al. *[64], where the underlying assumptions and limitations of the model are discussed in detail, and a validation against experimental data is made. In summary, the model is based on the following assumptions:

Dilute solution theory is used to describe the electrolyte

• 2.

Phase separation in the electrodes is not explicitly considered

• 3.

The electrostatic double layer that forms at the matrix-pore interface follows the Helmholtz model.

We now extend this model to account for temperature-dependent material properties and discuss how the microstucture of the electrode can be incorporated into the volume-averaged macroscale equations using the P2D framework proposed by Doyle *et al. *[29]. Temperature induces an autocatalytic effect in the system whereby the reactions generate heat which speed up the reaction rates. Temperature-dependent parameters allow us to explicitly model this effect.

2.1.1 Charge conservation in the electrode solid phase

Ohm’s law relates the current density to the gradient in the electrical potential in the solid phase of the electrode. Following Li *et al. *[22], we account for the capacitance of electronic double layers that form at the surface of the solid matrix of the electrodes. Thus, the conservation of electric charge in the solid phase of the electrodes (i = p, n) is governed by

[TABLE]

where $t$ is time, $T_{\text{{i}}}$ is temperature, $i_{a,\text{{i}}}$ is the current density in the active solid phase, and $\Phi_{a,\text{{i}}}$ and $\Phi_{e,\text{{i}}}$ are the electric potential in the active solid and electrolyte, respectively. Furthermore, $\sigma_{s,\text{{i}}}$ is the electrical conductivity, $C_{\Gamma,\text{{i}}}$ is the capacitance per unit area of the active solid, $\phi_{a,\text{{i}}}$ is the volume fraction of active solid material, and $a_{\text{{i}}}$ is the specific area of active electrode material per unit volume, $a_{\text{{i}}}$ is the surface area of the interface formed between active solid material and the electrolyte. Surface-averaged electrochemical currents produced at the interface between electrolyte and active material are denoted $g_{\text{{i}}}$ and will be specified in Sec. 2.1.6 using Butler–Volmer kinetics.

2.1.2 Charge conservation in the electrode liquid phase

Charge transport in the liquid phase of the electrode (i.e., the electrolyte) occurs via molecular diffusion of lithium ions and anions as well as drift due to the presence of an electric field. Thus, conservation of electronic charge in the liquid phase of the electrodes is described by

[TABLE]

where $F$ is Faraday’s constant; $i_{e,\text{{i}}}$ is the current density; and $c_{L,\text{{i}}}$ and $D_{L}$ are the concentration and diffusivity of lithium ions, respectively; and $D_{A}$ is the diffusivity of the anion. The quantity $\phi_{e,\text{{i}}}$ represents the porosity of the electrode, which is also the volume fraction of electrolyte, and is defined as the ratio of the electrolyte volume to the total volume. The ionic conductivity $\sigma_{e}$ can be written as

[TABLE]

where $\mu_{L}$ and $\mu_{A}$ are the mobility of the lithium ions and anions. Assuming the Nernst–Einstein equation applies, the mobility and the diffusivity are linked through

[TABLE]

where $R$ is the ideal gas constant.

Within the electrode there must be global conservation of charge, which arises from adding (1b) and (2b) to find

[TABLE]

From this point forward, (2b) will be replaced with (5).

2.1.3 Lithium conservation in the electrode solid phase

We consider two approaches for describing lithium conservation in the solid phase of the electrodes. The first is based on volume averaging [64]. In this case, the (volume-averaged) concentration of intercalated lithium $c_{a,\text{{i}}}$ changes due to electrochemical reactions according to

[TABLE]

where diffusion of intercalated lithium across the macroscale has been neglected, which is a valid approximation [64].

The second approach is based on the P2D framework, which treats the electrode as a regular array of spherical particles of radius $R_{p,\text{{i}}}$. During lithiation, lithium enters the surface of each particle and is then transported into the bulk via diffusion. Similarly, during delithiation, the lithium near the surface of the particles is depleted and then replenished by diffusion. Thus, accounting for conservation of lithium in the P2D framework involves solving the radial diffusion equation at each point $x$ in the spatial domain of the electrodes:

[TABLE]

The volume-averaged model (6) and the P2D model (7) are related through the limit of fast microscale diffusion. Integrating (7a) across the particle domain, using the flux conditions (7b) and (7c), and assuming that $c_{a,\text{{i}}}^{*}$ is spatially uniform leads to the simplified equation

[TABLE]

Volume averaging (8) and then using the relation $a_{\text{{i}}}=3\phi_{a,\text{{i}}}/R_{p,\text{{i}}}$ as well as the conservation law (1b) yields (6).

Asymptotic solutions will be derived from the volume-averaged form of the mass conservation given by (6). We then compare the asymptotic solutions to those obtained via numerical simulations of the volume-averaged (VA) and a P2D model based on (7).

2.1.4 Lithium conservation in the electrode liquid phase

Lithium conservation in the liquid phase of the electrode leads to a reaction-diffusion equation:

[TABLE]

In this equation, the effective ionic diffusivity $D_{e}$ in the electrolyte is given by

[TABLE]

The parameter $\theta$ plays the role of a transference number and is given by

[TABLE]

where the second equality arises from using the Nernst–Einstein equations (4).

2.1.5 Charge transport and mass conservation in the separator

Conservation of charge and lithium in the separator is given by

[TABLE]

2.1.6 Electrochemical kinetics

The electrochemical reactions that are responsible for lithiation and delithiation are described using the Bulter–Volmer kinetics [66] in the Helmholtz limit of thin electric double layers. The electrochemical current produced at the electrode-electrolyte interface is given by

[TABLE]

The reaction constants are assumed to have an Arrhenius form

[TABLE]

where $E_{y,\text{{i}}}$ is the activation energy of the reaction. By inserting (14) into (13c), the effective reaction constant can also be written in Arrhenius form,

[TABLE]

where $K_{\text{eff},\text{{i}}}^{0}=(K_{y,\text{{i}}}^{0})^{\beta_{\text{{i}}}}(K_{L,\text{{i}}}^{0})^{1-\beta_{\text{{i}}}}$ and $E_{\text{eff},\text{{i}}}=\beta_{\text{{i}}}E_{a,\text{{i}}}+(1-\beta_{\text{{i}}})E_{L,\text{{i}}}$. Aside from the terms involving $c_{a,\text{{i}}}^{\text{max}}$, the electrochemical kinetics described by (13a)–(14) are consistent with those derived by Ferguson and Bazant [19] using non-equilibrium thermodynamics. The theoretical open-circuit potential, which is consistent with the Butler–Volmer kinetics (13), is given by

[TABLE]

Extended forms of the open-circuit potential can be used to account for additional physics such as phase separation, relevant for LFP electrodes, and multiple lithiation stages, relevant for graphite electrodes [19, 20, 21, 67].

2.2 Thermal model

Conservation of energy in each electrode takes the form

[TABLE]

where

[TABLE]

are the phase-averaged volumetric heat capacity and thermal conductivity, respectively. We note that (17) includes the inactive material which can still transport heat. Furthermore, the equation is derived by assuming that (i) there is local thermal equilibrium in the three phases so that a common temperature is well defined and (ii) the inactive and active solid material have the same thermal properties. The source terms in (17) are given by

[TABLE]

representing Ohmic heating in the active solid material, ionic Ohmic heating in the electrolyte, and heat generation due to surface electrochemical reactions, respectively. Full derivations of (17) can be found in [68, 48] while discussion of the electrochemical surface reaction heat can be found in [22, 40, 66]. Often, the electrochemical heat generation (19c) is decomposed into irreversible and reversible contributions given by

[TABLE]

respectively. The reversible contribution $q_{r,\text{{i}}}^{\text{rev}}$ can cool the battery whereas the irreversible contribution $q_{r,\text{{i}}}^{\text{irr}}$ always heats the battery. By inserting the expression for the open-circuit potential (16) into the equation for the electrochemical heat generation (19c), we find that

[TABLE]

The electrochemically inert separator can also transport heat where conservation of energy is given by

[TABLE]

where the phase-averaged volumetric heat capacity and thermal conductivity are given by (18) and $q_{e,\text{{s}}}$ is given by (19b).

2.3 Boundary and initial conditions

The electrolyte is free to flow between the voids of the electrodes and separator. Therefore, we impose continuity of concentration, molar flux of lithium ions, electrolyte current density, electrolyte potential at the separator-electrode interfaces $x=x_{p}$ and $x=x_{n}$. Furthermore, we also require the temperature and heat flux to be continuous. These stipulations result in the following boundary conditions

[TABLE]

The process of volume averaging microscopic boundary conditions for flux-like quantities introduces volume fractions in (23b), (23c), and (23e) accounting for porosity difference in the materials. The volume fractions are absent in (23g) because the temperature has been phase averaged under the conditions of local thermal equilibrium.

No current can pass through the solid material of the separator because it does not conduct electricity. Therefore, we impose

[TABLE]

The electrode surfaces at $x=0$ and $x=L$ are in contact with current collectors which allow for the flow of electrons into and out of the battery. Ionic flow through the current collectors is not permitted and therefore at the electrode-collector interface, the ionic molar fluxes and electrolyte current must vanish, leading to

[TABLE]

We require a grounding condition and without loss of generality, set the electrolyte potential at the electrode-collector interface of the negative electrode to zero, yielding

[TABLE]

There are two common modes of operation: constant current and constant voltage. The constant current operation involves drawing a prescribed current from the cell and measuring the change in electric potential across the cell as a function on time. To model this scenario, we impose a draw current at the positive electrode given by

[TABLE]

where $\mathcal{I}$ is the C-rate of the battery and $i_{0}$ the draw current density at a discharge rate of 1C (corresponding to $\mathcal{I}=1$). The C-rate is a commonly used measure of how quickly a battery is being charged or discharged. The negative on the right-hand side of (27) means that positive C-rates correspond to a discharge process. Alternatively, in the constant voltage operation, the cell voltage $\Delta V$ defined by

[TABLE]

is prescribed instead. When (27) is prescribed then the cell voltage becomes an output of the system. When (28) is used, then (27) is instead used to determine the C-rate $\mathcal{I}(t)$.

We assume that the current collectors are in contact with the exterior environment and admit heat transfer with an ambient temperature, $T_{a}$. We model this using Newton cooling conditions

[TABLE]

where denote heat transfer coefficients.

The initial concentrations are assumed to be spatially uniform throughout each of the battery components and are given by $c_{a,\text{{i}}}(x,0)=c_{a,\text{{i}}}^{0}$ and $c_{L,\text{{i}}}(x,0)=c_{L}^{0}$. The electric potentials satisfy $\Phi_{a,\text{{i}}}(x,0)-\Phi_{e,\text{{i}}}(x,0)=U_{\text{{i}}}$, that is, the overpotential is initially zero. The initial temperature is $T_{\text{{i}}}=T_{a}$.

2.4 Parameter determination

The parameter values are adapted from Refs. [65, 22], in which thermal P2D models are developed for ANR26650m1-a lithium-ion cells [69]. Tables of parameters values along with descriptions of how they were obtained are given in A.

2.5 Numerical implementation

The VA and P2D models are solved numerically using a semi-implicit finite-difference method. At each time step, a fixed-point method is used to solve the full nonlinear system, which is decomposed into simpler subproblems that are solved sequentially. Each fixed-point iteration begins by solving for the solid-phase potentials and current densities using Newton’s method. The liquid-phase variables are then obtained by first calculating the current density, followed by the lithium concentration, and then the electrical potential. The last step of the fixed-point iteration is to update the concentration of intercalated lithium using either the VA or P2D approach. The fixed-point iterations repeat until a stopping condition is reached. For a prescribed applied current, the stopping condition is met when the change in cell potential is less than 1 mV. In the case of a prescribed cell potential, the stopping condition is met when then change in the C-rate is less than $10^{-4}$. Upon satisfying the stopping condition, the electrical problem is assumed to be solved, the temperature is updated, and then the process repeats.

2.6 Non-dimensionalisation

We introduce characteristic scales for each of the variables to write the model in dimensionless terms. The motivation behind the choice of scales for the electrochemical model is discussed in detail in Moyles *et al. *[64]. The electric potentials are written as $\Phi_{a,\text{{i}}}=(RT_{a}/F)[\log\mathcal{U}_{\text{{i}}}+\Phi_{a,\text{{i}}}^{\prime}]$, $\Phi_{e,\text{{i}}}=(RT_{a}/F)\Phi_{e,\text{{i}}}^{\prime}$, and $U_{\text{{i}}}=(RT_{a}/F)[\log\mathcal{U}_{\text{{i}}}+U_{\text{{i}}}^{\prime}]$, where $\log\mathcal{U}_{\text{{i}}}$ is the initial value of the open-circuit potential normalised by the thermal voltage. Similarly, the temperature is written as $T_{\text{{i}}}=T_{a}+(\Delta T)T_{\text{{i}}}^{\prime}$. The temperature scale $\Delta T$ will be discussed below. Each temperature-dependent parameter $\psi(T_{\text{{i}}})$ is rescaled according to $\psi(T_{\text{{i}}})=\psi^{0}\psi^{\prime}(T_{\text{{i}}}^{\prime})$ where $\psi^{0}=\psi(T_{a})$. The only exception to this is the transference number $\theta$, which is written as $\theta(T_{\text{{i}}})=\theta^{\prime}(T_{\text{{i}}}^{\prime})$. Spatial variables are scaled according to $x=Lx^{\prime}$, $x_{p}=Lx_{p}^{\prime}$, and $x_{n}=Lx_{n}^{\prime}$. Time is written as $t=(L^{2}/D_{e}^{0})t^{\prime}$. The current densities are scaled as $i_{a,\text{{i}}}=i_{0}i_{a,\text{{i}}}^{\prime}$ and $i_{e,\text{{i}}}=i_{0}i_{e,\text{{i}}}^{\prime}$. The solid-phase and liquid-phase lithium concentrations are written as $c_{a,\text{{i}}}=c_{a,\text{{i}}}^{0}+(\Delta c)c_{a,\text{{i}}}^{\prime}$ and $c_{L,\text{{i}}}=c_{L,\text{{i}}}^{0}+(\Delta c)c_{L,\text{{i}}}^{\prime}$, where $\Delta c=(i_{0}L)/(FD_{e}^{0})$. The overpotential is written in terms of the thermal voltage, $\eta_{\text{{i}}}=(RT_{a}/F)\eta_{\text{{i}}}^{\prime}$. The electrochemical surface current is scaled as $g_{\text{{i}}}=g_{\text{{i}}}^{0}g_{\text{{i}}}^{\prime}$ where

[TABLE]

We select a common, constant porosity $\phi_{e,\text{{n}}}=\phi_{e,\text{{p}}}=\phi_{e,\text{{s}}}=\phi_{\text{{e}}}$ for each of the cell components. We also assume that the volume fractions of active solid are constant in time and uniform in space.

2.6.1 Electrochemical model

Upon dropping the primes, the following bulk equations for the electrode active material are obtained:

[TABLE]

The governing equations for the electrolyte in the electrode are given by

[TABLE]

where the non-dimensional effective diffusivity and ionic conductivity are given by

[TABLE]

The bulk equations in the separator are

[TABLE]

The non-dimensional Butler–Volmer kinetics can be written as

[TABLE]

with $\Gamma_{\text{eff},\text{{i}}}=\gamma_{T}E_{\text{eff},\text{{i}}}/(RT_{a})$. The definitions of $\Gamma_{\text{eff},\text{{i}}}$ and $K_{\text{eff},\text{{i}}}$ lead to a simplified form in the limit $\gamma_{T}T_{\text{{i}}}\ll 1$ as shown in (50c). The non-dimensional overpotential can be written as

[TABLE]

where the constant ($\mathcal{U}_{\text{{i}}}$) and composition-dependent ($\mathcal{V}_{\text{{i}}}$) contributions to the open-circuit potential are defined according to

[TABLE]

The non-dimensional numbers are described in [64] and given by

[TABLE]

where $\sigma_{e}^{0}=F^{2}(\mu_{L}^{0}+\mu_{A}^{0})c_{L}^{0}$ and $D_{e}^{0}=D_{L}^{0}-\theta^{0}(D_{L}^{0}-D_{A}^{0})$.

2.6.2 Thermal model

We choose a temperature scale based on the heat generated from the electrochemical reactions in the positive electrode,

[TABLE]

anticipating this to be the dominant heat generation mechanism. Indeed, using parameters from Tables 3 and 3, we estimate the heat generation from Ohmic heating to be $q_{a,\text{{p}}}\sim 86$ W m*-3*, $q_{a,\text{{n}}}\sim 1.8$ W m*-3*, and $q_{e,\text{{i}}}\simeq 73$ W m*-3* while for the electrochemical heat, we find that $q_{r,\text{{p}}}\sim 1.4\times 10^{4}$ W m*-3* and $q_{r,\text{{n}}}\sim 1.3\times 10^{4}$ W m*-3*. Evaluating (40) gives a temperature scale of $\Delta T\sim 0.65$ K.

The resulting dimensionless heat equations become (upon dropping primes):

[TABLE]

with additional non-dimensional numbers

[TABLE]

where $\mathrm{Le}$ is the Lewis number which represents the ratio of thermal diffusion to mass diffusion and $\mathcal{Q}_{i}$ is the ratio of Ohmic heat generated to that from surface reactions. Furthermore, $\varrho_{\text{{i}}}$, $\mathcal{K}_{\text{{i}}}$, and $\alpha_{\text{{n}}}$ are the ratios of electrode and ionic conductivities, phase-averaged volumetric heat capacities, phase-averaged thermal conductivities, and surface electrode currents respectively. The non-dimensional heat sources $q_{\text{{i}}}$ are given by

[TABLE]

2.6.3 Boundary and initial conditions

The dimensionless boundary conditions at the positive electrode-collector interface are

[TABLE]

where $\mathcal{H}_{\text{{n}}}=h_{\text{{n}}}/h_{\text{{p}}}$ is the ratio of negative and positive electrode heat transfer coefficients. Under galvanostatic or galvanodynamic conditions, we impose

[TABLE]

Alternatively, the dimensionless cell potential $\Delta v=(F\Delta V)/(RT_{a})-\log\mathcal{U}_{\text{{p}}}+\log\mathcal{U}_{\text{{n}}}$ can be prescribed and $\mathcal{I}(t)$ solved for.

Finally, the non-dimensional initial conditions are given by $c_{a,\text{{i}}}(x,0)=c_{L,\text{{i}}}(x,0)=\Phi_{a,\text{{i}}}(x,0)-\Phi_{e,\text{{i}}}(x,0)=T_{\text{{i}}}=0$.

3 Reduction of the cell model

Using the parameters in Tables 3–3 results in the following non-dimensional numbers associated with the electrochemical model:

[TABLE]

The corresponding non-dimensional parameters for the thermal model are

[TABLE]

The cell-scale model will now be simplified by exploiting the smallness of these parameters.

3.1 Simplification of the thermal problem

The heat equation (41) can be simplified by neglecting terms relating to Ohmic heating ($\mathcal{Q}_{a,\text{{i}}}\ll 1$ and $\mathcal{Q}_{\text{{e}}}\ll 1$). The large Lewis number, $\mathrm{Le}$, implies that rapid thermal diffusion will lead to a temperature profile that is roughly uniform across the cell. Thus, the temperatures are written as $T_{\text{{i}}}(x,t)=T(t)+\mathrm{Le}^{-1}\tilde{T_{\text{{i}}}}(x,t)$. Integrating the heat equations across their respective domains, adding them together, and using the boundary conditions results in a single heat equation for the dimensionless cell temperature $T$ given by

[TABLE]

where $\langle\varrho\rangle_{\text{cell}}=x_{p}+(x_{n}-x_{p})\varrho_{\text{{s}}}+(1-x_{n})\varrho_{\text{{n}}}$ is the cell-averaged dimensionless heat capacity and

[TABLE]

are the component-averaged volumetric heating terms due to electrochemical reactions. In deriving (48), we have taken $\mathrm{Le}\,\mathrm{Bi}=O(1)$ in size, as justified by the parameter values in (47), which is a distinguished limit of the model. For $\mathrm{Le}\,\mathrm{Bi}\gg 1$ there is rapid heat exchange with the environment, leading to small temperature changes of size $O(\mathrm{Le}^{-1}\mathrm{Bi}^{-1})$. For $\mathrm{Le}\,\mathrm{Bi}\ll 1$, heat exchange occurs slowly and the temperature can undergo substantial variations. The small terms of size $O(\mathrm{Le}^{-1})$ in the heat equation (48) will be neglected in the analysis below.

3.2 Simplification of the electrochemical kinetics

The dimensionless Bulter–Volmer kinetics (36) can be simplified using the fact that $\gamma_{T}\ll 1$. We also set $\beta_{\text{{i}}}=1/2$ in order to write asymptotic solutions in terms of hyperbolic trigonometric functions. This is equivalent to the common assumption of symmetric reactions at the electrode surface [22, 70, 35, 18]. Thus, the simplified electrochemical kinetics are given by

[TABLE]

To simplify the overpotential given by (37), we first note that the combination $\log\mathcal{U}_{\text{{i}}}+\Delta\Gamma_{\text{{i}}}$ is $O(1)$ in magnitude despite the individual contributions being large for the positive electrode. The term $\gamma_{T}T(\log\mathcal{U}_{\text{{i}}}+\Delta\Gamma_{\text{{i}}})$ will therefore be small and can be neglected, leading to

[TABLE]

with $\mathcal{V}_{\text{{i}}}$ as in (38). The simplified surface heating terms are given by

[TABLE]

which is exactly (43c), but with $g_{\text{{i}}}$ defined by (50).

Equation (50c) implies that Arrhenius effects will become relevant if the dimensionless temperature reaches values given by $T\sim\Gamma_{\text{eff},\text{{n}}}^{-1}\sim 38$, which is equivalent to a dimensional temperature rise of $\Delta T\sim 25$ K. As shown in Sec. 5, a temperature rise of this magnitude is achievable by coupling several LIB cells together.

3.3 Simplification of the electrochemical model

Due to the spatial uniformity of the temperature, the problem for the concentration of lithium ions in the electrolyte simplifies to

[TABLE]

together with the boundary conditions (44b), (44d), (44e), and (44n).

Reductions also apply to the electrochemical problem for the solid phase (32) but these depend on the type of discharge condition (fixed current or fixed voltage) used and will be discussed appropriately.

4 Analysis of the cell problem

The reduced model is now used to understand the thermal behaviour of an LIB cell under various operation modes. Matched asymptotic expansions are used to calculate approximate solutions in key time regimes. The parameters $\nu_{a,\text{{i}}}$ and $\gamma_{c}$ will be assumed to have the same order of magnitude as $\nu_{\text{{e}}}$. By doing so, a regular asymptotic expansion in powers $\nu_{\text{{e}}}$ can be used in each time regime to systematically capture the effects of gradients in the lithium concentration and electrical potentials.

4.1 Galvanostatic discharge

We first consider the case where the discharge current $\mathcal{I}$ is constant in time. The isothermal version of this problem has been solved using asymptotic methods by Moyles *et al. *[64], who found four distinct time regimes, the most relevant of which for this mode of battery operation occurs when $t=O(\gamma_{c}^{-1})$. On this time scale, large changes in the concentration of intercalated lithium occur, whereas the liquid-phase lithium concentration remains $O(1)$ in size due to a balance between diffusion and electrochemical consumption and generation. Following Ref. [64], we rescale the variables as $t=\gamma_{c}^{-1}\hat{t}$, $c_{a,\text{{i}}}=\gamma_{c}^{-1}\hat{c}_{a,\text{{i}}}$. It is also convenient to introduce a rescaled temperature $\hat{u}$ defined by $T=\mathcal{T}\hat{u}$, where $\mathcal{T}=\mathrm{Le}^{-1}\mathrm{Bi}^{-1}(1+\mathcal{H}_{\text{{n}}})^{-1}\mathcal{G}_{\text{{p}}}^{-1}$. All of the variables on this time scale are denoted with a hat and written as asymptotic expansions of the form $\hat{f}=\hat{f}^{(0)}+\nu_{\text{{e}}}\hat{f}^{(1)}+O(\nu_{\text{{e}}}^{2})$. Capacitance effects are negligible on this time scale and do not need to be considered in (32c).

4.1.1 The leading-order problem

The leading-order problem is obtained by setting $\nu_{\text{{e}}}=0$, $\nu_{a,\text{{i}}}=0$, and $\gamma_{c}=0$ in the rescaled system. Equations (32b), (33b), and (35b) imply that the electric potentials are uniform in space to leading order. The grounding condition on the electrolyte potential gives that $\hat{\Phi}_{e,\text{{i}}}^{(0)}(x,t)=0$. Since the concentration of intercalcated lithium is initially uniform in space, the spatially independent potentials further imply that the leading-order contribution $\hat{c}_{a,\text{{i}}}^{(0)}$ remains uniform for all time. Thus, integration of (32a) over its respective domain shows that the solid-phase concentration changes linearly in time according to

[TABLE]

Similarly, the leading-order parts of (32c) can be integrated in space to yield nonlinear equations for the solid-phase potentials,

[TABLE]

where the electrochemical kinetics are given by

[TABLE]

with $\hat{\eta}_{\text{{i}}}^{(0)}=\hat{\Phi}_{a,\text{{i}}}^{(0)}-\log\hat{\mathcal{V}_{\text{{i}}}}^{(0)}$ and

[TABLE]

The current densities in the solid-phase of the electrodes are given by

[TABLE]

and the current densities in the electrolyte are

[TABLE]

The leading-order liquid-phase lithium concentration is given by its quasi-steady profile, $\hat{c}_{L,\text{{i}}}^{(0)}=\mathcal{I}[1-\theta(\mathcal{T}\hat{u}^{(0)})]/[D_{e}(\mathcal{T}\hat{u}^{(0)})\phi_{\text{{e}}}]c_{L,\text{{i}}}^{\infty}$, where

[TABLE]

and $\mathcal{B}=(1/3)[(1-x_{n})^{2}-x_{p}^{2}]-1$.

We are now in a position to solve for the leading-order contribution to the temperature. The electrochemical heating terms (52) can be simplified using (55). Upon taking the limit $\gamma_{c}\to 0$ in the rescaled heat equation, we find that the leading-order temperature evolves in a quasi-static manner so that

[TABLE]

Due to the temperature dependence of the reaction constant in (56c), the expression in (61) must be regarded as a nonlinear equation for the temperature that must be solved at each point in time. If the temperature dependence is neglected, then the leading-order electrochemical problem completely decouples from the thermal problem. In this case, there is no positive feedback mechanism that could lead to thermal runaway and unsafe overheating of the LIB cell. Even if there was a strong dependence of $K_{\text{eff},\text{{i}}}$ on temperature, then from the electrochemical kinetics (56) it can be seen that a large increase in $K_{\text{eff},\text{{i}}}$ only acts to push the electric potentials closer to the open-circuit potentials ($\hat{\Phi}_{a,\text{{i}}}^{(0)}\to\log\mathcal{V}_{\text{{i}}}^{(0)}$). Consequently, there is no mechanism for auto-catalytic heat generation and thermal runaway. However, it is important to point out that the terms associated with electrochemical heat generation become singular as the battery becomes fully charged or fully drained. This singularity is due to an incompatibility with the constant-current boundary condition and results in the solid-phase potentials blowing up at the end of the (dis)charge process [64].

4.1.2 The first-order problem

The corresponding equations for the electrode solid phase are

[TABLE]

which are supplemented with the boundary and initial conditions $\hat{\imath}_{a,\text{{p}}}^{(1)}(0,\hat{t})=\hat{\imath}_{a,\text{{p}}}^{(1)}(x_{p},\hat{t})=0$, $\hat{\imath}_{a,\text{{n}}}^{(1)}(x_{n},\hat{t})=\hat{\imath}_{a,\text{{n}}}^{(1)}(1,\hat{t})=0$, and $\hat{c}_{a,\text{{i}}}^{(1)}(x,0)=0$. The electrolyte potential can be obtained from

[TABLE]

with $\hat{\Phi}_{e,\text{{n}}}^{(1)}(1,t)=0$ and continuity at $x=x_{p}$ and $x=x_{n}$. The correction to the Butler–Volmer reaction current is given by

[TABLE]

where $\eta_{\text{{i}}}^{(1)}=\Phi_{a,\text{{i}}}^{(1)}-\Phi_{e,\text{{i}}}^{(1)}-\mathcal{V}_{\text{{i}}}^{(1)}/\mathcal{V}_{\text{{i}}}^{(0)}$ and

[TABLE]

In order to solve this problem, we first note that integration of (62a) over the electrode domain and using the boundary and initial conditions reveals that the net corrections to the solid-phase concentrations are zero:

[TABLE]

The correction to the steady-state electric potential of the electrolyte can be written as

[TABLE]

The corrections to the solid-phase potential are given by $\hat{\Phi}_{a,\text{{p}}}^{(1)}=\Phi_{a,\text{{p}}}^{\infty}+P_{p}$ and $\hat{\Phi}_{a,\text{{n}}}^{(1)}=\Phi_{a,\text{{n}}}^{\infty}+P_{n}$, where

[TABLE]

and $P_{i}$ are the constants of integration. The constants $P_{i}$ are determined by solvability conditions that are obtained through the integration of (62c), yielding

[TABLE]

Due to the conditions (66) and the linear dependence of $c_{a,\text{{i}}}^{(1)}$ in $g_{\text{{i}}}^{(1)}$, the contributions from $c_{a,\text{{i}}}^{(1)}$ to the integrals in (69) become zero, implying that the correction to the concentration of intercalated lithium does not influence the correction to the electrical problem. Having determined the constants $P_{i}$, the correction to the (dimensionless) cell potential is simply given by $\Delta v^{(1)}=P_{p}-P_{n}$.

The correction to the temperature simplifies due to the contributions from $g_{\text{{i}}}^{(1)}$ integrating to zero due to (69). Thus, we find that

[TABLE]

4.1.3 Comparison with numerics

The asymptotic solutions are compared with numerical simulations of the VA and P2D models by considering discharge processes of 1C, 2C, and 4C. These correspond to $\mathcal{I}=1$, $2$, and $4$, respectively. The resulting discharge curves (cell potential as a function of time) and cell temperatures are shown in Fig. 2. From the parameters in (47), $\mathcal{T}\simeq 3.8\times 10^{-2}$ and therefore we neglect the temperature dependence of parameters other than the reaction rate. The asymptotic solutions that are constructed from only the leading-order contributions (see panels (a)–(b)) are in good agreement with the numerical solution of the VA model across all of the C-rates. When the first-order contributions are included in the asymptotic solutions (see panels (c)–(d)), the agreement is nearly perfect. There is also good agreement between the cell potentials and temperatures computed from VA and P2D models, even at a high discharge rate of 4C. The VA-P2D model divergence occurs near the end of the discharge process, which induces a singularity in the open-circuit potential. This divergence can be explained in terms of the different ways in which the models treat the solid concentrations. The VA model uses a bulk concentration which is equivalent to a uniform particle concentration. The P2D model however captures particle-scale diffusion which induces a concentration gradient between the surface and center of the particle. For slow diffusion in the particles, the surface concentration and hence the particle concentration gradients are amplified. Near a singularity in the open-circuit potential, a change in the surface concentration due to slow diffusion can have a large impact on the electrochemical kinetics and drive the P2D model away from the VA model. A divergence between the models at the start of the discharge process, which is also close to singular point, is not observed in Fig. 2 because the time scale over which the concentration rapidly leaves the singular zone is fast. The electrode-averaged surface concentrations $c_{a,\text{{i}}}^{\text{surf}}$ are plotted as functions of time in Fig. 2 (e)–(f). In the VA model, the surface and bulk solid-phase concentrations are the same. Furthermore, when time is normalised by the discharge time $t_{\text{dis}}=(3600/\mathcal{I})$ s, the results lie on a master curve (shown as the black dashed line) that is independent of the C-rate. The positive (negative) electrode (de)lithiates on discharge. The surface concentrations in the P2D model are therefore greater (smaller) compared to the bulk, which adds resistance and causes a greater voltage drop relative to the VA model.

As predicted by the asymptotic analysis, the temperature change that occurs during discharge is small. This is due to the smallness of the Biot number, which implies there is excellent heat exchange with the environment. The strong agreement between asymptotic and numerical solutions also confirms the dominant mechanism of heat generation is due to electrochemical reactions rather than Ohmic heating.

4.2 Potentiostatic hold

We now consider the behaviour of the battery during a potentiostatic hold, whereby the dimensionless cell potential is held at a constant value of $\Delta v$. There are four time regimes that occur which are similar to case of the galvanostatic discharge. The first two, defined by $t=O(\mathcal{C}_{\text{{n}}}\mathcal{G}_{\text{{n}}}\nu_{a,\text{{n}}}/\phi_{a,\text{{n}}})$ and $t=O(\mathcal{C}_{\text{{p}}})$, are extremely short and capture the formation of double layers and the onset of electrochemical reactions. The third, defined by $t=O(1)$, involves the dynamics on the diffusive time scale and captures the transient behaviour of the lithium concentration and voltage profiles in the electrolyte. The fourth regime occurs on the saturation time scale, $t=O(\gamma_{c}^{-1})$, and describes the approach to a quiescent state.

The majority of heat generation occurs in the third and fourth time regimes, which will be the focus of our asymptotic analysis. The governing equations in these regimes are solved using matched asymptotic expansions in time and regular expansions in powers of $\nu_{\text{{e}}}$. The two capacitance regimes are briefly discussed and further details can be found in B. The temperature dependence of the parameters can be ignored due to the small amount of heat that is generated over these short time regimes.

4.2.1 Dynamics on the capacitance time scales

The first capacitance time regime occurs when $t=O(\mathcal{C}_{\text{{n}}}\mathcal{G}_{\text{{n}}}\nu_{a,\text{{n}}}/\phi_{a,\text{{n}}})$. A small-time analysis in this regime reveals that instantaneously changing the cell potential results in a large initial C-rate given by

[TABLE]

where $\langle\nu\rangle_{\text{cell}}$ is the cell-averaged resistivity, which is based on the equivalent resistivities of the cell components

[TABLE]

The battery resistivity, $\langle\nu\rangle_{\text{cell}}$, shows that the electrodes and separator behave as if connected in series, whereas the phase resistivities, $\nu_{\mathrm{eq},\text{{i}}}$, show that the porous phases behave as if connected in parallel. The initial current is thus set by the most resistive contribution to $\langle\nu\rangle_{\text{cell}}$, which in this case comes from the positive electrode. For larger times, the electrical resistance in the negative electrode significantly increases due to the onset of electrochemical reactions. Consequently, the magnitude of the C-rate undergoes a marked decrease and plateaus at an approximate value given by the solution of the nonlinear equation

[TABLE]

The second capacitance time regime is defined by $t=\mathcal{O}\left(\mathcal{C}_{\text{{p}}}\right)$. This regime marks the onset of electrochemical reactions in the positive electrode which leads to further increases in the resistance and decreases in the C-rate. Eventually, the C-rate reaches another plateau and settles to a value given by the solution of

[TABLE]

Despite the large currents that can arise in these two regimes, the change in lithium concentration and temperature remains negligible. Therefore, we impose the original initial conditions in the third (diffusive) time regime.

4.2.2 Dynamics on the diffusive time scale

The concentration-dependent contribution to the open-circuit potential and the exchange current density are expanded in powers of $\nu_{\text{{e}}}$ and the leading- and first-order terms are given by $\mathcal{V}_{\text{{i}}}^{(0)}=1$, $j_{\text{{i}}}^{(0)}=1$,

[TABLE]

In the leading-order problem, the potentials and concentration of intercalated lithium are spatially uniform, implying that $\Phi_{e,\text{{i}}}^{(0)}=0$ and $\Phi_{a,\text{{p}}}^{(0)}=\Phi_{a,\text{{n}}}^{(0)}+\Delta v$. Following a similar strategy as in Sec. 4.1.1, we find that the leading-order electrical problem is given by

[TABLE]

The solution for $\mathcal{I}^{(0)}$ matches with the leading-order solution (in terms of $\nu_{\text{{e}}}$) to (74). The leading-order solid-phase potentials and C-rate are constant in time, implying that the solid-phase concentrations evolve linearly in time according to (54). Furthermore, the leading-order current densities are given by (58)–(59). Having solved the electrical problem, the leading-order thermal problem for the scaled temperature $u^{(0)}=\mathcal{T}^{-1}T^{(0)}$ is given by

[TABLE]

with $u^{(0)}(0)=0$, which has a simple analytical solution that tends to a constant value.

The first-order problem requires the leading-order solution for the liquid-phase lithium concentration, which satisfies (53). Since the electrolyte current density $i_{e,\text{{i}}}^{(0)}$ is constant in time, the solution is straightforwardly obtained using separation of variables. By writing the concentration in terms of the deviation from the quasi-steady profile, $c_{L,\text{{i}}}^{(0)}=[\mathcal{I}^{(0)}(1-\theta^{0})/\phi_{\text{{e}}}](c_{L,\text{{i}}}^{\infty}+C_{L})$, it is possible to consider a single liquid-phase concentration $C_{L}$ for the entire cell, which is given by

[TABLE]

where we recall that we have ignored the temperature dependence of the parameters. The correction to the electrolyte potential satisfies (33b) and (35b), and can be found in a similar manner by writing $\Phi_{e,\text{{i}}}^{(1)}=[(\mathcal{I}^{(0)}\mathcal{D}_{\text{{A}}})/(2\phi_{\text{{e}}})](\Phi_{e,\text{{i}}}^{\infty}+\Psi_{e})$, where $\Psi_{e}(x,t)=2(1-\theta^{0})(\mathcal{D}_{\text{{L}}}/\mathcal{D}_{\text{{A}}}-1)[C_{L}(1,t)-C_{L}(x,t)]$. The solutions for the correction to the solid-phase potential is given by (68); however, the constants of integration $P_{p}$ and $P_{n}$ that appear in these expressions must now be equal in order to satisfy the voltage condition $\Phi_{a,\text{{p}}}^{(1)}(0,t)-\Phi_{a,\text{{n}}}^{(1)}(1,t)=0$. Thus, we let $P_{p}=P_{n}=P$. Solvability conditions that are analogous to (69) can be derived by integration of (32c), and these provide two equations for the correction to the current $\mathcal{I}^{(1)}$ and the constant $P$ in the solid-phase potentials,

[TABLE]

where $g_{\text{{i}}}^{(1)}$ is given by (64) with $\eta_{\text{{i}}}^{(0)}=\Phi_{a,\text{{i}}}^{(0)}$, $\eta_{\text{{i}}}^{(1)}=\Phi_{a,\text{{i}}}^{(1)}-\Phi_{e,\text{{i}}}^{(1)}-\mathcal{V}_{\text{{i}}}^{(1)}$, and using the expressions in (75). Having, in principle, determined $\mathcal{I}^{(1)}$, it can be shown that the corrections to the (integrated) concentration of intercalated lithium satisfy

[TABLE]

and thus do not vanish like in the galvanostatic case. The correction to the temperature satisfies

[TABLE]

with $u^{(1)}(0)=0$.

4.2.3 Dynamics on the saturation time scale

Following the analysis in Sec. 4.1, on the saturation time scale we write $t=\gamma_{c}^{-1}\hat{t}$, $c_{a,\text{{i}}}=\gamma_{c}^{-1}\hat{c}_{a,\text{{i}}}$, and use hats on the all of the other variables. The analysis is very similar to that of the galvanostatic case; however, we must now account for the fact that the leading-order contribution to the current density, given by $\hat{\mathcal{I}}^{(0)}$, is now time dependent. As discussed in Sec. 4.4, this only leads to minor modifications in the solutions. At leading-order, the solid-phase potentials and lithium concentrations are uniform. We thus have $\hat{\Phi}_{e,\text{{i}}}^{(0)}=0$, $\hat{\Phi}_{a,\text{{p}}}^{(0)}=\hat{\Phi}_{a,\text{{n}}}^{(0)}+\Delta v$, and find that the leading-order electrical problem is

[TABLE]

where $g_{\text{{i}}}^{(0)}$ is given by (56) and the leading-order concentrations of intercalated lithium evolve according to

[TABLE]

The leading-order current densities are given by (58) and (59) and the corresponding solution to the thermal problem is

[TABLE]

To derive the correction, we first note that in this time regime, both the liquid-phase concentration $\hat{c}_{L,\text{{i}}}^{(0)}$ and potential $\hat{\Phi}_{e,\text{{i}}}^{(1)}$ are known and given by their quasi-steady profiles (60) and (67), respectively. The corrections to the solid-phase potential are again given by (68) with $P_{p}=P_{n}=P$ and the solvability conditions (79) still apply. However, in this case the electrochemistry is more tightly coupled and the solution to the electrical problem cannot be obtained without the correction to the solid-phase concentration. This is because the expressions for $\mathcal{V}_{\text{{i}}}^{(1)}$ and $j_{\text{{i}}}^{(1)}$ are given by (64) and the contributions to these from $\hat{c}_{a,\text{{i}}}^{(1)}$ do not integrate to zero unlike in the galvanostatic case analysed in Sec. 4.1. Thus, obtaining the correction amounts to simultaneously solving the differential-algebraic system given by (79)–(80), which determines the constant $P$ in the solid-phase potentials, the correction to the current $\mathcal{I}^{(1)}$, and the correction to the net concentration of intercalated lithium. After the electrochemical problem has been solved, the correction to the temperature can be evaluated via the expression

[TABLE]

4.2.4 Construction of the composite solution

The leading-order composite asymptotic solution to the electrochemical problem is straightforward to obtain because this is given by the leading-order solution on the saturation time scale, i.e. the solution to (82) along with (83). The leading-order composite solution to the thermal problem is given by

[TABLE]

where $u^{(0)}$ solves (77), $\hat{u}^{(0)}$ is given by (84), and $u_{\text{overlap}}^{(0)}$ corresponds to the temperature in the overlap region, which is given by the steady-state solution to (77),

[TABLE]

Constructing the composite solution for the leading- and first-order problems is more involved but can be done using the solutions in the overlap region given in C.

4.3 Comparison with numerics

We simulate the potentiostatic case using the parameters in (46) and (47) which corresponds to a (dimensional) open-circuit potential of approximately $\Delta V=3.47$ V. Therefore we take $\Delta V=3.45$ V and $\Delta V=3.49$ V to discharge and charge the cell respectively with the results in Fig. 3 for the C-rate and temperature. The asymptotic and numerical solutions are compared to the simulations, with the focus being on the diffusive and saturation time regimes. The dynamics in the capacitance time regimes are shown in Fig. 9 (a)–(b) of B.

Figure 3 (a)–(b) shows the C-rate and cell temperature using the leading-order composite solutions and numerical simulations of the VA and P2D models. As with the galvanostatic case at these parameters, the leading-order correction is quite accurate with a near-perfect match when the first-order correction is considered. However, the agreement between the VA and P2D models has now diminished. As previously stated, the two models can diverge when the system is close to a singularity in the open-circuit potential due to the sensitive dependence on the surface concentration. For the simulations in Fig. 3, the solid-phase concentrations do not vary significantly and thus the system remains close to the singular points of the open-circuit potential. However, agreement holds for $t\ll 1$ because the solid concentration stays at the initial value while agreement for $t\gg 1$ is due to the current decreasing to zero leading to small surface concentrations.

If we alter the initial conditions to be away from the singular points, the disagreement between the VA and P2D models should diminish. Therefore, we take the initial concentration of intercalated lithium to be $c_{a,\text{{p}}}^{0}=0.39\,c_{a,\text{{p}}}^{\text{max}}$ and $c_{a,\text{{n}}}^{0}=0.43\,c_{a,\text{{n}}}^{\text{max}}$. These values can be obtained by simulating a discharge at a low C-rate for half of the total discharge time. We note that by changing these initial concentrations, the parameters in (46) and (47) associated with $c_{a,\text{{i}}}^{0}$ and $g_{\text{{i}}}^{0}$ change as well. At these concentrations, the (dimensional) resting potential of the cell is about $\Delta V=3.33$ V. We then prescribe cell potentials of $\Delta V=3.30$ V and $3.35$ V. We once again compare the asymptotic and numerical solutions for the C-rate and temperature, displaying the results in Fig. 4 (a)–(b). The dynamics in the capacitance time regimes are shown in Fig. 9 (c)–(d). We see that the error between VA and P2D has indeed been reduced and excellent agreement is found.

Interestingly, the matching between the VA and P2D models has come at the expense of a disagreement in the leading-order asymptotics. This is due to the increased role of the electrolyte at small values of the open-circuit potential. If we instead compare with a composite solution including the leading- and first-order contributions (Fig. 4 (c)–(d)) then we recover excellent agreement between the three cases.

The necessary inclusion of first-order corrections in Fig. 3 for agreement between the asymptotic and numerical solutions can be rationalised by analysing (76a), which can be interpreted as current balance involving two resistors connected in parallel, where the parameters $\mathcal{G}_{\text{{i}}}$ act as inverse resistances (i.e. conductivities). Thus, (76a) describes how the electrical load is distributed across the two electrodes and determines where the drop in potential occurs. Recalculating the values of $\mathcal{G}_{\text{{i}}}$ using the new values for the initial solid-phase concentrations reveals that $\mathcal{G}_{\text{{p}}}x_{p}\simeq 9.5$ but $\mathcal{G}_{\text{{n}}}(1-x_{n})\simeq 1.8$, in contrast to the previous values $\mathcal{G}_{\text{{p}}}x_{p}\simeq 2.8$ and $\mathcal{G}_{\text{{n}}}(1-x_{n})=1.3$. For the new initial conditions, most of the potential drop occurs in the negative electrode and the electrical potential in the positive electrode remains close to zero. Consequently, the higher-order terms which have been neglected from the electrical problem in the positive electrode play much a stronger role. For the original parameters, the load is equally distributed across both electrodes.

4.4 Galvanodynamic charging/discharging

The analysis in Sec. 4.1 for a galvanostatic discharge can be extended to the case when the C-rate is a function of time. If the variations in current occur over the saturation time scale, then the asymptotic solutions presented in Sec. 4.1 can be used by replacing the expressions for the solid-phase concentrations (54) with the differential equations (83). If the C-rate varies on the diffusive time scale, then the framework developed in Sec. 4.2.2 can be applied. One caveat, however, is that the transient solution for the liquid-phase concentration (78) needs to be adapted to account for a time-dependent current if the first-order correction is used.

To compare the asymptotic solutions with numerical simulations in the case of a time-dependent applied current, we consider an oscillating C-rate of the form $\mathcal{I}(t)=A\sin(2\pi t/\tau)$, which is used in electrochemical impedance spectroscopy. We fix the amplitude at $A=2$ so that the maximum (dis)charge rate is 2C. We consider three sets of initial conditions for the solid-phase concentrations, each of which corresponds to a different initial state of charge (SOC). We define the SOC to be the fraction of total capacity that is available to use. To generate new initial conditions, we solve the leading-order galvanostatic model with the initial conditions in Table 3 to full discharge, which define the 100% and 0% SOC points, respectively. We then determine the solid concentrations in each electrode at a given SOC along the voltage curve and use these values to re-initialise the model. Since the solid concentrations at a given SOC are dependent on the initial conditions used in solving the galvanostatic model, there is not a one-to-one correspondance between a given SOC and the new initial solid concentrations. The new initial conditions we take are: $c_{a,\text{{p}}}^{0}=0.21\,c_{a,\text{{p}}}^{\text{max}}$, $c_{a,\text{{n}}}^{0}=0.64\,c_{a,\text{{n}}}^{\text{max}}$ (75% SOC); $c_{a,\text{{p}}}^{0}=0.39\,c_{a,\text{{p}}}^{\text{max}}$, $c_{a,\text{{n}}}^{0}=0.43\,c_{a,\text{{n}}}^{\text{max}}$ (50% SOC); and $c_{a,\text{{p}}}^{0}=0.58\,c_{a,\text{{p}}}^{\text{max}}$, $c_{a,\text{{n}}}^{0}=0.21\,c_{a,\text{{n}}}^{\text{max}}$ (25% SOC). We recall that each time the initial conditions are changed, the parameters associated with $c_{a,\text{{i}}}^{0}$ and $g_{\text{{i}}}^{0}$ in (46) and (47) are changed as well. The period of oscillation is chosen to be either $\tau=1200$ s or $\tau=60$ s. In the former, the current varies on the saturation time scale, in the latter, it varies on the diffusive time scale.

The results are shown in Fig. 5, which plots the cell potential and temperature for the different SOC values using the two different frequencies. In the low-frequency case ($\tau=1200$ s; panels (a)–(b)), the asymptotic solution is constructed using the leading- and first-order contributions. The agreement between all three models is excellent. The temperature exhibits an oscillatory behaviour as well as a result of the electrochemical heat generation being dominated by the reversible contribution. For the high-frequency case ($\tau=60$ s; panels (c)–(d)), the asymptotic solution is only based on the leading-order expressions for simplicity but it still leads to an accurate approximation of the numerical solution.

5 Upscaling to a multi-cell model

While a single-cell model provides important insights into the coupling that occurs between the electrochamical and thermal processes, a real battery consists of many cells that are coupled together. The heat that is generated within a cell can accumulate due to the cooled exterior of the battery being further away, which can cause a greater increase in temperature than predicted by the single-cell model. To capture this behaviour, we now consider a multi-cell model consisting of repeating LIB cells that are separated by metal current collectors (Fig. 6). Current collectors are required for a multi-cell model to prevent short-circuiting between adjacent electrodes. Using homogenisation theory, we can derive a single heat equation that holds across the array of cells, where the volumetric heating term is calculated using the cell-scale asymptotic solutions to the electrochemical problem derived in the previous section.

Homogenisation of the heat equation requires a consideration of the current collectors, which were neglected in the single-cell model. The non-dimensional half-thickness of the current collectors adjoining the positive and negative electrodes are denoted by $\lambda_{\text{{p}}}^{\text{cc}}$ and $\lambda_{\text{{n}}}^{\text{cc}}$, respectively. Each current collector is adjacent to the same type of electrode, as shown in Fig. 6. Thus, for periodicity on the cell scale, we take the cell problem to involve two half current collectors for the positive electrode, two positive and negative electrodes, two separators, and one current collector for the negative electrode. The total length of this geometry is denoted by $2L_{c}$ where $L_{c}=1+\lambda_{\text{{p}}}^{\text{cc}}+\lambda_{\text{{n}}}^{\text{cc}}$. The centerlines of the positive-electrode current collectors are taken to be at $x=-\lambda_{\text{{p}}}^{\text{cc}}$ and $x=2(1+\lambda_{\text{{n}}}^{\text{cc}})+\lambda_{\text{{p}}}^{\text{cc}}$. At the edge of the periodic array, each end is connected to a positive half current collector. The dimensionless volumetric heat capacity and thermal conductivities of the current collectors are denoted by $\varrho_{\text{{i}}}^{\text{cc}}$ and $\mathcal{K}_{\text{{i}}}^{\text{cc}}$. Although the current collectors can generate heat through Ohmic heating, their large electrical conductivities makes this negligible compared to the electrochemical reaction heat [22].

The LIB is assumed to consist of $N\gg 1$ (physical) cells, where $N$ is even so that there are $N/2$ unit cells. To carry out the multiple-scales expansion, we introduce the slow variable $X=\epsilon(x+2\lambda_{\text{{p}}}^{\text{cc}})$, where $\epsilon=(NL_{c}+2\lambda_{\text{{p}}}^{\text{cc}})^{-1}\ll 1$, which describes spatial variations across the battery scale. This definition of $X$ and $\epsilon$ is chosen so that $0\leq X\leq 1$. We are interested in the thermal diffusion over the length of the array and therefore the analysis is performed in the distinguished limit where $\mathrm{Le}=\epsilon^{-2}\overline{\mathrm{Le}}$ with $\overline{\mathrm{Le}}=O(1)$. However, the resulting model can also be used to study the subcases $\overline{\mathrm{Le}}\ll 1$ and $\overline{\mathrm{Le}}\gg 1$. The details of the homogenisation procedure are given in Ref. [71] so we will only briefly discuss them here. It is convenient to group the five heat capacities and thermal conductivies into piecewise-constant functions, $\varrho(x)$ and $\mathcal{K}(x)$, that vary on the cell scale. After neglecting sub-dominant Ohmic heat terms, the non-dimensional heat equation describing the temperature in the “unit cell” is given by

[TABLE]

with $q_{r}=q_{r}(x,X,t)$ denoting a piecewise function describing the electrochemical heat generation in each electrode. The Newton cooling conditions (44c) and (44o) are now applied at $X=0$ and $X=1$ while periodic boundary conditions are imposed at the edges of the unit cell. Continuity of temperature and flux is maintained at the electrode-collector boundaries. The temperature is expanded as $T=T^{(0)}+\epsilon T^{(1)}+\epsilon^{2}T^{(2)}+O(\epsilon^{3})$. At leading-order we find that $T^{(0)}=T^{(0)}(X,t)$ and thus the temperature is constant across individual cells, in agreement with the analysis in Sec. 3.1. From the $O(\epsilon)$ problem we can deduce that

[TABLE]

which corresponds to an effective form of Fourier’s law based on the harmonic mean of the thermal conductivity. The $O(\epsilon^{2})$ problem leads to the battery-scale model for the temperature $T^{(0)}$, which is given by

[TABLE]

where $\overline{\varrho}=L_{c}^{-1}\int_{-\lambda_{\text{{p}}}^{\text{cc}}}^{1+\lambda_{\text{{n}}}^{\text{cc}}}\varrho(x)\,\mathrm{d}x$. The averages of the electrochemical heating terms, $\langle q_{r,\text{{p}}}\rangle$ and $\langle q_{r,\text{{n}}}\rangle$, are defined in (49). However, these heating terms are now functions of the battery-scale variable $X$ through their dependence on the temperature $T^{(0)}$, which leads to all of the cell-scale electrochemical problems being coupled together.

As all of the cell components are effectively lumped together on the battery scale, and only the positive-electrode current collectors are in contact with the environment, the correct form of the Newton cooling conditions are

[TABLE]

where $\overline{\mathrm{Bi}}=\epsilon^{-1}\mathrm{Bi}$.

5.1 Heat generation during a galvanostatic discharge

To illustrate the thermal response of several connected LIB cells, we consider a situation where a constant current $\mathcal{I}_{\text{tot}}$ is drawn from the battery. The cells are assumed to be identical and connected in parallel. A similar configuration has been used by An *et al. *[43] to model heat generation in a commercial LIB. The potential difference $\Delta v(t)$ across all of the cells will be the same but the local applied current density (i.e. the local C-rate) $\mathcal{I}(X,t)$ can vary between cells. Thus, currents must satisfy

[TABLE]

Since the electrical potentials are approximately uniform on the cell scale, we also have that $\Phi_{a,\text{{p}}}(X,t)-\Phi_{a,\text{{n}}}(X,t)=\Delta v(t)$.

The leading-order analysis in Sec. 4.1 is used to construct and simplify the thermal-electrochemical problem. By writing $t=\gamma_{c}^{-1}\hat{t}$, $c_{a,\text{{i}}}=\gamma_{c}^{-1}\hat{c}_{a,\text{{i}}}$, and taking the limit as $\mathcal{C}_{\text{{i}}}\to 0$ and $\gamma_{c}\to 0$, we find that the corresponding problem is given by the following system of partial-differential-integral equations

[TABLE]

where the solid-phase concentrations are given by (83) and electrochemical kinetics are given by (56a)–(56b) with $K_{\text{eff},\text{{i}}}^{(0)}=\exp\left(\Gamma_{\text{eff},\text{{i}}}T^{(0)}\right)$. Equations (93a) and (93b) determine the local potential in the negative electrode and the local current density. Equation (93c) is a global constraint that determines the potential drop across the battery. Finally, Eqn (93d) describes the evolution of the local temperature.

The spatial variation that occurs in the electrical problem across the battery scale is due to gradients in the temperature. However, the effective Biot number $\overline{\mathrm{Bi}}$ will generally be small and thus temperature gradients will be weak. In this case, it is possible to carry out an asymptotic expansion in terms of $\overline{\mathrm{Bi}}\ll 1$ to find that the leading-order contributions are independent of $X$ and thus only functions of time. The leading-order battery temperature is given by

[TABLE]

Equation (94) is the battery-scale analogue of (61) and it shows that the temperature rise scales roughly with the number of cells $N$. Using $\mathcal{I}=\mathcal{I}_{\text{tot}}$, the battery potential $\Delta v$ can be obtained be solving (93a)–(93b) and the solid-phase lithium concentrations are given by (54).

We numerically solve the battery problem (93) for a 1C discharge and explore how the number of cells $N$ affects the evolution of the temperature. The results are shown in Fig. 7. In all the cases considered, the effective Biot number $\overline{\mathrm{Bi}}$ is small and so the temperature is approximately uniform across the battery pack. As expected, a battery that consists of more cells becomes hotter during discharge. For a 60-cell battery, the temperature rise is roughly 20 K, which is close to the threshold at which Arrhenius reaction kinetics begin to play a role.

In Fig. 8, we focus our attention on the heat generation in a 60-cell battery at higher discharge rates extending up to 4C. The temperature increase in these cases, shown in Fig. 8 (a), can be substantial and exceed 50 K. The corresponding voltage curves are plotted as solid lines in Fig. 8 (b), with the isothermal curves shown as dashed lines for reference. The increase in battery temperature leads to a greater cell potential relative to the isothermal case. This is due to the increase in the effective reaction rate. There is less electrical resistance caused by the finite rate of electrochemical reactions driving (de)lithiation and thus the cell potentials are pushed closer to the open-circuit potentials. Such behaviour was predicted by analysing the thermal-electrochemical problem at the cell level in Sec. 4.1.1. Importantly, we see that even with a temperature rise of roughly 60 K and the onset of Arrhenius-dominated reaction kinetics, the model does not predict the occurrence of thermal runaway. We emphasise that the rapid temperature rise near the end of the discharge process is associated with the sharp change in potential caused by the electrodes becoming fully saturated or fully depleted of lithium.

6 Discussion and conclusion

We present thermo-electrochemical models for an LIB which treat the electrode microstructure in two different ways. In one case, volume averaging is used to homogenise the microstructure. In the second, the electrode is treated as a collection of spherical particles following the P2D framework. The two approaches become equivalent in the limit of fast diffusion in the solid particles. The VA model is mathematically simpler and, consequently, asymptotic solutions can be obtained from it for common modes of battery operation. The VA model for a single LIB cell is then homogenised to obtain a model for a battery consisting of $N\gg 1$ cells.

In the case of a single LIB cell, numerical simulations show strong agreement between the VA and P2D models for the moderate discharge rates considered here. The asymptotic solutions give excellent approximations to the VA model and thus to the P2D model as well. The P2D framework has become a standard modelling approach in the battery community and the results presented here show that the simpler and computationally advantageous VA model can yield the same predictions of macroscopic quantities such as cell temperature and potential. However, the loss of information about the microstructure can prevent the VA model from being used in applications involving phase-separating electrodes, degradation, and active material utilisation. Furthermore, the VA model is unable to capture the effective loss of capacity that occurs in some cells at high C-rates [72], which is attributed to non-uniform lithium distribution in the electrode particles [64]. Thus, the best model to use is situation dependent. Establishing the conditions under which the P2D model can be asymptotically reduced to the VA model can serve as a useful guide for model selection and is a key area of future research.

An important area of battery modelling concerns model validation. Typically, models are compared with experimental voltage curves that are obtained by subjecting a battery (or a cell) to a given current. The good agreement between the leading-order asymptotic solutions and the numerical simulations of the full models in the case of fixed or alternating currents at moderate C-rates indicates that spatial effects and the electrolyte play almost no role in determining the cell potential. Thus, these aspects of the model may be difficult to experimentally validate using galvanostatic or galvanodynamic modes of operation. However, these models are ideal for validating the thermal aspects of the model because a sustained current held at a moderate C-rate is required to raise the temperature of the battery to a detectable level. In the case of a potentiostatic hold, whereby the cell potential is fixed, obtaining good agreement between the asymptotic and numerical solutions for the current required the first-order contributions to be included, which captures effects relating the spatial gradients and the electrolyte. The current obtained from potentiostatic experiments therefore encodes useful information about the dynamics of the electrolyte which can assist with model validation.

The scaling analysis revealed that the electrochemical reactions are the dominant source of heat generation. Given the Arrhenius dependence of the reaction coefficients and the large temperature rises that can occur when many LIB cells are coupled together, it is surprising to find that thermal runaway is not possible in our model. On one hand, this is reassuring, as it shows that LIBs are generally safe and thermal runaway must be due to mechanical abuse or poor design. On the other hand, our model has neglected the fact that above a critical temperature, the electrolyte, active materials, and the solid electrolyte interphase will decompose. These reactions are exothermic and generate flammable gases, which can lead to catastrophic battery failure. In other LIB designs, Ohmic heating in the electrolyte can be a strong contributor to the total heat generation [22]. The framework developed here can easily accommodate this case: the asymptotic solutions for the liquid-phase concentrations and current densities can be used to calculate the component-averaged Ohmic heating terms in reduced heat equations such as (90). Our reduction framework can be extended even further to account for empirical expressions for the open-circuit voltage and its derivative with respect to temperature. Asymptotically reduced and homogenised models such as those derived here are not only adaptable, but computationally inexpensive, making them ideal for use in on-board thermal management systems, thus paving the way towards safer battery operation.

Acknowledgements

IRM gratefully acknowledges funding from the Charlemont scholar program of the Royal Irish Academy as well as an NSERC Discovery Grant (2019-06337). MGH has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No. 707658. The authors thank Brian Wetton, Tim Myers, Jon Chapman, Colin Please, and Mohit Dalwadi for insightful discussions relating to this work.

Appendix A Parameter determination

The parameters values are adapted from Refs. [65, 22], in which thermal P2D models are developed for ANR26650m1-a lithium-ion cells [69]. The negative and positive electrodes of these cells are constructed from graphite and LiFePO4 (LFP), respectively. The physical dimensions of the cell, along with its capacity and common parameters, are given in Table 3. Values for the electrochemical and thermal parameters are listed in Tables 3 and 3, respectively. The parameters $I_{\rm app}$ and $A_{\rm cell}$ define the current density $i_{0}=I_{\rm app}/A_{\rm cell}$.

The reaction constants at ambient temperature, $K_{a,\text{{i}}}^{0}$ and $K_{L,\text{{i}}}^{0}$, are determined by fitting the analytical expression for the open-circuit potential (16) to the empirical expressions for graphite and LFP given by

[TABLE]

where $y_{\text{{i}}}$ is the state of charge of the battery defined as $y_{\text{{i}}}=c_{a,\text{{i}}}/c_{a,\text{{i}}}^{\text{max}}$. Having determined $K_{a,\text{{i}}}^{0}$ and $K_{L,\text{{i}}}^{0}$, the difference in activation energies $E_{L,\text{{i}}}-E_{a,\text{{i}}}$ is then obtained by fitting the analytical expression for $\partial U_{\text{{i}}}/\partial T$ obtained from (16) to empirical expressions given by

[TABLE]

The diffusivitity of lithium ions and anion in the electrolyte, $D_{L}^{0}$ and $D_{A}^{0}$, are obtained by simultaneously solving (10) and (11) given known values of the effective diffusivity $D_{e}^{0}$ and the transference number $\theta^{0}$ at ambient temperature.

Appendix B Analysis of capacitance time regimes for the potentiostatic case

The first capacitance regime for the case of a galvanostatic discharge is defined by $t=\mathcal{C}_{\text{{n}}}\tilde{t}$ for C-rates that are $O(1)$ in size due to the singular nature of the time derivative in (32). For the potentiostatic problem, a very small-time analysis demonstrates that $\Phi_{a,\text{{i}}}(x)=\Phi_{e,\text{{i}}}(x)$, corresponding to zero initial overpotential. If the C-rate is again $O(1)$ then these potentials are spatially independent. The grounding condition implies that $\Phi_{a,\text{{i}}}=\Phi_{e,\text{{i}}}=0$, which prevents the condition on the cell potential, $\Phi_{a,\text{{p}}}-\Phi_{a,\text{{n}}}=\Delta v$, from being satisfied. The potentials at very small times must therefore be spatially inhomogeneous and this suggests that there is a distinguished small-time limit with large current densities. This large current is due to the lack of charge-transfer resistance at extremely small times. Therefore, to capture the dynamics in the first capacitance regime, we write $t=\tau_{n}\tilde{t}$ where $\tau_{n}=\mathcal{C}_{\text{{n}}}\mathcal{G}_{\text{{n}}}\nu_{a,\text{{n}}}/\phi_{a,\text{{n}}}$. The C-rate, and thus the current densities, scale with $\langle\nu\rangle_{\text{cell}}^{-1}\Delta v$, where $\langle\nu\rangle_{\text{cell}}\ll 1$ is defined by (71). On this time scale, the change in solid- and liquid-phase lithium concentration scales like $\Delta c\sim\tau_{n}\langle\nu\rangle_{\text{cell}}^{-1}\Delta v$ and will be small provided that $\Delta v=O(1)$. The large current amplifies Ohmic heating, which leads to temperature changes that scale like $\Delta T\sim\tau_{n}\mathcal{Q}_{a,\text{{i}}}(\langle\nu\rangle_{\text{cell}}^{-1}\Delta v)^{2}$ for the electrode contributions and $\Delta T\sim\tau_{n}\mathcal{Q}_{\text{{e}}}(\langle\nu\rangle_{\text{cell}}^{-1}\Delta v)^{2}$ for the electrolyte. Electrochemical reactions lead to changes in temperature that scale like $\Delta T\sim\tau_{n}$. Based on these scaling estimates, we can assume that neither the concentration nor the temperature changes in this time regime, and thus we only need to consider the electrical problem. By eliminating the current, the electrical problem takes the form

[TABLE]

where $\ell=[(\phi_{a,\text{{p}}}\nu_{a,\text{{n}}}\mathcal{G}_{\text{{n}}})/(\phi_{a,\text{{n}}}\nu_{a,\text{{p}}}\mathcal{G}_{\text{{p}}})]^{1/2}\ll 1$, $\mathcal{C}=\mathcal{C}_{\text{{p}}}/\mathcal{C}_{\text{{n}}}\gg 1$, and $\varepsilon_{n}=\mathcal{G}_{\text{{n}}}\nu_{a,\text{{n}}}/\phi_{a,\text{{n}}}\ll 1$. Boundary conditions are given by

[TABLE]

In the limit $\tilde{t}\to 0$, we find that the solid- and liquid-phase potentials are equal and independent of time, $\Phi_{a,\text{{i}}}(x)=\Phi_{e,\text{{i}}}(x)$, satisfying the zero initial overpotential condition. The equality of solid- and liquid-phase potentials can be used in (99) to show that

[TABLE]

where the equivalent resistances are defined in (72). Note that (104) also applies to the separator. Integration of (104) across each of the respective domains, adding the results, and using continuity of the potential along with $\Phi_{e,\text{{p}}}(0)-\Phi_{e,\text{{n}}}(1)=\Delta v$, leads to the expression for the initial C-rate given by (71). As time increases (i.e. for $\tilde{t}=O(1)$), the gradient in the solid-phase current in the negative electrode (the left-hand side of (98)) drives an overpotential (the second term on the right-hand side of (98)), resulting in the onset of electrochemical reactions. The increase in electrical resistance due to charge transfer across the electrode-electrolyte interface leads to a marked decrease in the current. For even longer times (i.e. for $1\ll t\ll\varepsilon_{n}^{-1}\mathcal{C}$), the current will relax an $O(1)$ value which reaches a plateau when the time derivative in (98) vanishes. When this happens, the current in the negative electrode is once again obtained from the electrochemistry through the equation $g_{\text{{n}}}=\mathcal{I}/[\mathcal{G}_{\text{{n}}}(1-x_{n})]$. We can isolate the cell potential from this expression:

[TABLE]

To leading order in $\nu_{a,\text{{p}}}$ and $\nu_{\text{{e}}}$, we find $\Phi_{e,\text{{p}}}=\Phi_{a,\text{{p}}}=0$. However, we can easily obtain the corrections by noting that (104) still holds for the positive electrode and separator which can once again be integrated across each domain. We couple this with $\partial\Phi_{e,\text{{n}}}/\partial x=-\nu_{\text{{e}}}i_{e,\text{{n}}}=\nu_{\text{{e}}}\mathcal{I}[1-(x-x_{n})/(1-x_{n})]$ to produce $\Phi_{a,\text{{p}}}(0,t)=-\mathcal{I}[\nu_{\mathrm{eq},\text{{p}}}x_{p}+\nu_{\mathrm{eq},\text{{s}}}(x_{n}-x_{p})+\nu_{\mathrm{eq},\text{{s}}}/2(1-x_{n})]$ yielding the plateau value (73).

In the second capacitance regime, we write $t=\tau_{p}\check{t}$, where $\tau_{p}=\mathcal{C}_{\text{{p}}}$ and take $\mathcal{C}\to\infty$ to obtain

[TABLE]

Due to the smallness of $\nu_{a,\text{{i}}}$ and $\nu_{\text{{e}}}$, the potentials in both electrodes are approximately uniform in space. Thus, the solid- and liquid-phase electrical potentials in the electrodes are related through the equations $g_{\text{{n}}}=\mathcal{I}/[\mathcal{G}_{\text{{n}}}(1-x_{n})]$ and $g_{\text{{p}}}=\mathcal{I}/(\mathcal{G}_{\text{{p}}}x_{p})$ while the currents are given by the expressions in (58)–(59). Taking the limit $\check{t}\gg 1$ yields

[TABLE]

The positive electrolyte potential satisfies $\Phi_{e,\text{{p}}}(0,t)=0$ to leading order, but the correction can be computed similarly to $\Phi_{a,\text{{p}}}(0,t)$ in the first capacitive regime producing

[TABLE]

resulting in the plateau value given by (74).

The transient evolution of the C-rate in the two capacitance time regimes is shown in Fig. 9. The curves are obtained by solving (99) using a fully implicit method based on a spectral discretisation in space and Newton iterations at each time step. The initial conditions are the same as those discussed in Sec. 4.3. Therefore, the curves in Fig. 9 capture the early dynamics that are not shown in Fig. 4. The capacitance time scales are $\tau_{n}\simeq 3\times 10^{-5}$ s and $\tau_{p}\simeq 0.3$ s for panels (a)–(b) and $\tau_{p}=0.1$ for panels (c)–(d). The three plateaus described by (71), (73), and (74) as indicated by dashed lines.

Appendix C The composite solution in the case of a potentiostatic hold

We first construct the composite solution to the leading- and first-order electrical problem and then consider the solution to the thermal problem. To simplify the notation, we write the steady-state liquid-phase concentration and potential as

[TABLE]

The large-time expansions of the corrections to the open-circuit potential $\mathcal{V}_{\text{{i}}}^{(1)}$ and exchange current density $j_{\text{{i}}}^{(1)}$ on the diffusive time scale are given by $\mathcal{V}_{\text{{i}}}^{(1)}\sim M_{V,\text{{i}}}t+N_{V,\text{{i}}}(x)$, $j_{\text{{i}}}^{(1)}\sim M_{j,\text{{i}}}t+N_{j,\text{{i}}}(x)$ as $t\to\infty$, where

[TABLE]

where we have used the fact that the leading-order solid-phase lithium concentrations $c_{a,\text{{i}}}^{(0)}$ are linear functions of time. Consequently, the constant $P$ that appears in the correction to the solid-phase potentials $\Phi_{a,\text{{i}}}^{(1)}\sim\Phi_{a,\text{{i}}}^{\infty}+P$ can be written as $P\sim M_{P}t+N_{P}$ as $t\to\infty$, where

[TABLE]

The correction to the C-rate can also be expanded as $\mathcal{I}^{(1)}\sim M_{I}t+N_{I}$ as $t\to\infty$, where

[TABLE]

The solution to the C-rate in the overlap region is found by matching to the solutions on the saturation time scale in the limit as $\hat{t}\to 0$. We find that the leading-order contribution is given by $\mathcal{I}^{(0)}_{\text{overlap}}=\mathcal{I}^{(0)}$, where $\mathcal{I}^{(0)}$ is determined from (76). The first-order contribution is given by $\mathcal{I}^{(1)}_{\text{overlap}}=M_{I}t+N_{I}$.

To construct the composite solution for the temperature, we use the fact that $u^{(1)}\sim M_{u}t+N_{u}$ as $t\to\infty$, where

[TABLE]

By matching to the asymptotic solution for the tempetaure on the saturation time scale in the limit as $\hat{t}\to 0$, we find that the first-order contribution in the overlap region is given by $u^{(1)}_{\text{overlap}}=M_{u}t+N_{u}$.