1 Introduction

The role of rechargeable batteries has become more and more important in recent years due to the increase in the use of electronic devices and electric vehicles. In particular, most of these applications rely on lithium ion batteries, a type of porous electrode battery, and therefore, research on this type of battery has received a lot of attention in the past few years. The first electrochemical model for porous electrode batteries was put forward in the seminal paper of Doyle et al. [1]. In that article, an effective model was presented which described the mass and charge transport in both the electrode and the electrolyte. The porous structure of the electrode was modelled by assuming that at each point of the electrodes there is a representative particle in which lithium intercalates and diffuses. This model is commonly known as the pseudo-two-dimensional (P2D) model, because of the two scales involved (electrode particles and cell scales) or as the Doyle, Fuller & Newman (DFN) model. In this paper, we will refer to it as the DFN model. This model has been since the backbone of lithium ion battery modelling and, more generally, of porous electrode batteries. As most of the modelling efforts on porous electrode batteries have been driven by lithium ion batteries, most of the references and notation in this article refer to that particular application, despite that the results presented are more general.

In the models for porous electrode batteries, the equations derived from the thermodynamics are at the microscale. This means that they need to be solved in a very complex geometry in order to capture the porous electrode, which requires a lot of computational power. For this reason effective models, such as the DFN model, have been popular: they are posed at the macroscale and therefore are computationally cheaper, but still retain most features of the microstructure and the microscale dynamics. One way to obtain effective macroscale equations from the microscale equations is using the volume averaging method, like in [1,2,3]. This method provides good results, but it requires to either define the effective model ad hoc [1] or neglect some terms in the averaging [3]. In both cases, the parameters of the macroscale model need to be defined as effective parameters and, therefore, fitted from data. This means that the macroscale parameters cannot be directly and systematically derived from the microstructure, and has led to the use of empirical correlations which are still a subject of debate in the field, such as the Bruggeman correlation [4]. For more details on volume averaging we refer the reader to [5].

Another approach to upscale the microscale equations is to use asymptotic homogenisation. The main difference between volume averaging and homogenisation is that the latter allows the effective parameters from the microstructure to be derived in a systematic manner. The key idea of asymptotic homogenisation is to assume that the material is composed of a periodic structure of a size much smaller compared to that of the material. Then, one can exploit the disparity of scales and define space variables at each level that can be treated as independent, so all the differential operators can be split into operators at each scale and a small parameter (ratio of length scales) appears. By performing an asymptotic expansion in this small parameter, an equation at the large scale can be determined, which accounts for the effects occurring at the small scale. Homogenisation is a well-known technique: the technical details can be found in the handbooks of the subject [6, 7] and a detailed comparison between volume averaging and asymptotic homogenisation is provided in [8].

In our case, we have a problem in which there is transport of mass around the particles (i.e. the microstructure), with a chemical reaction at the surface of the particles and diffusion inside them. This type of problem appears in many different applications, apart from porous electrode batteries, such as filtration [9, 10], biology [11], metallurgical furnaces [12] or even coffee roasting [13]. In order to capture the diffusion in the particles that form the microstructure, it is necessary to use high-contrast homogenisation as the diffusion coeffient in the particles is much smaller than the diffusion coefficient in the medium surrounding them [14, 15]. High-contrast homogenisation is not as developed as standard homogenisation, and it presents numerous challenges, especially on the analysis side [16,17,18].

Thermal effects can have a notable impact on the behaviour of batteries, and therefore, they are an important aspect of our model. The thermal model is coupled to the electrochemical model in the following way: the currents in the battery generate heat in several ways, while the temperature of the battery affects the parameters of the electrochemical model. The coupled thermal electrochemical model can then be upscaled from the microscale equations from a few different ways, similarly to the electrochemical model. The most common approaches in the literature are posing the model ad hoc [19] or using volume averaging [3, 20], which present the same advantages and disadvantages discussed for the purely electrochemical model. The goal of this paper is to upscale and derive the coupled thermal-electrochemical model systematically, using asymptotic homogenisation, and including some features that extend the models present in the literature. First, we keep the relation between the microstructure and the effective parameters and also consider a more general microstructure than the single particle at each point considered in the classic DFN model. We also retain all the features in the electrochemistry of the DFN model, in particular diffusion in the particles, which is important in the behaviour of the batteries, especially during the relaxation period after the current is switched off. Finally, we include the double layer capacitance effects in the kinetics between the electrode and the electrolyte, similarly to the model in [2].

Homogenisation was first used to derive an effective electrochemical model for porous electrode batteries by Ciucci and Lai in [21]. Assuming that the microstructure is formed by packed spheres they derive the DFN model from the microscale equations. They use high-contrast homogenisation so diffusion in the particles is retained.

A similar approach was taken by Richardson et al. in [22]. Here, the authors homogenise the electrochemical microscale equations to obtain effective macroscale equations and use the following assumptions: fast diffusion in the electrode particles, high electronic conductivity in the electrodes and dilute electrolyte. They present a very detailed analysis of the homogenisation and then perform an asymptotic analysis of the effective model to obtain analytical solutions.

In [23, 24], Arunachalam et al. perform asymptotic homogenisation of the electrochemical microscale model. The authors consider the diffusivity in the particles to be the same order of magnitude as the diffusivity in the electrolyte. Therefore, the effective model includes lithium diffusion over the whole thickness of the electrode, rather than within the particle. In [24] the authors also provide an analysis of the region in the parameter space in which the DFN model is valid.

Focusing on another aspect within batteries, Hennessy and Moyles [25] used homogenisation to derive the battery heat equation for a double coated electrode from the cell heat equation. The thermal properties of the cell problem are taken to be scalars and the connection to the microscale bulk properties is not considered.

The outline of this paper is as follows. In Sect. 2 we provide the homogenised dimensionless model for thermal and electrochemical behaviour of porous electrode batteries. The microscale equations which are the starting point of the derivation of the homogenised model are presented in Sect. 3 and the details of the homogenisation procedure are given in Sect. 4. In Sect. 5 we compare our homogenised model with the widely used DFN model. Finally, in Sect. 6 we discuss the results.

2 Dimensionless homogenised model

We use asymptotic homogenisation to derive a thermal-electrochemical model for porous electrode batteries both at cell and battery levels, as this method allows us to derive macroscale equations from the microscale ones in a systematic way. In this section the effective model is presented. This model can be seen as a generalised version of the Doyle, Fuller & Newman (DFN) model [1] and the physical laws underpinning it are discussed at the beginning of Sect. 3. Notice that, even though we use lithium ion battery terminology in the article, the model presented here can be applied to other types of porous electrode batteries.

The variables of the problem are the concentration of lithium in the electrode \(c_\mathrm {s}\) and ions in the electrolyte \(c_\mathrm {e}\), the electric potentials \(\varPhi _\mathrm {s}\) and \(\varPhi _\mathrm {e}\), and the temperature T. The subscripts \(\mathrm {s}\) and \(\mathrm {e}\) distinguish the quantities in the solid electrode and the electrolyte, respectively. For convenience, the currents \({{\varvec{i}}}_\mathrm {s}\) and \({{\varvec{i}}}_\mathrm {e}\) in the electrode and the electrolyte are defined. These are quantities derived from the concentrations and potentials.

Fig. 1
figure 1

Sketch of the geometry at the three length scales that are involved in the effective model: battery, cell, and particles. We label as well the regions at each scale. At the macroscale, where the macroscale heat equation (5) is defined, the domain is the battery \(\varOmega _{\mathrm {batt}}\), which is composed of several cells. At the mesoscale, where the electrochemical model (2)–(3) and the cell thermal model (4) are defined, the domains are the negative electrode \(\varOmega _\mathrm {n}\), the separator \(\varOmega _{\mathrm {sep}}\), and the positive electrode \(\varOmega _\mathrm {p}\) (apart from the collectors which have not been labelled as we do not consider them in the problem). The cell is defined as the union of the three parts so \(\varOmega _{\mathrm {cell}} = \varOmega _\mathrm {n} \cup \varOmega _{\mathrm {sep}} \cup \varOmega _\mathrm {p}\). At the microscale, where the diffusion in the particles (1) is defined, the domain is the representative microstructure of the porous structure \(\varOmega _\mathrm {s}\)

As shown in Fig. 1, the effective model accounts for phenomena occurring at three different scales. At the microscale, denoted by z, there are the electrode particles \(\varOmega _\mathrm {s}\). The mesoscale, denoted by y, represents a cell, which is composed of the negative electrode \(\varOmega _\mathrm {n}\), the separator \(\varOmega _\mathrm {sep}\), and the positive electrode \(\varOmega _\mathrm {p}\). The boundary of the domain needs to be separated into the part in contact with the current collectors (denoted by \(\partial \varOmega _\mathrm {n}^\mathrm {collector}\) and \(\partial \varOmega _\mathrm {p}^\mathrm {collector}\)) and the rest of the boundary, as different conditions apply to each one. At each point of the cell, there is the solid electrode and the electrolyte and thus different variables are used to measure them. The cell domain, defined as the union of the three parts, is \(\varOmega _\mathrm {cell}\). Finally, there is the macroscale, denoted by x, which is the battery, composed of several cells and represented by \(\varOmega _\mathrm {batt}\). In the differential operator \(\nabla \) subscripts x, y, and z are used to denote at which scale the operator is applied. Similarly, we use these subscripts in the volume and area integration variables to denote the scale.

Then, the homogenised dimensionless problem is the following. For the concentration of lithium in the electrodes the problem is defined at the microscale domain. It reads

$$\begin{aligned}&\partial _tc_\mathrm {s} =\nabla _z \cdot \left( D_\mathrm {s} \nabla _z c_\mathrm {s} \right) \quad \text { in }\quad x\varOmega _\mathrm {s}, \end{aligned}$$
(1a)
$$\begin{aligned}&-D_\mathrm {s} \nabla _z c_\mathrm {s} \cdot {{\varvec{n}}}_\mathrm {s} = \left( g + C \partial _t\left( \varPhi _\mathrm {s} - \varPhi _\mathrm {e} \right) \right) \quad \text { at }\quad \partial \varOmega _\mathrm {s}^{\mathrm {in}}, \end{aligned}$$
(1b)

where \(D_\mathrm {s}\) is the diffusion coefficient of lithium in the electrode (which may depend on the concentration itself, space variables, and temperature), G is the ratio between applied and exchange current, and C is the double layer capacitance. The exchange current g and the overpotential \(\eta \) are defined as

$$\begin{aligned} g&= c_\mathrm {s}^{\beta } \left( 1 - \frac{c_\mathrm {s}}{c_\mathrm {s}^{\max }} \right) ^{1-\beta } c_{\mathrm {e}}^{1-\beta } \left[ \exp \left( (1-\beta ) \lambda \frac{\eta }{1 + \gamma T} \right) - \exp \left( - \beta \lambda \frac{\eta }{1 + \gamma T} \right) \right] , \end{aligned}$$
(1c)
$$\begin{aligned} \eta&= \varPhi _\mathrm {s}- \varPhi _\mathrm {e} - U_{\mathrm {ocp}}\left( c_\mathrm {s}\right) , \end{aligned}$$
(1d)

where \(\beta \) is the charge-transfer coefficient of the interacalation reaction, \(c_\mathrm {s}^{\max }\) is the maximum concentration of the electrode, \(\lambda \) is the ratio between the typical and thermal potentials, \(\gamma \) is the ratio between temperature variation and reference temperature, and \(U_{\mathrm {ocp}}\) is the open circuit potential as a function of the electrode lithium concentration at the interface. All these parameters are dimensionless quantities, and more details about the non-dimensionalisation are provided in Appendix A. In addition, \(c_s\) must be periodic in \(z \in {\mathbb {R}}^3\). Even though here we consider the electrode particles to be isotropic, the same model would apply to anisotropic materials just taking the diffusion coefficients \(D_\mathrm {s}\) to be tensors.

At the mesoscale, the conservation of charge equation is posed in each electrode separately

$$\begin{aligned} \nabla _y \cdot {{\varvec{i}}}_\mathrm {s}&= -J \quad \text { in } \varOmega _\mathrm {p}, \end{aligned}$$
(2a)
$$\begin{aligned} {{\varvec{i}}}_\mathrm {s} \cdot {{\varvec{n}}}_\mathrm {p}&= - i_\mathrm {app}(y,t) \quad \text { at } \, \partial \varOmega _\mathrm {p}^{\mathrm {collector}}, \end{aligned}$$
(2b)
$$\begin{aligned} {{\varvec{i}}}_\mathrm {s} \cdot {{\varvec{n}}}_\mathrm {p}&= 0 \quad \text { at } \partial \varOmega _\mathrm {p} \setminus \partial \varOmega _\mathrm {p}^{\mathrm {collector}}, \end{aligned}$$
(2c)

with the current \({{\varvec{i}}}_\mathrm {s}\) and the exchange current per unit of volume J defined as

$$\begin{aligned} {{\varvec{i}}}_\mathrm {s}&= -{\mathcal {S}} \nabla _y \varPhi _\mathrm {s}, \end{aligned}$$
(2d)
$$\begin{aligned} J&= \frac{1}{|\varOmega |} \int _{\partial \varOmega _\mathrm {s}^{\mathrm {in}}} G \left( g + C \partial _t\left( \varPhi _\mathrm {s} - \varPhi _\mathrm {e} \right) \right) {\mathrm{d}}A_z, \end{aligned}$$
(2e)

In the previous equations, \(i_\mathrm {app}\) is the current density applied to the battery (defined to be positive during discharge), and \({\mathcal {S}}\) is the mesoscale electronic conductivity tensor in the electrode, which may vary in space and depend on temperature. The domain \(\partial \varOmega _{\mathrm {s}}^\mathrm {in}\) is the electrode–electrolyte contact surface as shown in Fig. 2. Because in this model we do not consider the current collectors, we assume that the applied current density \(i_\mathrm {app}\) is known. However, in most of the cell configurations, in order to determine \(i_\mathrm {app}\) an additional problem at the macroscale for the specific geometry of the cell. This has been addressed in the literature mostly from a numerical point of view [26,27,28]. The incorporation of the current collectors into this homogenisation analysis is an area of future work.

Notice that two instances of problems (1) and (2) are needed to describe a cell, one for each electrode. The problem in (2) is for the positive electrode, while the problem in the negative electrode is the same but with the boundary condition (2b) having the opposite sign. In each electrode the parameters can take different values, but they are constant in each electrode. Also, the microstructure, defined by \(\varOmega _s\), is different in each electrode and the concentration fields \(c_\mathrm {s}\) may have significantly different behaviours in each electrode.

The electrolyte problem can be posed across the electrodes and the separator, and it is given by

$$\begin{aligned}&\varphi _\mathrm {e} \partial _tc_\mathrm {e} - \nabla _y \cdot \left( {\mathcal {D}}_\mathrm {L} \nabla _y c_\mathrm {e} + \lambda {\mathcal {M}}_\mathrm {L} c_\mathrm {e} \nabla _y \varPhi _\mathrm {e} \right) = J \quad \text { in } \varOmega _\mathrm {cell}, \end{aligned}$$
(3a)
$$\begin{aligned}&\nabla _y \cdot {{\varvec{i}}}_\mathrm {e}= J \quad \text { in } \varOmega _\mathrm {cell}, \end{aligned}$$
(3b)
$$\begin{aligned}&- \left( {\mathcal {D}}_\mathrm {L} \nabla _y c_\mathrm {e} + \lambda {\mathcal {M}}_\mathrm {L} c_\mathrm {e} \nabla _y \varPhi _\mathrm {e} \right) \cdot {{\varvec{n}}}_\mathrm {cell} = 0 \quad \text { at } \partial \varOmega _\mathrm {cell},\end{aligned}$$
(3c)
$$\begin{aligned}&{{\varvec{i}}}_\mathrm {e}\cdot {{\varvec{n}}}_\mathrm {cell} = 0 \quad \text { at } \partial \varOmega _\mathrm {cell}, \end{aligned}$$
(3d)

with the current in the electrolyte \({{\varvec{i}}}_\mathrm {e}\) defined as

$$\begin{aligned}&{{\varvec{i}}}_\mathrm {e}= - \left( ({\mathcal {D}}_\mathrm {L} - {\mathcal {D}}_\mathrm {A}) \nabla _y c_\mathrm {e} + \lambda ({\mathcal {M}}_\mathrm {L} + {\mathcal {M}}_\mathrm {A}) c_\mathrm {e} \nabla _y \varPhi _\mathrm {e} \right) . \end{aligned}$$
(3e)

The parameters \({\mathcal {D}}_\mathrm {L}\) and \({\mathcal {D}}_\mathrm {A}\) are the mesoscale diffusivity tensors for lithium ions and anions in the electrolyte, respectively (that is, accounting for the microstructure), \({\mathcal {M}}_\mathrm {L}\) and \({\mathcal {M}}_\mathrm {A}\) are the mesoscale ion mobilities (which are tensors too), and \(\varphi _\mathrm {e}\) is the volume fraction occupied by the electrolyte. As we show during the homogenisation process, \({\mathcal {D}}_\mathrm {L}\) and \({\mathcal {D}}_\mathrm {A}\) must be multiples of each other as any anisotropy arising in the electrolyte can be caused only by the geometry, which is the same for both types of ions, and the same applies to the pair \({\mathcal {M}}_\mathrm {L}\) and \({\mathcal {M}}_\mathrm {A}\). Both diffusivities and mobilities may depend on the ion concentration in the electrolyte and temperature. All these parameters, even if they are constant in each part (electrodes and separator) like the electrolyte volume fraction, have different values in each part. In particular, J is equal to zero in the separator as there is no reaction taking place. Notice as well that the electrolyte equations only consider the concentration of lithium ions because electroneutrality has been assumed in the electrolyte, and therefore the concentration of cations and anions must be the identical.

Because this is a homogenised model, we assume that \(D_\mathrm {s}\), \({\mathcal {D}}_\mathrm {L}\) and \({\mathcal {D}}_\mathrm {A}\) are all of the same order of magnitude. Even though diffusion in the electrode is much slower than in the electrolyte, this effect has already been captured by the fact that diffusion in the electrode occurs at a much smaller length scale. In particular, \(D_\mathrm {s}\) is the same as in the microscale model in Sect. 3, while \({\mathcal {D}}_\mathrm {L} = D_\mathrm {L} {\mathcal {B}}\) and \({\mathcal {D}}_\mathrm {A} = D_\mathrm {A} {\mathcal {B}}\), where \(D_\mathrm {L}\) and \(D_\mathrm {A}\) are the microscale diffusivities and \({\mathcal {B}}\) is the tensor accounting for the geometry of the porous material as defined in Sect. 4.7.

The mesoscale heat equation is given by

$$\begin{aligned} \theta \partial _tT&= \nabla _y \cdot \left( {\mathcal {K}} \nabla _y T \right) + Q, \end{aligned}$$
(4a)

where \(\theta \) is the volumetric heat capacity and \({\mathcal {K}}\) is the thermal conductivity tensor. The heat source term Q accounts for four different heat generation mechanisms: Joule heating in the electrode, Joule heating in the electrolyte, irreversible reaction heating and reversible reaction heating. These terms can be written as

$$\begin{aligned} Q_\mathrm {s}&= - \lambda {{\varvec{i}}}_\mathrm {s} \cdot \left( {\mathcal {Q}}_\mathrm {s} \nabla _y \varPhi _\mathrm {s} \right) , \end{aligned}$$
(4b)
$$\begin{aligned} Q_\mathrm {e}&= - \lambda {{\varvec{i}}}_\mathrm {e} \cdot \left( {\mathcal {Q}}_\mathrm {e} \nabla _y \varPhi _\mathrm {e} \right) , \end{aligned}$$
(4c)
$$\begin{aligned} Q_{\mathrm {irr}}&= \frac{1}{|\varOmega |} \int _{\partial \varOmega _\mathrm {s}^{\mathrm {in}}} \lambda G g \eta {\mathrm{d}}A_z, \end{aligned}$$
(4d)
$$\begin{aligned} Q_{\mathrm {rev}}&= \frac{1}{|\varOmega |} \int _{\partial \varOmega _\mathrm {s}^{\mathrm {in}}} \lambda G g \varPi {\mathrm{d}}A_z, \end{aligned}$$
(4e)

where \({\mathcal {Q}}_\mathrm {s}\) and \({\mathcal {Q}}_\mathrm {e}\) are tensors that account for the microstructure effects in the heat generation, and \(\varPi \) is the Peltier term. The Peltier term is defined as

$$\begin{aligned} \varPi = T \frac{\partial U_\mathrm {ocp}}{\partial T}, \end{aligned}$$
(4f)

so it is temperature dependent. However, given that in many practical applications the \(U_\mathrm {ocp}\) is provided from experimental data (and so is its derivative with respect to temperature), we treat the Peltier term as a parameter function of the model.

Then, the heat generation term is defined as

$$\begin{aligned} Q = {\left\{ \begin{array}{ll} Q_\mathrm {s} + Q_\mathrm {e} + Q_{\mathrm {irr}} + Q_{\mathrm {rev}} &{}\quad \text { in }\,\, \varOmega _\mathrm {p} \text { and } \varOmega _\mathrm {n}, \\ Q_\mathrm {e} &{} \quad \text { in }\,\, \varOmega _\mathrm {sep}. \end{array}\right. } \end{aligned}$$
(4g)

The boundary conditions are not specified at this point as they depend on the problem of interest. To study a single cell, heat exchange conditions at the boundary could be used. If, instead, this equation is to be homogenised to obtain the battery heat equation, then periodic boundary conditions are required.

In a similar way, the heat equation at the battery level can be defined as

$$\begin{aligned} \theta _{\mathrm {batt}} \partial _tT&= \nabla _x \cdot \left( {\mathcal {K}}_{\mathrm {batt}} \nabla _x T \right) + Q_{\mathrm {batt}}\, \text { in }\, \varOmega _\mathrm {batt}, \end{aligned}$$
(5a)

with suitable boundary conditions, where \(\theta _{\mathrm {batt}}\) is the average volumetric heat capacity of a cell, and \({\mathcal {K}}_{\mathrm {batt}}\) is the thermal conductivity tensor of the battery. The heat source term \(Q_{\mathrm {batt}}\) accounts for the heat generation in each cell and it is defined as

$$\begin{aligned} \begin{aligned} Q_{\mathrm {batt}}&= \frac{1}{|\varOmega _{\mathrm {cell}}|} \left( \int _{\varOmega _p} \left( Q_{\mathrm {s}} + Q_{\mathrm {irr}} + Q_{\mathrm {rev}} \right) {\mathrm{d}}V_y \right. \left. + \int _{\varOmega _n} \left( Q_{\mathrm {s}} + Q_{\mathrm {irr}} + Q_{\mathrm {rev}} \right) {\mathrm{d}}V_y + \int _{\varOmega _{\mathrm {cell}}} Q_e {\mathrm{d}}V_y \right) \\&= \frac{1}{|\varOmega _{\mathrm {cell}}|} \int _{\varOmega _{\mathrm {cell}}} Q {\mathrm{d}}V_y. \end{aligned} \end{aligned}$$
(5b)

An extra term could be added to \(Q_\mathrm {batt}\) to account for the heat generated in the current collectors, but for simplicity of the homogenisation problem we have not considered it.

All the tensors that appear in the homogenised model can be calculated from material properties and microstructure, as detailed in Sect. 4.7.

3 Dimensionless microscale model

We now consider the dimensionless microscale model, from which we derive the effective model. The electrochemical model accounts for conservation of mass and charge in both the electrodes and the electrolyte, with Butler-Volmer kinetics for the intercalation reaction. The transport of lithium in the electrodes is by diffusion only, while in the electrolyte there is diffusion and migration due to the electric field. We use a dilute electrolyte model, and thus use Nernst-Planck equations, but analogous results can be found for concentrated electrolytes (see [21, 24]). We assume as well that charge transport in the electrodes follows Ohm’s law. The thermal model imposes conservation of energy in both electrode and electrolyte accounting only by diffusion effects, with a heat source term at the interface between them due to the chemical reaction. For a complete discussion of the microscale equations we refer the reader to the books [3, 29].

The variables for the microscale problem are still the concentration of lithium in the electrodes \(c_\mathrm {s}\) and lithium ions in the electrolyte \(c_\mathrm {e}\), the potentials \(\varPhi _\mathrm {s}\) and \(\varPhi _\mathrm {e}\), and the temperature T. To be rigorous, these variables are not the same as those defined in Sect. 2, which are the leading order term in the asymptotic expansions of the microscale variables. However, to keep the notation simple, we do not use any symbol to distinguish them because for the rest of the paper we will refer to the microscale variables and their asymptotic expansions. For fluxes, on the other hand, we need to distinguish the homogenised ones from the bulk ones observed at the microscale; therefore, we use tilde for the fluxes at the microscale (i.e. \(\tilde{{{\varvec{N}}}}_\mathrm {s}\), \(\tilde{{{\varvec{i}}}}_\mathrm {s}\) and \(\tilde{{{\varvec{K}}}}\)).

The microscale equations are defined in a domain \(\varOmega \), which is divided into the electrode particles \(\varOmega _\mathrm {s}\) and the electrolyte \(\varOmega _\mathrm {e}\). In the homogenisation process, we consider the electrode and electrolyte separately (except for the thermal model); therefore, we need to carefully define the notation in each domain. The thermal problem is defined as a single equation in the joint domain \(\varOmega \) because the thermal problem is the same in both electrode and electrolyte, just with different parameters. A detailed sketch of each of the subdomains and their boundaries is shown in Fig. 2. A key assumption for the homogenisation is that the microstructure repeats periodically; therefore, a finite domain that is representative of the geometry can be considered. The length scale of \(\varOmega \) is much smaller than the length scale of the porous electrode. Therefore, the ratio between length scales \(\delta \), which is a small parameter, arises in the equations and can be used in the homogenisation from microscale to mesoscale.

The details on the non-dimensionalisation of the microscale model are provided in Appendix A. The scalings are chosen so all the parameters in the dimensionless model are of \(O\left( 1 \right) \) except for \(\delta \) which is small. The dimensionless equations for the solid phase are

$$\begin{aligned}&\,\partial _tc_\mathrm {s}+ \nabla \cdot \tilde{{{\varvec{N}}}}_\mathrm {s}= 0 \quad \text { in }\,\varOmega _\mathrm {s}, \end{aligned}$$
(6a)
$$\begin{aligned}&\,\nabla \cdot \tilde{{{\varvec{i}}}}_\mathrm {s}= 0 \quad \text { in }\, \varOmega _\mathrm {s},\end{aligned}$$
(6b)
$$\begin{aligned}&\,\tilde{{{\varvec{N}}}}_\mathrm {s}\cdot {{\varvec{n}}}_\mathrm {s}= \delta G \left( g + C \partial _t\left( \varPhi _\mathrm {s}- \varPhi _\mathrm {e}\right) \right) \quad \text { at } \,\partial \varOmega _\mathrm {s}^{\mathrm {in}},\end{aligned}$$
(6c)
$$\begin{aligned}&\,\tilde{{{\varvec{i}}}}_\mathrm {s}\cdot {{\varvec{n}}}_\mathrm {s}= \delta G \left( g + C \partial _t\left( \varPhi _\mathrm {s}- \varPhi _\mathrm {e}\right) \right) \quad \text { at } \,\partial \varOmega _\mathrm {s}^{\mathrm {in}},\end{aligned}$$
(6d)
$$\begin{aligned}&\,\text {periodic} \quad \text { at } \,\partial \varOmega _\mathrm {s}^{\mathrm {out}}, \end{aligned}$$
(6e)

with

$$\begin{aligned}&\tilde{{{\varvec{N}}}}_\mathrm {s}= -\delta ^2 D_\mathrm {s}\nabla c_\mathrm {s},\end{aligned}$$
(6f)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_\mathrm {s}= - \sigma _\mathrm {s}\nabla \varPhi _\mathrm {s}, \end{aligned}$$
(6g)

and

$$\begin{aligned} g = c_\mathrm {s}^\beta \left( 1 - \frac{c_\mathrm {s}}{c_\mathrm {s}^{\max }} \right) ^{1-\beta } c_\mathrm {e}^{1-\beta } \left[ \exp \left( (1-\beta ) \lambda \frac{\eta }{1 + \gamma T} \right) - \exp \left( - \beta \lambda \frac{\eta }{1 + \gamma T} \right) \right] , \end{aligned}$$
(6h)
Fig. 2
figure 2

Schematic of the domain definition at the microscale. The microscale periodic cell \(\varOmega \) is composed of the electrolyte domain \(\varOmega _\mathrm {e}\) and the solid electrode domain \(\varOmega _\mathrm {s}\). For the homogenisation problem we define \(\partial \varOmega \) as the boundary of \(\varOmega \), which is composed by \(\partial \varOmega _\mathrm {s}^{\mathrm {out}}\) and \(\partial \varOmega _\mathrm {e}^{\mathrm {out}}\) depending on whether the boundary is at the solid or the electrolyte domain. The interface between electrolyte and solid is defined by both \(\partial \varOmega _\mathrm {e}^{\mathrm {in}}\) and \(\partial \varOmega _\mathrm {s}^{\mathrm {in}}\)

where the only new parameters are the ratio between length scales at micro- and mesoscale \(\delta \), and the microscale electronic conductivity \(\sigma _\mathrm {s}\). The latter may depend on temperature and the spacial variables in order to account for inhomogeneities in the electrode material. Diffusion of lithium in the solid is much slower than diffusion of ions in the electrolyte; therefore, we consider the limit in which the diffusivity in the solid is of \(O\left( \delta ^2 \right) \). As pointed out in [12], this particular limit allows diffusion to happen only at the microscale, as we expect in a porous electrode battery.

The equations for the electrolyte are

$$\begin{aligned}&\partial _tc_\mathrm {e}+ \nabla \cdot \tilde{{{\varvec{N}}}}_\mathrm {e}= 0 \quad \text { in } \varOmega _\mathrm {e},\end{aligned}$$
(7a)
$$\begin{aligned}&\nabla \cdot \tilde{{{\varvec{i}}}}_\mathrm {e}= 0 \quad \text { in } \varOmega _\mathrm {e},\end{aligned}$$
(7b)
$$\begin{aligned}&\tilde{{{\varvec{N}}}}_\mathrm {e}\cdot {{\varvec{n}}}_\mathrm {e}= -\delta G \left( g + C \partial _t\left( \varPhi _\mathrm {s}- \varPhi _\mathrm {e}\right) \right) \quad \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {in}},\end{aligned}$$
(7c)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_\mathrm {e}\cdot {{\varvec{n}}}_\mathrm {e}= -\delta G \left( g + C \partial _t\left( \varPhi _\mathrm {s}- \varPhi _\mathrm {e}\right) \right) \quad \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {in}},\end{aligned}$$
(7d)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {out}}, \end{aligned}$$
(7e)

with

$$\begin{aligned}&\tilde{{{\varvec{N}}}}_\mathrm {e}= -\left( D_\mathrm {L} \nabla c_\mathrm {e}+ \lambda \mu _\mathrm {L} c_\mathrm {e}\nabla \varPhi _\mathrm {e}\right) ,\end{aligned}$$
(7f)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_\mathrm {e}= -\left( (D_\mathrm {L} - D_\mathrm {A}) \nabla c_\mathrm {e}+ \lambda (\mu _\mathrm {L} + \mu _\mathrm {A}) c_\mathrm {e}\nabla \varPhi _\mathrm {e}\right) . \end{aligned}$$
(7g)

The new parameters here are \(D_\mathrm {L}\) and \(D_\mathrm {A}\) which are the diffusion coefficients of lithium ions and anions in the electrolyte, and \(\mu _\mathrm {L}\) and \(\mu _\mathrm {A}\) which are the ion mobilities. These are scalars as they are bulk parameters and do not account for the porous structure of the medium. We assume that diffusivities and mobilities of ions are homogeneous in space and do not depend on the concentration of ions or the voltage. However, they may depend on temperature. Notice the \(\delta \) factor in the fluxes (6c)–(6d) and (7c)–(7d). This scaling arises from the fact that the surface area of the particles is of \(O\left( \delta ^{-1} \right) \), and given that the total exchange current over each electrode is of \(O\left( 1 \right) \) this implies that the exchange currents (and ion fluxes) must be of \(O\left( \delta \right) \). To pose the electrolyte microscale equations we have assumed electroneutrality, which is true except in a thin layer at the electrode–electrolyte interface. The size of this layer is defined by the Debye length, which is of the order of nanometres. Therefore, even at the microscale (which is of the order of micrometers) the electroneutrality assumption is reasonable.

The dimensionless thermal model is

$$\begin{aligned}&\rho c_p \partial _tT + \nabla \cdot \tilde{{{\varvec{K}}}} = - \lambda \tilde{{{\varvec{i}}}} \cdot \nabla \varPhi \quad \text { in } \varOmega ,\end{aligned}$$
(8a)
$$\begin{aligned}&\left( \tilde{{{\varvec{K}}}}_\mathrm {s}- \tilde{{{\varvec{K}}}}_\mathrm {e}\right) \cdot {{\varvec{n}}}_\mathrm {e}= \delta \lambda G (g \eta + g \varPi ) \quad \text { in } \partial \varOmega ^{\mathrm {in}}_\mathrm {e},\end{aligned}$$
(8b)
$$\begin{aligned}&\text {periodic} \text { at } \partial \varOmega , \end{aligned}$$
(8c)

where

$$\begin{aligned} \tilde{{{\varvec{K}}}} = -k \nabla T, \end{aligned}$$
(8d)

and \(\rho \) is the density, \(c_p\) the specific heat capacity, and k the thermal conductivity. All these parameters depend on z as they are different in each material, and they can vary within each material too. We also have that \(\tilde{{{\varvec{i}}}}\) and \(\varPhi \) correspond to \(\tilde{{{\varvec{i}}}}_\mathrm {e}\) and \(\varPhi _\mathrm {e}\), or \(\tilde{{{\varvec{i}}}}_\mathrm {s}\) and \(\varPhi _\mathrm {s}\) depending on whether the point in the domain is in \(\varOmega _\mathrm {e}\) or \(\varOmega _\mathrm {s}\). We use \(\tilde{{{\varvec{K}}}}_\mathrm {s}\) to denote the flux on the electrode side of the interface and \(\tilde{{{\varvec{K}}}}_\mathrm {e}\) to denote the flux on the electrolyte side.

These microscale equations account for the classic conservation of mass, charge and energy laws, so they are physically consistent (see [3, 29] for details). However, given the complexity of the porous structure and the multiple scales involved in the problem, solving this model numerically is very challenging. The motivation to derive the homogenised mesoscale equations is to obtain a model of a lower complexity level that still captures most of the microstructure effects.

4 Derivation of the effective equations

We now proceed to homogenise the microscale equations defined in Sect. 3 in order to derive the effective equations presented in Sect. 2. To do so, we need to homogenise the equations at the microscale (i.e. porous material structure) to obtain the mesoscale equations (i.e. cell level). Then, we can assemble the mesoscale equations for a cell and homogenise them to obtain the macroscale equations (i.e. battery level).

In our problem, we have a particle microstructure of length scale \(\ell \) composing a porous material of length scale L which makes a cell. A large number of these cells are put together to make a battery of length scale NL, where N is number of cells in the battery. Then, we define the following dimensionless numbers

$$\begin{aligned} \delta = \frac{\ell }{L} \quad \text { and } \quad \epsilon = \frac{1}{N}, \end{aligned}$$
(9)

and we have that both \(\delta ,\epsilon \ll 1\). These are the small numbers that we exploit for the homogenisation. In particular, for lithium ion batteries, using as typical values those in [14, 15] and that the radius of a typical cylindrical cells is around 1 cm [30], we estimate values of \(\delta \) and \(\epsilon \) of the order of \(10^{-2}\).

Taking the limit of small \(\delta \) and \(\epsilon \), we can apply the chain rule to split the differential operators into the multiple scales. Recall, that the variable at the macroscale is x, at the mesoscale is y, and at the microscale is z. Then, when deriving the mesoscale equations, the operator acts on both the meso and microscales, so for a general function f we have

$$\begin{aligned} \nabla f\left( y,z\right) = \nabla f\left( y,\frac{y}{\delta }\right)&= \nabla _y f(y,z) + \frac{1}{\delta } \nabla _z f(y,z). \end{aligned}$$
(10a)

When we assemble the cell model and upscale it to the battery level, the macro- and mesoscales are involved, so we have

$$\begin{aligned} \nabla f\left( x,y\right) = \nabla f\left( x,\frac{x}{\epsilon }\right)&= \nabla _x f(x,y) + \frac{1}{\epsilon } \nabla _y f(x,y), \end{aligned}$$
(10b)

where the subscripts in \(\nabla \) denote with respect to which variable the operator is applied.

The homogenisation analysis presented here assumes that all the parameters except for \(\delta \) or \(\epsilon \) are of \(O\left( 1 \right) \). In practise, for the typical lithium ion chemistries, some of these parameter are big or small but this does not affect the homogenisation analysis. Having defined the different scales and the operators, we can now proceed to homogenise the equations.

4.1 Conservation of mass in the electrode

We start by homogenising the conservation of mass in the electrode, which written accounting for the multiple scales, is given by

$$\begin{aligned}&\partial _tc_\mathrm {s}+ \nabla _y \cdot \tilde{{{\varvec{N}}}}_\mathrm {s}+ \frac{1}{\delta } \nabla _z \cdot \tilde{{{\varvec{N}}}}_\mathrm {s}= 0 \quad \text { in } \varOmega _\mathrm {s}, \end{aligned}$$
(11a)
$$\begin{aligned}&\tilde{{{\varvec{N}}}}_\mathrm {s}\cdot {{\varvec{n}}}_\mathrm {s}= \delta G \left( g + C \partial _t\left( \varPhi _\mathrm {s}- \varPhi _\mathrm {e}\right) \right) \quad \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {in}},\end{aligned}$$
(11b)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {out}}, \end{aligned}$$
(11c)

with

$$\begin{aligned} \tilde{{{\varvec{N}}}}_\mathrm {s}= -\delta ^2 \quad D_\mathrm {s}\left( \nabla _y c_\mathrm {s}+ \frac{1}{\delta } \nabla _z c_\mathrm {s}\right) . \end{aligned}$$
(11d)

We expand the concentration and the flux as

$$\begin{aligned} c_\mathrm {s}&= c_{\mathrm {s},0} + \delta c_{\mathrm {s},1} + \delta ^2 c_{\mathrm {s},2} + O\left( \delta ^3 \right) ,\end{aligned}$$
(12a)
$$\begin{aligned} \tilde{{{\varvec{N}}}}_\mathrm {s}&= \delta \tilde{{{\varvec{N}}}}_{\mathrm {s},1} + \delta ^2 \tilde{{{\varvec{N}}}}_{\mathrm {s},2} + O\left( \delta ^3 \right) , \end{aligned}$$
(12b)

so now we can substitute these expansions into (11) and linearise the problem.

At leading order we find

$$\begin{aligned}&\partial _tc_{\mathrm {s},0} + \nabla _{z} \cdot \tilde{{{\varvec{N}}}}_{\mathrm {s},1} = 0 \text { in } \varOmega _\mathrm {s},\end{aligned}$$
(13a)
$$\begin{aligned}&\tilde{{{\varvec{N}}}}_{\mathrm {s},1} \cdot {{\varvec{n}}}_\mathrm {s}= G \left( g_0 + C \partial _t\left( \varPhi _{\mathrm {s},0} - \varPhi _{\mathrm {e},0} \right) \right) \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {in}},\end{aligned}$$
(13b)
$$\begin{aligned}&\text {periodic} \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {out}}, \end{aligned}$$
(13c)

with

$$\begin{aligned} \tilde{{{\varvec{N}}}}_{\mathrm {s},1} = -D_\mathrm {s}\nabla _z c_{\mathrm {s},0}, \end{aligned}$$
(13d)

where \(g_0\), \(\varPhi _{\mathrm {s},0}\) and \(\varPhi _{\mathrm {e},0}\) are the leading order terms in the expansion for small \(\delta \) of g, \(\varPhi _\mathrm {s}\) and \(\varPhi _\mathrm {e}\), respectively. Then, we conclude that the governing equation for the conservation of mass in the electrodes is

$$\begin{aligned}&\partial _tc_{\mathrm {s},0} = \nabla _{z} \cdot \left( D_\mathrm {s}\nabla _z c_{\mathrm {s},0}\right) \text { in } \varOmega _\mathrm {s},\end{aligned}$$
(14a)
$$\begin{aligned}&-D_\mathrm {s}\nabla _z c_{\mathrm {s},0} \cdot {{\varvec{n}}}_\mathrm {s}= G \left( g_0 + C \partial _t\left( \varPhi _{\mathrm {s},0} - \varPhi _{\mathrm {e},0} \right) \right) \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {in}},\end{aligned}$$
(14b)
$$\begin{aligned}&\text {periodic} \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {out}}, \end{aligned}$$
(14c)

which is the problem stated in (1). Notice that this problem is still at the microscale because the diffusion coefficient is of \(O\left( \delta ^2 \right) \) and, therefore, diffusion is so slow that can only be observed at the microscale. Because the homogenisation process does not change the diffusion coefficients for lithium in the electrode particles, the same result holds for anisotropic materials, just replacing the diffusion coefficient by a tensor.

4.2 Conservation of charge in the electrode

Next we consider the conservation of charge in the electrode, which is given by

$$\begin{aligned}&\nabla _y \cdot \tilde{{{\varvec{i}}}}_\mathrm {s}+ \frac{1}{\delta } \nabla _z \cdot \tilde{{{\varvec{i}}}}_\mathrm {s}= 0 \quad \text { in } \varOmega _\mathrm {s},\end{aligned}$$
(15a)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_\mathrm {s}\cdot {{\varvec{n}}}_\mathrm {s}= \delta G \left( g + C \partial _t\left( \varPhi _\mathrm {s}- \varPhi _\mathrm {e}\right) \right) \quad \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {in}},\end{aligned}$$
(15b)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {out}}, \end{aligned}$$
(15c)

with

$$\begin{aligned} \tilde{{{\varvec{i}}}}_\mathrm {s}&= -\,\sigma _\mathrm {s}\left( \nabla _y \varPhi _\mathrm {s}+ \frac{1}{\delta } \nabla _z \varPhi _\mathrm {s}\right) . \end{aligned}$$
(15d)

We expand the potential and the current as

$$\begin{aligned} \varPhi _\mathrm {s}&= \varPhi _{\mathrm {s},0} + \delta \varPhi _{\mathrm {s},1} + \delta ^2 \varPhi _{\mathrm {s},2} + O\left( \delta ^3 \right) ,\end{aligned}$$
(16a)
$$\begin{aligned} \tilde{{{\varvec{i}}}}_\mathrm {s}&= \delta ^{-1} \tilde{{{\varvec{i}}}}_{\mathrm {s},-1} + \tilde{{{\varvec{i}}}}_{\mathrm {s},0} + \delta \tilde{{{\varvec{i}}}}_{\mathrm {s},1} + O\left( \delta ^2 \right) . \end{aligned}$$
(16b)

Notice that the expansion for current starts at \(\delta ^{-1}\) because of the definition (15d). Substituting these expansions into (15), multiplying by a power of \(\delta \) so the leading order term is \(O\left( 1 \right) \), and expanding for small \(\delta \) we find the following problems.

4.2.1 O(1) problem

At leading order we have the problem

$$\begin{aligned}&\nabla _z \cdot \tilde{{{\varvec{i}}}}_{\mathrm {s},-1} = 0 \quad \text { in } \varOmega _\mathrm {s}, \end{aligned}$$
(17a)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_{\mathrm {s},-1} \cdot {{\varvec{n}}}_\mathrm {s}= 0 \quad \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {in}},\end{aligned}$$
(17b)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {out}}, \end{aligned}$$
(17c)

with

$$\begin{aligned} \tilde{{{\varvec{i}}}}_{\mathrm {s},-1}&= - \sigma _\mathrm {s}\nabla _z \varPhi _{\mathrm {s},0}. \end{aligned}$$
(17d)

The first step in the homogenisation process is to show that \(\varPhi _{\mathrm {s},0}\) is independent of z. Multiplying (17a) by \(\varPhi _{\mathrm {s},0}\) and integrating it over the \(\varOmega _\mathrm {s}\), we have that

$$\begin{aligned} \int _{\varOmega _\mathrm {s}} \nabla _z \cdot \left( \sigma _\mathrm {s}\nabla _z \varPhi _{\mathrm {s},0} \right) \varPhi _{\mathrm {s},0} \,{\mathrm{d}}V_z = 0, \end{aligned}$$
(18)

and applying the divergence theorem to the left hand side we find

$$\begin{aligned} -\int _{\varOmega _\mathrm {s}} \sigma _\mathrm {s}|\nabla _z \varPhi _{\mathrm {s},0}|^2\, {\mathrm{d}}V_z + \int _{\partial \varOmega _\mathrm {s}} \varPhi _{\mathrm {s},0} \sigma _\mathrm {s}\nabla _z \varPhi _{\mathrm {s},0} \cdot {{\varvec{n}}}_s \,{\mathrm{d}}A_z = 0. \end{aligned}$$
(19)

The boundary integral over \(\partial \varOmega _\mathrm {s}\) vanishes because of the conditions (17b) and (17c), and given that \(\sigma _s\) is positive, \(|\nabla _z \varPhi _{\mathrm {s},0}|^2 = 0\) over the entire domain \(\varOmega _\mathrm {s}\) which implies that \(\varPhi _{\mathrm {s},0}\) is independent of z. This also means that \(\tilde{{{\varvec{i}}}}_{\mathrm {s},-1} = {{\varvec{0}}}\) and, therefore, the expansion of the current starts at \(O\left( 1 \right) \) as one would expect.

4.2.2 \(O(\delta )\) problem

Now we consider the problem at \(O\left( \delta \right) \), which, using the results from the \(O\left( 1 \right) \) problem, is given by

$$\begin{aligned}&\nabla _z \cdot \tilde{{{\varvec{i}}}}_{\mathrm {s},0} = 0 \quad \text { in } \varOmega _\mathrm {s},\end{aligned}$$
(20a)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_{\mathrm {s},0} \cdot {{\varvec{n}}}_\mathrm {s}= 0 \quad \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {in}},\end{aligned}$$
(20b)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \,\varOmega _\mathrm {s}^{\mathrm {out}}, \end{aligned}$$
(20c)

with

$$\begin{aligned} \tilde{{{\varvec{i}}}}_{\mathrm {s},0}&= - \,\sigma _\mathrm {s}\left( \nabla _y \varPhi _{\mathrm {s},0} + \nabla _z \varPhi _{\mathrm {s},1}\right) . \end{aligned}$$
(20d)

Following the usual procedure for homogenisation (see [7] for details), we now write \(\varPhi _{\mathrm {s},1} = {{\varvec{W}}}_\mathrm {s}\cdot \nabla _y \varPhi _{\mathrm {s},0}\), where \({{\varvec{W}}}_\mathrm {s}\) depends only on z. Substitution into (20) gives the cell problem

$$\begin{aligned}&\nabla _z \cdot \left( \sigma _\mathrm {s}\left( {\mathcal {I}} + \nabla _z {{\varvec{W}}}_\mathrm {s}\right) \right) = {{\varvec{0}}} \quad \text { in } \varOmega _\mathrm {s}, \end{aligned}$$
(21a)
$$\begin{aligned}&\sigma _\mathrm {s}\left( {\mathcal {I}} + \nabla _z {{\varvec{W}}}_\mathrm {s}\right) {{\varvec{n}}}_\mathrm {s}= {{\varvec{0}}} \quad \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {in}},\end{aligned}$$
(21b)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {out}},\end{aligned}$$
(21c)
$$\begin{aligned}&\int _{\varOmega _\mathrm {s}} {{\varvec{W}}}_\mathrm {s}{\mathrm{d}}V_z = {{\varvec{0}}}, \end{aligned}$$
(21d)

which determines \({{\varvec{W}}}_\mathrm {s}\).

4.2.3 \(O(\delta ^2)\) problem

We finally address the \(O\left( \delta ^2 \right) \) problem which is given by

$$\begin{aligned}&\nabla _y \cdot \tilde{{{\varvec{i}}}}_{\mathrm {s},0} + \nabla _z \cdot \tilde{{{\varvec{i}}}}_{\mathrm {s},1} = 0 \quad \text { in } \varOmega _\mathrm {s},\end{aligned}$$
(22a)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_{\mathrm {s},1} \cdot {{\varvec{n}}}_\mathrm {s}= G \left( g_0 + C \partial _t\left( \varPhi _{\mathrm {s},0} - \varPhi _{\mathrm {e},0} \right) \right) \quad \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {in}},\end{aligned}$$
(22b)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \varOmega _\mathrm {s}^{\mathrm {out}}, \end{aligned}$$
(22c)

with

$$\begin{aligned} \tilde{{{\varvec{i}}}}_{\mathrm {s},1}&= - \sigma _\mathrm {s}\left( \nabla _y \varPhi _{\mathrm {s},1} + \nabla _z \varPhi _{\mathrm {s},2}\right) . \end{aligned}$$
(22d)

We now average (22a) over the domain \(\varOmega \) to determine the homogenised equations. Averaging the first term we find

$$\begin{aligned} \frac{1}{|\varOmega |}\int _\varOmega \nabla _y \cdot \tilde{{{\varvec{i}}}}_{\mathrm {s},0} {\mathrm{d}}V_z = \nabla _y \cdot {{\varvec{i}}}_{\mathrm {s},0}, \end{aligned}$$
(23)

where

$$\begin{aligned} {{\varvec{i}}}_{\mathrm {s},0} = -{\mathcal {S}} \nabla _y \varPhi _{\mathrm {s},0}, \end{aligned}$$
(24)

is the homogenised current, and

$$\begin{aligned} {\mathcal {S}} = \frac{1}{|\varOmega |}\int _\varOmega \sigma _\mathrm {s}\left( {\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_\mathrm {s}\right) ^\mathrm {T} \right) {\mathrm{d}}V_z, \end{aligned}$$
(25)

is the electric conductivity tensor.

We can apply the divergence theorem to the second term of (22a) jointly with the conditions (22b) and (22c) to obtain

$$\begin{aligned} \frac{1}{|\varOmega |}\int _\varOmega \nabla _z \cdot \tilde{{{\varvec{i}}}}_{\mathrm {s},1} {\mathrm{d}}V_z = \frac{1}{|\varOmega |}\int _{\partial \varOmega _\mathrm {s}^\mathrm {in}} G \left( g_0 + C \partial _t\left( \varPhi _{\mathrm {s},0} - \varPhi _{\mathrm {e},0} \right) \right) {\mathrm{d}}A_z. \end{aligned}$$
(26)

Therefore, we conclude

$$\begin{aligned} \nabla _y \cdot {{\varvec{i}}}_{\mathrm {s},0} = -J, \end{aligned}$$
(27a)

where

$$\begin{aligned} J = \frac{1}{|\varOmega |}\int _{\partial \varOmega _\mathrm {s}^\mathrm {in}} G \left( g_0 + C \partial _t\left( \varPhi _{\mathrm {s},0} - \varPhi _{\mathrm {e},0} \right) \right) {\mathrm{d}}A_z. \end{aligned}$$
(27b)

4.3 Conservation of mass and charge in the electrolyte

We now focus on the equations in the electrolyte presented in (7). Splitting the differential operators into the two different scales we find

$$\begin{aligned}&\partial _tc_\mathrm {e}+ \nabla _y \cdot \tilde{{{\varvec{N}}}}_\mathrm {e}+ \frac{1}{\delta } \nabla _z \cdot \tilde{{{\varvec{N}}}}_\mathrm {e}= 0 \text { in } \varOmega _\mathrm {e}, \end{aligned}$$
(28a)
$$\begin{aligned}&\nabla _y \cdot \tilde{{{\varvec{i}}}}_\mathrm {e}+ \frac{1}{\delta } \nabla _z \cdot \tilde{{{\varvec{i}}}}_\mathrm {e}= 0 \text { in } \varOmega _\mathrm {e},\end{aligned}$$
(28b)
$$\begin{aligned}&\tilde{{{\varvec{N}}}}_\mathrm {e}\cdot {{\varvec{n}}}_\mathrm {e}= -\,\delta G \left( g + C \partial _t\left( \varPhi _\mathrm {s}- \varPhi _\mathrm {e}\right) \right) \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {in}},\end{aligned}$$
(28c)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_\mathrm {e}\cdot {{\varvec{n}}}_\mathrm {e}= -\,\delta G \left( g + C \partial _t\left( \varPhi _\mathrm {s}- \varPhi _\mathrm {e}\right) \right) \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {in}},\end{aligned}$$
(28d)
$$\begin{aligned}&\text {periodic} \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {out}}, \end{aligned}$$
(28e)

with

$$\begin{aligned} \tilde{{{\varvec{N}}}}_\mathrm {e}&= -\left[ D_\mathrm {L} \left( \nabla _y c_\mathrm {e}+ \frac{1}{\delta } \nabla _z c_\mathrm {e}\right) + \lambda \mu _\mathrm {L} c_\mathrm {e}\left( \nabla _y \varPhi _\mathrm {e}+ \frac{1}{\delta } \nabla _z \varPhi _\mathrm {e}\right) \right] ,\end{aligned}$$
(28f)
$$\begin{aligned} \tilde{{{\varvec{i}}}}_\mathrm {e}&= -\left[ (D_\mathrm {L} - D_\mathrm {A}) \left( \nabla _y c_\mathrm {e}+ \frac{1}{\delta } \nabla _z c_\mathrm {e}\right) + \lambda (\mu _\mathrm {L} + \mu _\mathrm {A}) c_\mathrm {e}\left( \nabla _y \varPhi _\mathrm {e}+ \frac{1}{\delta } \nabla _z \varPhi _\mathrm {e}\right) \right] . \end{aligned}$$
(28g)

Now, we expand the following variables and fluxes as

$$\begin{aligned} c_\mathrm {e}&= c_{\mathrm {e},0} + \delta c_{\mathrm {e},1} + \delta ^2 c_{\mathrm {e},2} + O\left( \delta ^3 \right) ,\end{aligned}$$
(29a)
$$\begin{aligned} \varPhi _\mathrm {e}&= \varPhi _{\mathrm {e},0} + \delta \varPhi _{\mathrm {e},1} + \delta ^2 \varPhi _{\mathrm {e},2} + O\left( \delta ^3 \right) ,\end{aligned}$$
(29b)
$$\begin{aligned} \tilde{{{\varvec{N}}}}_\mathrm {e}&= \tilde{{{\varvec{N}}}}_{\mathrm {e},0} + \delta \tilde{{{\varvec{N}}}}_{\mathrm {e},1} + O\left( \delta ^2 \right) ,\end{aligned}$$
(29c)
$$\begin{aligned} \tilde{{{\varvec{i}}}}_\mathrm {e}&= \tilde{{{\varvec{i}}}}_{\mathrm {e},0} + \delta \tilde{{{\varvec{i}}}}_{\mathrm {e},1} + O\left( \delta ^2 \right) . \end{aligned}$$
(29d)

Even though the fluxes should have a term of \(O\left( \delta ^{-1} \right) \), following a similar analysis as in Sect. 4.2.1 we find that these components vanish. Then, the linearised problem yields the following equations.

4.3.1 O(1) problem

At leading order we have

$$\begin{aligned}&-\left( D_\mathrm {L} \nabla _z c_{\mathrm {e},0} + \lambda \mu _\mathrm {L} c_{\mathrm {e},0} \nabla _z \varPhi _{\mathrm {e},0} \right) = {{\varvec{0}}},\end{aligned}$$
(30a)
$$\begin{aligned}&-\left( (D_\mathrm {L} - D_\mathrm {A}) \nabla _z c_{\mathrm {e},0} + \lambda (\mu _\mathrm {L} + \mu _\mathrm {A}) c_{\mathrm {e},0} \nabla _z \varPhi _{\mathrm {e},0} \right) = {{\varvec{0}}}. \end{aligned}$$
(30b)

If the diffusion coefficients and mobilities are independent of \(c_{\mathrm {e},0}\), \(\varPhi _{\mathrm {e},0}\), and z, using a similar method as we did for \(\varPhi _{\mathrm {s},0}\) we can show that \(c_{\mathrm {e},0}\) and \(\varPhi _{\mathrm {e},0}\) do not depend on z.

4.3.2 \(O(\delta )\) problem

Using the results at \(O\left( 1 \right) \) to simplify the equations, at \(O\left( \delta \right) \) we have

$$\begin{aligned}&\nabla _z \cdot \tilde{{{\varvec{N}}}}_{\mathrm {e},0} = 0 \quad \text { in } \varOmega _\mathrm {e},\end{aligned}$$
(31a)
$$\begin{aligned}&\nabla _z \cdot \tilde{{{\varvec{i}}}}_{\mathrm {e},0} = 0 \quad \text { in } \varOmega _\mathrm {e},\end{aligned}$$
(31b)
$$\begin{aligned}&\tilde{{{\varvec{N}}}}_{\mathrm {e},0} \cdot {{\varvec{n}}}_\mathrm {e}= 0 \quad \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {in}},\end{aligned}$$
(31c)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_{\mathrm {e},0} \cdot {{\varvec{n}}}_\mathrm {e}= 0 \quad \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {in}},\end{aligned}$$
(31d)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {out}}, \end{aligned}$$
(31e)

with

$$\begin{aligned}&\tilde{{{\varvec{N}}}}_{\mathrm {e},0} = -\left( D_\mathrm {L} \left( \nabla _y c_{\mathrm {e},0} + \nabla _z c_{\mathrm {e},1} \right) + \lambda \mu _\mathrm {L} c_{\mathrm {e},0} \left( \nabla _y \varPhi _{\mathrm {e},0} + \nabla _z \varPhi _{\mathrm {e},1} \right) \right) ,\end{aligned}$$
(31f)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_{\mathrm {e},0} = -\left( (D_\mathrm {L} - D_\mathrm {A}) \left( \nabla _y c_{\mathrm {e},0} + \nabla _z c_{\mathrm {e},1} \right) + \lambda (\mu _\mathrm {L} + \mu _\mathrm {A}) c_{\mathrm {e},0} \left( \nabla _y \varPhi _{\mathrm {e},0} + \nabla _z \varPhi _{\mathrm {e},1} \right) \right) . \end{aligned}$$
(31g)

We now write \(c_{\mathrm {e},1} = {{\varvec{W}}}_\mathrm {e}\cdot \nabla _y c_{\mathrm {e},0}\) and \(\varPhi _{\mathrm {e},1} = {{\varvec{V}}}_\mathrm {e}\cdot \nabla _y \varPhi _{\mathrm {e},0}\) so, using the fact that diffusivities and mobilities do not depend on z, (31) can be rearranged into the cell problems for each one. We find that both problems are identical, so we conclude that \({{\varvec{W}}}_\mathrm {e}\equiv {{\varvec{V}}}_\mathrm {e}\). In terms of notation we use \({{\varvec{W}}}_\mathrm {e}\), which is determined by the cell problem

$$\begin{aligned}&\nabla _z \cdot \left( {\mathcal {I}} + \nabla _z {{\varvec{W}}}_\mathrm {e}\right) = {{\varvec{0}}} \quad \text { in } \varOmega _\mathrm {e},\end{aligned}$$
(32a)
$$\begin{aligned}&\left( {\mathcal {I}} + \nabla _z {{\varvec{W}}}_\mathrm {e}\right) {{\varvec{n}}}_\mathrm {e}= {{\varvec{0}}} \quad \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {in}},\end{aligned}$$
(32b)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {out}},\end{aligned}$$
(32c)
$$\begin{aligned}&\int _{\varOmega _\mathrm {e}} {{\varvec{W}}}_\mathrm {e}{\mathrm{d}}V_z = {{\varvec{0}}}. \end{aligned}$$
(32d)

4.3.3 \(O(\delta ^2)\) problem

Now we can consider the \(O\left( \delta ^2 \right) \) to determine the homogenised equations. The equations read

$$\begin{aligned}&\partial _tc_{\mathrm {e},0} + \nabla _y \cdot \tilde{{{\varvec{N}}}}_{\mathrm {e},0} + \nabla _z \cdot \tilde{{{\varvec{N}}}}_{\mathrm {e},1} = 0 \text { in } \varOmega _\mathrm {e},\end{aligned}$$
(33a)
$$\begin{aligned}&\nabla _y \cdot \tilde{{{\varvec{i}}}}_{\mathrm {e},0} + \nabla _z \cdot \tilde{{{\varvec{i}}}}_{\mathrm {e},1} = 0 \text { in } \varOmega _\mathrm {e},\end{aligned}$$
(33b)
$$\begin{aligned}&\tilde{{{\varvec{N}}}}_{\mathrm {e},1} \cdot {{\varvec{n}}}_\mathrm {e}= - G \left( g + C \partial _t\left( \varPhi _\mathrm {s}- \varPhi _\mathrm {e}\right) \right) \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {in}},\end{aligned}$$
(33c)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_{\mathrm {e},1} \cdot {{\varvec{n}}}_\mathrm {e}= - G \left( g + C \partial _t\left( \varPhi _\mathrm {s}- \varPhi _\mathrm {e}\right) \right) \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {in}},\end{aligned}$$
(33d)
$$\begin{aligned}&\text {periodic} \text { at } \partial \varOmega _\mathrm {e}^{\mathrm {out}}, \end{aligned}$$
(33e)

with

$$\begin{aligned}&\tilde{{{\varvec{N}}}}_{\mathrm {e},1} = -\left( D_\mathrm {L} \left( \nabla _y c_{\mathrm {e},1} + \nabla _z c_{\mathrm {e},2} \right) \right. \nonumber \\&\quad \quad \quad \quad \left. + \lambda \mu _\mathrm {L} \left( c_{\mathrm {e},0} \left( \nabla _y \varPhi _{\mathrm {e},1} + \nabla _z \varPhi _{\mathrm {e},2} \right) + c_{\mathrm {e},1} \left( \nabla _y \varPhi _{\mathrm {e},0} + \nabla _z \varPhi _{\mathrm {e},1} \right) \right) \right) , \end{aligned}$$
(33f)
$$\begin{aligned}&\tilde{{{\varvec{i}}}}_{\mathrm {e},1} = -\left( (D_\mathrm {L} - D_\mathrm {A}) \left( \nabla _y c_{\mathrm {e},0} + \nabla _z c_{\mathrm {e},1} \right) \right. \nonumber \\&\quad \quad \quad \quad \left. + \lambda (\mu _\mathrm {L} + \mu _\mathrm {A}) \left( c_{\mathrm {e},0} \left( \nabla _y \varPhi _{\mathrm {e},1} + \nabla _z \varPhi _{\mathrm {e},2} \right) + c_{\mathrm {e},1} \left( \nabla _y \varPhi _{\mathrm {e},0} + \nabla _z \varPhi _{\mathrm {e},1} \right) \right) \right) . \end{aligned}$$
(33g)

We average (33a) and (33b) over \(\varOmega \) and use the divergence theorem with conditions (33c)–(33e) to obtain the homogenised equations

$$\begin{aligned}&\varphi _\mathrm {e}\partial _tc_{\mathrm {e},0} + \nabla _y \cdot {{\varvec{N}}}_{\mathrm {e},0} = J, \end{aligned}$$
(34a)
$$\begin{aligned}&\nabla _y \cdot {{\varvec{i}}}_{\mathrm {e},0} = J, \end{aligned}$$
(34b)

where

$$\begin{aligned}&{{\varvec{N}}}_{\mathrm {e},0} = -\left( D_{\mathrm {L}} {\mathcal {B}} \nabla _y c_{\mathrm {e},0} + \lambda \mu _L {\mathcal {B}} c_{\mathrm {e},0} \nabla _y \varPhi _{\mathrm {e},0} \right) ,\end{aligned}$$
(34c)
$$\begin{aligned}&{{\varvec{i}}}_{\mathrm {e},0} = -\left( (D_{\mathrm {L}} - D_{\mathrm {A}}) {\mathcal {B}} \nabla _y c_{\mathrm {e},0} + \lambda (\mu _L + \mu _{\mathrm {A}}) {\mathcal {B}} c_{\mathrm {e},0} \nabla _y \varPhi _{\mathrm {e},0} \right) ,\end{aligned}$$
(34d)
$$\begin{aligned}&{\mathcal {B}} = \frac{1}{|\varOmega |}\int _\varOmega \left( {\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_\mathrm {e}\right) ^\mathrm {T} \right) {\mathrm{d}}V_z. \end{aligned}$$
(34e)

4.4 Conservation of heat at the mesoscale

The temperature equation, splitting the microscale from larger scales, is given by

$$\begin{aligned}&\rho c_p \partial _tT + \nabla _{y} \cdot \tilde{{{\varvec{K}}}} + \frac{1}{\delta } \nabla _z \cdot \tilde{{{\varvec{K}}}} = - \lambda \tilde{{{\varvec{i}}}} \cdot \left( \nabla _y \varPhi + \frac{1}{\delta } \nabla _z \varPhi \right) \quad \text { in } \varOmega ,\end{aligned}$$
(35a)
$$\begin{aligned}&\left( \tilde{{{\varvec{K}}}}_\mathrm {s}- \tilde{{{\varvec{K}}}}_\mathrm {e}\right) \cdot {{\varvec{n}}}_\mathrm {e}= \delta \lambda G (g \eta + g \varPi ) \quad \text { in } \partial \varOmega ^{\mathrm {in}}_\mathrm {e},\end{aligned}$$
(35b)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \varOmega , \end{aligned}$$
(35c)

where

$$\begin{aligned} \tilde{{{\varvec{K}}}} = -k \left( \nabla _{y} T + \frac{1}{\delta } \nabla _z T \right) . \end{aligned}$$
(35d)

Now we expand the following quantities as

$$\begin{aligned} T&= T_0 + \delta T_1 + \delta ^2 T_2 + O\left( \delta ^3 \right) ,\end{aligned}$$
(36a)
$$\begin{aligned} \varPhi&= \varPhi _0 + \delta \varPhi _1 + \delta ^2 \varPhi _2 + O\left( \delta ^3 \right) ,\end{aligned}$$
(36b)
$$\begin{aligned} \tilde{{{\varvec{K}}}}&= \tilde{{{\varvec{K}}}}_0 + \delta \tilde{{{\varvec{K}}}}_1 + O\left( \delta ^2 \right) ,\end{aligned}$$
(36c)
$$\begin{aligned} \tilde{{{\varvec{i}}}}&= \tilde{{{\varvec{i}}}}_0 + \delta \tilde{{{\varvec{i}}}}_1 + O\left( \delta ^2 \right) , \end{aligned}$$
(36d)

where again it can be shown that the \(O\left( \delta ^{-1} \right) \) contributions in the fluxes vanish.

4.4.1 O(1) problem

At leading order the problem reads

$$\begin{aligned} -k \nabla _z T_0 = {{\varvec{0}}}. \end{aligned}$$
(37)

Therefore, following a similar argument to the one for \(\varPhi _{\mathrm {s},0}\) we conclude that \(T_0\) is independent of z.

4.4.2 \(O(\delta )\) problem

We now consider the \(O\left( \delta \right) \) problem, and using that neither \(T_0\) nor \(\varPhi _0\) depend on z, we have

$$\begin{aligned}&\nabla _z \cdot \tilde{{{\varvec{K}}}}_0 = 0 \quad \text { in } \varOmega ,\end{aligned}$$
(38a)
$$\begin{aligned}&\left( \tilde{{{\varvec{K}}}}_{\mathrm {s},0} - \tilde{{{\varvec{K}}}}_{\mathrm {e},0} \right) \cdot {{\varvec{n}}}_\mathrm {e}= 0 \quad \text { in } \partial \varOmega ^{\mathrm {in}}_\mathrm {e},\end{aligned}$$
(38b)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \varOmega , \end{aligned}$$
(38c)

where

$$\begin{aligned} \tilde{{{\varvec{K}}}}_0 = -k \left( \nabla _y T_0 + \nabla _z T_1 \right) . \end{aligned}$$
(38d)

Defining \(T_1 = {{\varvec{W}}}_T \cdot \nabla _y T_0\) we obtain the cell problem for the thermal problem

$$\begin{aligned}&\nabla _z \cdot \left( k \left( {\mathcal {I}} + \nabla _z {{\varvec{W}}}_T \right) \right) = {{\varvec{0}}} \text { in } \varOmega ,\end{aligned}$$
(39a)
$$\begin{aligned}&\text {periodic} \text { at } \partial \varOmega ,\end{aligned}$$
(39b)
$$\begin{aligned}&\int _\varOmega {{\varvec{W}}}_T {\mathrm{d}}V_z = {{\varvec{0}}}, \end{aligned}$$
(39c)

and notice that we do not need to account for the condition at the internal boundary as (39a) takes care of it.

4.4.3 \(O(\delta ^2)\) problem

We finally consider the \(O\left( \delta ^2 \right) \) problem

$$\begin{aligned}&\rho c_p \partial _tT_0 + \nabla _y \cdot \tilde{{{\varvec{K}}}}_0 + \nabla _z \cdot \tilde{{{\varvec{K}}}}_1 = - \lambda \tilde{{{\varvec{i}}}}_0 \cdot \left( \nabla _y \varPhi _0 + \nabla _z \varPhi _1 \right) \text { in } \varOmega ,\end{aligned}$$
(40a)
$$\begin{aligned}&\left( \tilde{{{\varvec{K}}}}_{\mathrm {s},1} - \tilde{{{\varvec{K}}}}_{\mathrm {e},1} \right) \cdot {{\varvec{n}}}_\mathrm {e}= \lambda G (g_0 \eta _0 + g_0 \varPi _0) \text { in } \partial \varOmega ^{\mathrm {in}}_\mathrm {e},\end{aligned}$$
(40b)
$$\begin{aligned}&\text {periodic} \text { at } \partial \varOmega , \end{aligned}$$
(40c)

where

$$\begin{aligned} \tilde{{{\varvec{K}}}} = -\,k \left( \nabla _y T_1 + \nabla _z T_2 \right) . \end{aligned}$$
(40d)

We can now average (40a) over \(\varOmega \) to obtain the homogenised equation, doing the same type of manipulations as detailed in Sect. 4.2, The integral can be split up into one integral over \(\varOmega _{s}\) and one over \(\varOmega _{e}\), so it can be written as

$$\begin{aligned} \begin{aligned}&\int _{\varOmega } {{\varvec{i}}}_{0} \cdot (\nabla _{y}\varPhi _{0}+\nabla _{z}\varPhi _{1}) {\mathrm{d}}V_z = \int _{\varOmega _\mathrm {s}} {{\varvec{i}}}_{\mathrm {s},0} \cdot (\nabla _{y}\varPhi _{\mathrm {s},0} + \nabla _{z}\varPhi _{\mathrm {s},1}) {\mathrm{d}}V_z + \int _{\varOmega _\mathrm {e}} {{\varvec{i}}}_{\mathrm {e},0} \cdot (\nabla _{y}\varPhi _{\mathrm {e},0} + \nabla _{z}\varPhi _{\mathrm {e},1}) {\mathrm{d}}V_z \\&= \int _{\varOmega _\mathrm {s}} {{\varvec{i}}}_{\mathrm {s},0} \cdot \left( ({\mathcal {I}} + \nabla _{z} {{\varvec{W}}}_\mathrm {s}) \nabla _{y}\varPhi _{\mathrm {s},0} \right) {\mathrm{d}}V_z \quad + \int _{\varOmega _\mathrm {e}}{{\varvec{i}}}_{\mathrm {e},0} \cdot \left( ({\mathcal {I}} + \nabla _{z} {{\varvec{W}}}_\mathrm {e}) \nabla _{y} \varPhi _{\mathrm {e},0} \right) {\mathrm{d}}V_z. \end{aligned} \end{aligned}$$
(41)

The homogenised equation reads

$$\begin{aligned} \theta \partial _tT_0 = \nabla _y \cdot \left( {\mathcal {K}} \nabla _y T_0 \right) + Q_\mathrm {s}+ Q_\mathrm {e}+ Q_{\mathrm {irr}} + Q_{\mathrm {rev}}, \end{aligned}$$
(42a)

where

$$\begin{aligned} Q_{\mathrm {s}}&= -\, \lambda {{\varvec{i}}}_{\mathrm {s},0} \cdot \left( {\mathcal {Q}}_\mathrm {s}\nabla _y \varPhi _{\mathrm {s},0} \right) , \end{aligned}$$
(42b)
$$\begin{aligned} Q_{\mathrm {e}}&= -\, \lambda {{\varvec{i}}}_{\mathrm {e},0} \cdot \left( {\mathcal {Q}}_\mathrm {e}\nabla _y \varPhi _{\mathrm {e},0} \right) ,\end{aligned}$$
(42c)
$$\begin{aligned} Q_{\mathrm {irr}}&= \frac{1}{|\varOmega |} \int _{\partial \varOmega _\mathrm {e}^\mathrm {in}} \lambda G g_0 \eta _0\, {\mathrm{d}}A_z,\end{aligned}$$
(42d)
$$\begin{aligned} Q_{\mathrm {rev}}&= \frac{1}{|\varOmega |} \int _{\partial \varOmega _\mathrm {e}^\mathrm {in}} \lambda G g_0 \varPi _0\, {\mathrm{d}}A_z, \end{aligned}$$
(42e)

and

$$\begin{aligned} \theta&= \frac{1}{|\varOmega |} \int _\varOmega \rho c_p\, {\mathrm{d}}V_z,\end{aligned}$$
(42f)
$$\begin{aligned} {\mathcal {K}}&= \frac{1}{|\varOmega |}\int _\varOmega k \left( {\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_T\right) ^\mathrm {T} \right) {\mathrm{d}}V_z,\end{aligned}$$
(42g)
$$\begin{aligned} {\mathcal {Q}}_\mathrm {s}&= \frac{\left( {\mathcal {S}}^{-1}\right) ^\mathrm {T}}{|\varOmega |} \int _\varOmega \sigma _\mathrm {s}\left( {\mathcal {I}} + \nabla _z {{\varvec{W}}}_\mathrm {s}\right) \left( {\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_\mathrm {s}\right) ^\mathrm {T} \right) {\mathrm{d}}V_z, \end{aligned}$$
(42h)
$$\begin{aligned} {\mathcal {Q}}_\mathrm {e}&= \frac{\left( {\mathcal {B}}^{-1}\right) ^\mathrm {T}}{|\varOmega |} \int _\varOmega \left( {\mathcal {I}} + \nabla _z {{\varvec{W}}}_\mathrm {e}\right) \left( {\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_\mathrm {e}\right) ^\mathrm {T} \right) {\mathrm{d}}V_z. \end{aligned}$$
(42i)

4.5 Conservation of heat at the macroscale

We finally homogenise the heat equation at the mesoscale (cell level) to obtain the heat equation at the macroscale (battery level). The mesoscale equation comes from assembling three instances of (42) to account for the electrodes and the separator, or even consider a more complicated structure to account for double coated electrodes as done in [25]. However, given that temperature and heat flux must be continuous across the interface between parts we can write it as a single problem with space varying parameters. Then, we can split the scales as in (10b) and write the mesoscale equation for temperature as

$$\begin{aligned}&\theta \partial _tT + \nabla _{x} \cdot {{\varvec{K}}} + \frac{1}{\epsilon } \nabla _{y} \cdot {{\varvec{K}}} = Q \text { in } \varOmega _\mathrm {cell},\end{aligned}$$
(43a)
$$\begin{aligned}&\text {periodic} \text { at } \partial \varOmega _\mathrm {cell}, \end{aligned}$$
(43b)

where

$$\begin{aligned} {{\varvec{K}}} = -{\mathcal {K}} \left( \nabla _{x} T + \frac{1}{\epsilon } \nabla _y T \right) , \end{aligned}$$
(43c)

and Q is the heat generation term defined as

$$\begin{aligned} Q = {\left\{ \begin{array}{ll} Q_\mathrm {s} + Q_\mathrm {e} + Q_{\mathrm {irr}} + Q_{\mathrm {rev}} &{} \text { in } \varOmega _\mathrm {p} \text { and } \varOmega _\mathrm {n}, \\ Q_\mathrm {e} &{} \text { in } \varOmega _\mathrm {sep}. \end{array}\right. } \end{aligned}$$
(43d)

Notice that the gradients appearing in Q are at the mesoscale only, so when we expand Q in powers of \(\epsilon \) no term of \(O\left( \epsilon ^{-1} \right) \) arises.

We now expand the variables as

$$\begin{aligned} T&= T_0 + \epsilon T_1 + \epsilon ^2 T_2 + O\left( \epsilon ^3 \right) ,\end{aligned}$$
(44a)
$$\begin{aligned} {{\varvec{K}}}&= {{\varvec{K}}}_0 + \epsilon {{\varvec{K}}}_1 + O\left( \epsilon ^2 \right) , \end{aligned}$$
(44b)

substitute them into (43) and linearise.

4.5.1 O(1) problem

At leading order we have

$$\begin{aligned} -{\mathcal {K}} \nabla _y T_0 = {{\varvec{K}}}_{0}, \end{aligned}$$
(45)

therefore, by the same argument as in Sect. 4.4, we conclude that \(T_0\) does not depend on y.

4.5.2 \(O(\epsilon )\) problem

Using the results found at leading order, the \(O\left( \epsilon \right) \) problem reads

$$\begin{aligned}&\nabla _{y} \cdot {{\varvec{K}}}_0 = 0 \quad \text { in } \varOmega _\mathrm {cell}, \end{aligned}$$
(46a)
$$\begin{aligned}&\text {periodic} \quad \text { at } \partial \varOmega _\mathrm {cell}, \end{aligned}$$
(46b)

where

$$\begin{aligned} {{\varvec{K}}}_0 = -{\mathcal {K}} \left( \nabla _{x} T_0 + \nabla _y T_1 \right) . \end{aligned}$$
(46c)

To obtain the cell problem, we make the substitution \(T_1 = {{\varvec{W}}}_\mathrm {cell} \cdot \nabla _x T_0\) finding

$$\begin{aligned}&\nabla _y \cdot \left( {\mathcal {K}} \left( {\mathcal {I}} + \nabla _y {{\varvec{W}}}_\mathrm {cell} \right) \right) = {{\varvec{0}}} \text { in } \varOmega _\mathrm {cell},\end{aligned}$$
(47a)
$$\begin{aligned}&\text {periodic} \text { at } \partial \varOmega _\mathrm {cell},\end{aligned}$$
(47b)
$$\begin{aligned}&\int _{\varOmega _\mathrm {cell}} {{\varvec{W}}}_\mathrm {cell} {\mathrm{d}}V_y = {{\varvec{0}}}. \end{aligned}$$
(47c)

4.5.3 \(O(\epsilon ^2)\) problem

Finally, we consider the \(O\left( \epsilon ^2 \right) \) problem given by

$$\begin{aligned}&\theta \partial _tT_0 + \nabla _{x} \cdot {{\varvec{K}}}_0 + \nabla _{y} \cdot {{\varvec{K}}}_1 = Q_0 \text { in } \varOmega _\mathrm {cell},\end{aligned}$$
(48a)
$$\begin{aligned}&\text {periodic} \text { at } \partial \varOmega _\mathrm {cell}, \end{aligned}$$
(48b)

where

$$\begin{aligned} {{\varvec{K}}}_1 = -{\mathcal {K}} \left( \nabla _{x} T_1 + \nabla _y T_2 \right) . \end{aligned}$$
(48c)

Averaging over \(\varOmega _\mathrm {cell}\) we obtain the homogenised equation

$$\begin{aligned} \theta _\mathrm {batt} \partial _tT_0 = \nabla _x \cdot \left( {\mathcal {K}}_\mathrm {batt} \nabla _x T_0 \right) + Q_\mathrm {batt}, \end{aligned}$$
(49a)

where

$$\begin{aligned} \theta _\mathrm {batt}&= \frac{1}{|\varOmega _\mathrm {cell}|} \int _{\varOmega _\mathrm {cell}} \theta {\mathrm{d}}\, V_y,\end{aligned}$$
(49b)
$$\begin{aligned} {\mathcal {K}}_\mathrm {batt}&= \frac{1}{|\varOmega _\mathrm {cell}|} \int _{\varOmega _\mathrm {cell}} {\mathcal {K}} \left( {\mathcal {I}} + \left( \nabla _y {{\varvec{W}}}_\mathrm {cell}\right) ^\mathrm {T} \right) {\mathrm{d}}V_y,\end{aligned}$$
(49c)
$$\begin{aligned} Q_\mathrm {batt}&= \frac{1}{|\varOmega _\mathrm {cell}|} \int _{\varOmega _\mathrm {cell}} Q {\mathrm{d}}\, V_y. \end{aligned}$$
(49d)

4.6 Discussion of the effective equations

We observe a few key differences between the microscale model presented in Sect. 3 and the homogenised one derived in this section. The first one is that, even though the microscale equations assumed isotropic materials and therefore all the transport properties were defined as scalars, in the homogenised equations the same transport properties become tensors and thus allow for anisotropy. This anisotropy arises from the porous structure. In some cases, such as with the electronic conductivity \({\mathcal {S}}\), the anisotropy is caused by both geometric effects and the inhomogeneities in the material properties. However, in the electrolyte properties, the anisotropy arises purely from the geometry, which allows us to define the tensor \({\mathcal {B}}\) which accounts for the anisotropy of all electrolyte properties. Therefore, all the electrolyte transport tensors at the mesoscale must be multiples of each other.

Another difference is the appearance of source terms in the electrochemical equations. These source terms capture the chemical reaction that, at the microscale, occurs at the boundary between electrode and electrolyte. However, due to the homogenisation process, this boundary is lost at the mesoscale but the reactions are captured by the source term. The same happens with the reaction contribution to the heat source terms \(Q_\mathrm {irr}\) and \(Q_\mathrm {rev}\).

In the thermal equations, both at mesoscale and macroscale, we observe that what was the \(\rho c_p\) term in the microscale equation became the averaged volumetric heat capacities \(\theta \) and \(\theta _\mathrm {batt}\), which can no longer be split into the averaged density and the averaged specific heat capacity. We also notice that the heat generation terms at the mesoscale have a similar structure to those at the microscale, but with the introduction of two tensors \({\mathcal {Q}}_\mathrm {s}\) and \({\mathcal {Q}}_\mathrm {e}\) that account for microstructure effects in Joule heating.

The macroscale thermal equation holds for any periodic cell structure, however notice that normally we encounter layered materials. In that case, as it is well known from the literature [7], we find that the tensor \({\mathcal {K}}_\mathrm {batt}\) is a diagonal tensor, and the diagonal values are the arithmetic average of the conductivities of each layer in the directions parallel to the layers, and the harmonic average of the conductivities in the direction perpendicular to the layers. Therefore, asymptotic homogenisation provides a rigorous proof to a result that is commonly used in the literature (e.g. [19]).

For the sake of clarity, when performing the homogenisation we have not explicitly accounted for the temperature dependence on the parameters, even though the parameters are allowed to depend on temperature. Given that we have determined that, at leading order, temperature does not depend on z, the analysis for temperature-dependent parameters follows exactly the same way with some extra parameters arising from the derivatives of the parameters with respect to temperature.

4.7 Properties of the tensors

In the previous sections, we have derived tensors that account for the material properties in the homogenised problem. These tensors have certain properties that show the type of behaviour to expect from the homogenised equations.

4.7.1 Symmetry

We start showing that the tensors are symmetric, which can be used to simplify some of the expressions. We show it for \({\mathcal {S}}\) because it has the most complicated definition in our model, and any other tensor that we have defined can be thought to be, from a mathematical point of view, a particular case of \({\mathcal {S}}\). The method used here follows closely the one used in [22]. For simplicity in the notation, we define \(W_\mathrm {s}^{(i)}\) to be the i-th component of the vector \({{\varvec{W}}}_\mathrm {s}\) and \({{\varvec{e}}}_i\) is the basis vector in the i-th direction. Then, we define the integral

$$\begin{aligned} \int _{\varOmega _\mathrm {s}} \nabla _z \cdot \left( \sigma _\mathrm {s}W_\mathrm {s}^{(i)} \left( {{\varvec{e}}}_j +\nabla _z W_\mathrm {s}^{(j)} \right) - \sigma _\mathrm {s}W_\mathrm {s}^{(j)} \left( {{\varvec{e}}}_i +\nabla _z W_\mathrm {s}^{(i)} \right) \right) {\mathrm{d}}V_z = 0, \end{aligned}$$
(50)

which we can show to be zero using the divergence theorem and the boundary conditions in (21). Expanding the divergence in the integral, and using (21a) to eliminate some of the terms, we find

$$\begin{aligned} \int _{\varOmega _\mathrm {s}} \sigma _\mathrm {s}\nabla _z W_\mathrm {s}^{(i)} {{\varvec{e}}}_j - \sigma _\mathrm {s}\nabla _z W_\mathrm {s}^{(j)} {{\varvec{e}}}_i {\mathrm{d}}V_z = 0, \end{aligned}$$
(51)

and therefore we have

$$\begin{aligned} \int _{\varOmega _\mathrm {s}} \sigma _\mathrm {s}\frac{\partial W_\mathrm {s}^{(i)}}{\partial z_j} {\mathrm{d}}V_z = \int _{\varOmega _\mathrm {s}} \sigma _\mathrm {s}\frac{\partial W_\mathrm {s}^{(j)}}{\partial z_i} {\mathrm{d}}V_z, \end{aligned}$$
(52)

from which we deduce that \({\mathcal {S}}\) is symmetric. Using similar procedures, we can show that the other tensors we defined (\({\mathcal {B}}\), \({\mathcal {K}}\), \({\mathcal {K}}_\mathrm {batt}\), \({\mathcal {M}}_\mathrm {s}\) and \({\mathcal {M}}_\mathrm {e}\)) are symmetric as well.

4.7.2 Reduction of the transport properties

We now want to show that the presence of the microstructure reduces the transport properties. Mathematically, this is equivalent to showing that the diagonal elements of the tensor of a given transport phenomenon are smaller than the average of the microscale value over the cell for that same property. By considering the integral

$$\begin{aligned} \int _{\varOmega _\mathrm {s}} \nabla _z \cdot \left( \sigma _\mathrm {s}\left( {{\varvec{e}}}_i + \nabla _z W_\mathrm {s}^{(i)} \right) \right) {\mathrm{d}}V_z = 0, \end{aligned}$$
(53)

and manipulating it in a similar way it can be shown that the diagonal entries of the tensor satisfy

$$\begin{aligned} {\mathcal {S}}_{ii} < \frac{1}{|\varOmega _\mathrm {s}|} \int _{\varOmega _\mathrm {s}} \sigma _\mathrm {s}{\mathrm{d}}\, V_z. \end{aligned}$$
(54)

4.7.3 Analytical representation of the tensors

Exploiting the symmetry of the tensors we can simplify the way we defined some of them. Thus, we present the simpler expressions here under the same section for convenience of the reader.

The electric conductivity tensor in the electrode is given by

$$\begin{aligned} {\mathcal {S}} = \frac{1}{|\varOmega |}\int _\varOmega \sigma _\mathrm {s}\left( {\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_\mathrm {s}\right) ^\mathrm {T} \right) {\mathrm{d}}V_z, \end{aligned}$$
(55)

where the cell variable \({{\varvec{W}}}_\mathrm {s}\) solves (21).

The geometry of the electrolyte is captured by the tensor \({\mathcal {B}}\) which can then be used to define the transport properties in the electrolyte. The tensor is defined by

$$\begin{aligned} {\mathcal {B}} = \frac{1}{|\varOmega |}\int _\varOmega \left( {\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_\mathrm {e}\right) ^\mathrm {T} \right) {\mathrm{d}}V_z, \end{aligned}$$
(56)

where the cell variable \({{\varvec{W}}}_\mathrm {e}\) solves (32). Using \({\mathcal {B}}\) we can define the tensors for diffusivities and mobilities used in (2) as

$$\begin{aligned} {\mathcal {D}}_\mathrm {L} = D_\mathrm {L} {\mathcal {B}},\quad {\mathcal {D}}_\mathrm {A} = D_\mathrm {A} {\mathcal {B}},\quad {\mathcal {M}}_\mathrm {L} = \mu _\mathrm {L} {\mathcal {B}},\quad {\mathcal {M}}_\mathrm {A} = \mu _\mathrm {A} {\mathcal {B}}. \end{aligned}$$
(57)

The thermal conductivity of the cell is given by

$$\begin{aligned} {\mathcal {K}} = \frac{1}{|\varOmega |}\int _\varOmega k \left( {\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_T\right) ^\mathrm {T} \right) {\mathrm{d}}V_z, \end{aligned}$$
(58)

where the cell variable \({{\varvec{W}}}_T\) solves (39), and the Joule heating tensors in the electrode and the electrolyte are defined as

$$\begin{aligned} {\mathcal {Q}}_\mathrm {s}&= \frac{{\mathcal {S}}^{-1}}{|\varOmega |} \int _\varOmega \sigma _\mathrm {s}\Big ({\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_\mathrm {s}\right) \Big ) \left( {\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_\mathrm {s}\right) ^\mathrm {T} \right) {\mathrm{d}}V_z, \end{aligned}$$
(59)
$$\begin{aligned} {\mathcal {Q}}_\mathrm {e}&= \frac{{\mathcal {B}}^{-1}}{|\varOmega |} \int _\varOmega \Big ({\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_\mathrm {e}\right) \Big ) \left( {\mathcal {I}} + \left( \nabla _z {{\varvec{W}}}_\mathrm {e}\right) ^\mathrm {T} \right) {\mathrm{d}}V_z. \end{aligned}$$
(60)

Finally, the thermal conductivity of the battery is given by

$$\begin{aligned} {\mathcal {K}}_\mathrm {batt} = \frac{1}{|\varOmega |_\mathrm {cell}}\int _{\varOmega _\mathrm {cell}} {\mathcal {K}} \left( {\mathcal {I}} + \left( \nabla _y {{\varvec{W}}}_\mathrm {cell}\right) ^\mathrm {T} \right) {\mathrm{d}}V_y, \end{aligned}$$
(61)

where the cell variable \({{\varvec{W}}}_\mathrm {cell}\) solves (47).

5 Comparison with the DFN model

The model presented here can be regarded as a generalised version of the well-known DFN model [1]. In this section, we compare both models and point out the main differences. The DFN model does not include the battery level nor the thermal model, so the comparison will be only at the micro- and mesoscale and at constant temperature. Another difference betwen the model presented here and the DFN model is that the first considers a capacitance term in the exchange current between the electrode and the electrolyte; thus, we can reduce to the DFN model by setting \(C = 0\).

5.1 Comparison at the microscale

One of the key assumptions of the DFN model is that the electrode particles are spherical and that there is a single representative particle at each point of the mesoscale electrode (i.e. one particle in each homogenisation cell). With these assumptions, our model (1) reduces to

$$\begin{aligned}&\partial _tc_\mathrm {s}= \frac{1}{r^{2}}\frac{\partial }{\partial r}\left( r^{2}D_\mathrm {s}\frac{\partial c_\mathrm {s}}{\partial r}\right) \text { in } 0< r < R, \end{aligned}$$
(62a)
$$\begin{aligned}&- D_\mathrm {s}\frac{\partial c_ \mathrm {s}}{\partial r} = G g \text { at } r = R,\end{aligned}$$
(62b)
$$\begin{aligned}&\frac{\partial c_\mathrm {s}}{\partial r} = 0 \text { at } r = 0, \end{aligned}$$
(62c)

which is the same as in the DFN model. Here, r is the microscale space variable and takes the role of z in our model.

5.2 Comparison at the mesoscale

Note that the parameters in the DFN model are all scalar, whereas in the model presented in this paper some parameters are tensors, which can take into account anisotropy. However, it can be reduced to the isotropic case (as in the DFN model) by taking the tensors to be multiples of the identity tensor, as then they can be replaced by a scalar.

Another main difference is the difference in dilute and concentrated theory which will necessarily lead to different expressions. However, the analysis presented here can be applied to the concentrated electrolyte equations (see [21, 24] for details) to obtain analogous results.

5.2.1 Conservation of charge in the electrode

The conservation of charge for the solid in our model (2) reads

$$\begin{aligned} \nabla \cdot ({\mathcal {S}} \nabla \varPhi _\mathrm {s}) = J. \end{aligned}$$
(63)

Due to the spherical symmetry of the particles assumed in the DFN model, we can write the exchange current as

$$\begin{aligned} J = a G \left( g + C \partial _t(\varPhi _\mathrm {s}- \varPhi _\mathrm {e}) \right) , \end{aligned}$$
(64)

where a is the surface area of the particles per unit of volume. To reduce to the DFN model, we need to neglect the capacitance of the double layer and the anistropy of the material.

Notice that, even though we reduced our model to the DFN model, in fact, if we take the microstructure from the latter (a single particle surrounded by electrolyte) we will find that the conductivity is zero as the particles are not in contact with each other. In the DFN model, this is circumvented by taking an effective value for the conductivity instead of calculating it from the microscale problem.

5.2.2 Conservation of mass and charge in the electrolyte

Finally, we consider the equations in the electrolyte. Here is where the main difference arises: the DFN model uses concentrated electrolyte theory while we use dilute electrolyte theory (see [29] for details on both). Therefore, we cannot directly reduce the model presented here to the DFN model, but we can show that they have similar forms. This is useful as it points out that the concentrated electrolyte equations could be homogenised following a very similar method to the one presented here.

We start assuming isotropy, so we can take the coefficients to be scalars. We define \({\mathcal {D}}_{L}=D_{L} {\mathcal {B}}\), \({\mathcal {D}}_{A}=D_{A} {\mathcal {B}}\), \({\mathcal {M}}_{L}=\mu _{L} {\mathcal {B}}\), and \({\mathcal {M}}_{A}=\mu _{A} {\mathcal {B}}\), and due to isotropy we have \({\mathcal {B}} = B{\mathcal {I}} \) where B is a scalar that captures the microstructure effect on the transport properties. For example, if we used the Bruggeman correlation, then \(B = \varphi _\mathrm {e}^{1.5}\) (see [4] for details). The transference number of the lithium ions is defined as

$$\begin{aligned} t^+ = \frac{\mu _\mathrm {L}}{\mu _\mathrm {L} + \mu _\mathrm {A}}, \end{aligned}$$
(65)

and given that for dilute electrolyte theory the mobilities and diffusivities are proportional to each other, we can replace the mobilities in the definition above for diffusivities. Multiplying (3b) by \(t^+\) and substracting it from (3a) we find

$$\begin{aligned} \varphi _{e}\partial _{t}c_{e}-\nabla \cdot \left( 2 D_\mathrm {L} B (1-t^+) \nabla c_{e}\right) =(1-t^+) a G g. \end{aligned}$$
(66)

This is the same form as the DFN model by setting the effective ion diffusion coefficient to

$$\begin{aligned} D_\mathrm {e,eff} = 2 D_\mathrm {L} B (1-t^+). \end{aligned}$$
(67)

As explained earlier, B captures the microstructure effects. The \(2 (1-t^+)\) factor captures the migration effects in the effective diffusion, but for the concentrated electrolyte this term should be \((1 - 2t^+)\). This discrepancy has been deeply discussed and resolved in [31].

For the conservation of charge equation we can define

$$\begin{aligned}&\kappa _\mathrm {D,eff} = (D_\mathrm {L} - D_\mathrm {A}) c_\mathrm {e}\,\kappa _\mathrm {eff} = (\mu _\mathrm {L} + \mu _\mathrm {A}) c_\mathrm {e}, \end{aligned}$$
(68)

so (3b), together with (3e), can be rearranged into

$$\begin{aligned} \nabla \cdot \left( \kappa _\mathrm {D,eff} \nabla \log c_\mathrm {e}+ \lambda \kappa _\mathrm {eff} \nabla \varPhi _\mathrm {e}\right) = - a G g, \end{aligned}$$
(69)

which is the same form presented in [3]. However, this presents the same inconsistencies between dilute and concentrated electrolyte theories discussed in [31]. Therefore, in order to obtain the effective equations for a concentrated electrolyte one can use the same homogenisation method used in this article with the concentrated electrolyte equations.

6 Discussion

We derived an effective model for the thermal-electrochemical behaviour of porous electrode batteries (with a particular interest in lithium ion batteries) using the method of asymptotic homogenisation. We started from the governing equations at the microscale, which impose mass and charge conservation both in the electrode and the electrolyte, with an ion intercalation reaction at the electrode–electrolyte interface modelled by the Butler-Volmer equation. We assumed that transport of lithium in the electrode is governed by diffusion only and that Ohm’s law holds for the charge. In the electrolyte, we used Nernst-Planck equations to describe the transport phenomena, assuming thus a dilute electrolyte. For the thermal model we assumed heat is only driven by diffusion, and it is generated in the bulk of the material due to Joule heating and at the electrode–electrolyte interface due to the chemical reaction (we included both reversible and irreversible effects).

We exploited the disparity of length scales of the problem and took the limit of infinitely small electrode particles compared to electrode thickness. In this limit, we can derive the homogenised problem at the cell level (mesoscale). A particularity of the analysis is that we took the diffusion coefficient of lithium in the electrode to be small, in order to retrieve the diffusion of lithium only at the microscale, as observed in lithium ion batteries [14, 15]. After deriving the mesoscale model for an electrode, we assembled the cell model (the model at the separator is a particular case of electrode model) and homogenised it to obtain the battery model, exploiting the limit of an infinitely thin cell compared to the battery.

The homogenised model presented here is a generalised version of the DFN model [1], so widely used in the modelling literature. The main differences between our model and the DFN model are that our model includes thermal and capacitance effects, and can account for an arbitrary microstructure. Even though we have not explicitly introduced the temperature dependence of the parameters, as it would clutter the notation, the analysis presented here easily extends for that case and the same homogenised equations are obtained but with temperature-dependent tensors. With respect to the microstructure, in our model we can determine the effective mesoscale parameters directly from the microscopic properties, instead of using theoretical or empirical correlations, such as the Bruggeman correlation [4], which are a subject of debate within the battery modelling field. The model presented here also extends the literature on derivation of multiscale models for batteries using asymptotic homogenisation [21,22,23,24] because it includes thermal effects. In addition, it also captures the diffusion of lithium at the microscale.

The model presented here can be used to simulate electrodes with more realistic geometries, given that it can deal with microstructures more complex than the spherical particles used in the DFN model. For example, one could use SEM images of real electrodes to define the geometry of the cell problems and then solve them numerically. The method used here can also be used in a very similar way to homogenise a concentrated electrolyte model. Another area of future work is to include the current collectors into the model. Also, our model can be used as a stepping stone towards including more complex physical effects that occur in batteries which depend both on electrochemistry and temperature, such as degradation. Several degradation mechanisms that occur in batteries have been considered in the literature [32], but in many cases the degradation models have been coupled to the DFN model in a rather ad hoc way. The framework presented here could be extended to derive degradation models from the microscale models in a systematic and physically consistent way.