A modelling framework for efficient reduced order simulations of parametrised lithium-ion battery cells

Abstract

In this contribution, we present a modelling and simulation framework for parametrised lithium-ion battery cells. We first derive a continuum model for a rather general intercalation battery cell on the basis of non-equilibrium thermodynamics. In order to efficiently evaluate the resulting parameterised non-linear system of partial differential equations, the reduced basis method is employed. The reduced basis method is a model order reduction technique on the basis of an incremental hierarchical approximate proper orthogonal decomposition approach and empirical operator interpolation. The modelling framework is particularly well suited to investigate and quantify degradation effects of battery cells. Several numerical experiments are given to demonstrate the scope and efficiency of the modelling framework.

1. Introduction

Lithium-ion batteries (LIBs) are a key component of our modern society, with applications ranging from medical devices via consumer electronics to electric vehicles and aerospace industry. The further development of LIBs is based on various aspects, namely safety, cost, storage capacity and degradation stability. This research and development are assisted by mathematical models, which are capable of simulating the complex behavior of LIBs on various degrees of spatial and temporal resolution [9, 19, 40, 43, 50, 56, 59, 60, 61, 74, 75]. Mathematical models based on continuum thermodynamics allow, for example, the simulation of charging and discharging processes (cycling), yielding the cell voltage E as function of the capacity Q (or status of charge y) and the cycling rate $C_h$ .¹ These quantities are typically determined in galvanostatic electrochemical measurements, enabling a comparison between theoretical predictions and experimental data [50]. Such models can then be used to investigate and quantify degradation effects of LIBs [30], which are experimentally studied with periodic charging and discharging over N cycles. Experimental measurements seek to determine the dependency of the cell voltage E as function of cycling rate $C_h$ , cycle number n and status of charge, i.e. $E = E(y;\,C_h, n)$ , to systematically quantify different ageing aspects. This requires in general a huge batch of cells and massive amounts of measuring times, e.g. cycling a cell at $C_h=1$ for 1000 cycles requires at least 80 days of continuous electrochemical measurements. Variations of materials or material combinations subsequently increase measuring times exponentially.

Model-based simulations can help to reduce this lab time by reducing the number of batch cells, cycle numbers and material combinations. Further, numerical simulations provide a methodology to systematically deduce informations on a specific ageing effect from electrochemical data. While experimentally the superposition of various degradation effects is always observed, mathematical models allow to predict electrochemical data for single and multiple degradation phenomena. Such ageing effects can be represented in continuum models with different approaches [3, 71, 39], either by full spatial and temporal resolution of a specific effect, e.g. cracking due to intercalation stress, or by cycle number n dependent parameters, which is the focus of this work. This approach requires an evolution equation for some parameters with respect to the cycle number n, which can itself either be upscaled from some microscopic model or determined from experimental snapshots. However, for some prescribed cycle number dependent parameters, e.g. for the solid-state diffusion coefficient of lithium in the intercalation material (modelling the degradation of the intercalation particles) or the intercalation reaction rate (quantifying an increasing surface resistivity), the electrochemical data $E = E(y;\,C_h, n)$ of a battery cell can be numerically computed and provides a fingerprint of the underlying proposed degradation effect (forward problem).² Overall this requires repeated numerical computations of the underlying thermodynamic model and thus efficient numerical strategies. Modern mathematical tools allow to reduce the computational time of coupled non-linear partial differential equations significantly by setting up a reduced basis (RB) method. The combination of (continuum) mathematical modelling, parametrised degradation functions, numerical simulations and reduced basis methods yield a versatile toolbox for the investigation of battery ageing on the basis of electrochemical data.

1.1. Outline

The aim of the paper is to (i) derive a thermodynamically consistent homogenised battery model framework (Section 2), (ii) reduce the computation time using model order reduction techniques (Section 3) and (iii) provide some numerical examples of phenomenological degradation where reduced basis methods are useful (Section 4).

The mathematical model framework for a rather general intercalation battery is derived on the basis of non-equilibrium thermodynamics [16] in Section 2. Our scope is (a) to be fully thermodynamically consistent, (b) provide consistent boundary conditions and (c) to be material independent. The latter aspect allows the applicability of our model framework to various intercalation materials types and cell chemistries. The model considers three porous phases, namely a porous intercalation anode, a porous separator phase and a porous intercalation cathode. Each porous phase consists of an electrolyte phase, an active phase and a conductive additive phase. The electrolyte phase is based on a rigorously validated material model [52], accounting for solvation effects [23], incompressibility constraints, diffusion and conductivity. The active intercalation phase of the electrodes accounts for intercalation of (exemplarily) lithium in terms of a specific chemical potential function and solid-state diffusion with a concentration dependent diffusion coefficient [50]. The conductive additive phase accounts for the electron transport. All three phases are coupled through a reaction boundary condition, where a special emphasis is put on thermodynamic consistency [25, 48]. Subsequent spatial homogenisation techniques [1, 15, 43] for the porous structure yield an effective, non-linear coupled partial differential equation (PDE) system for the lithium concentration in each phase, the lithium-ion concentration in the electrolyte and the electrostatic potential in each phase.

Building on this, we present a reduced basis method for recurring numerical simulations of the parameterised PDE model in Section 3, which occur in the simulation repeated battery cycling under laboratory conditions. As a first step, we discretise the resulting fully nonlinear system of PDEs by the finite element method in space and the backward Euler method in time. The computational studies to identify critical parameters for estimating the ageing process of batteries require many evaluations of the finite element system with different parameter settings and thus involve a large amount of time and experimental effort. Therefore, we employ the RB method in order to further reduce the parameterised discretised battery model to obtain a reduced order model (ROM) that is cheap to evaluate. For an introduction and overview on recent developments in model order reduction, we refer to the monographs and collections [6, 7, 41, 73]. For time-dependent problems, the proper orthogonal decomposition (POD)-Greedy method [37] defines the Gold-Standard for systems, where rigorous and cheap to evaluate a posteriori error estimates are available. As this is not the case for the non-linear battery model at hand, we employ the hierarchical approximate proper orthogonal decomposition (HAPOD) [42] in this contribution. As RB methods rely on so called efficient offline/online splitting, they need to be combined with supplementary interpolation methods in case of non-affine parameter dependence or non-linear differential equations. The empirical interpolation method (EIM) [2] and its various generalisations are key technologies with this respect. In this contribution, we apply the empirical operator interpolation as introduced in [28, 58] to the full battery finite element operator. Several related model order reduction approaches have already been applied for battery simulation, as, e.g. in [12, 44, 55], where the authors apply model reduction techniques for Newman-type LIB models [19]. Further results on model order reduction in the context of battery models, including resolved electrode geometry and Butler–Volmer kinetics can be found in [31, 68, 69, 70, 82] and [84]. In addition, there are several solvers for Newman-type battery models, as e.g. [8, 46, 78].

We provide finally in section Section 4 some numerical examples of phenomenological degradation models, where the RB method yields the desired reduction of computation time. These phenomenological degradation models are expressed in terms of cycle number n dependent parameters, e.g. diffusion coefficients or exchange current densities, which can be linked to specific underlying degradation mechanisms. However, this is not the scope of this work, but rather showing the qualitative impact of the parameter variation on the electrochemical data, i.e. $E = E(y;\,C_h, n)$ , and its efficient and fast computation with RB methods. Our work is summarised in Section 5, and an overview of all variables and parameters is given Section 6 and Appendix B, respectively.

2. Derivation of the porous electrode model

In this section, we derive a continuum mathematical model for a porous battery cell. After setting up domains and proper definitions in Section 2.1, we state the chemical potential functions for all phases in Section 2.1.3 and briefly discuss their derivation as well as some consequence of the electro-neutrality condition. In Section 2.1.2, we state then the corresponding transport equations for each phase in the porous electrode, where Section 2.1.3 puts a special emphasis on the intercalation reaction boundary condition and its formulation on basis of non-equilibrium surface thermodynamics. Section 2.1.4 covers then the full PDE system of an unhomogenised porous electrode, setting the basis for spatial homogenisation in Section 2.2. Introducing proper scalings in Section 2.2.2 yields a general homogenised equation framework in Section 2.2.3 for a porous multi-phase electrode. Section 2.2.4 then covers the full homogenised battery model, where proper scalings are introduced on the basis of the $1-$ C current density. An overview of all parameters is given in Section 2.3.1.

2.1. Domains, definitions and initial scaling

We seek to model a porous electrochemical unit cell, consisting of a porous anode $\Omega^{\text{An}} \subset \mathbb{R}^3$ , a porous separator $\Omega^{\text{Sep}} \subset {\mathbb{R}}^3$ and a porous cathode $\Omega^{\text{Cat}} \subset {\mathbb{R}}^3$ , which are packed between two metallic deflectors $\Omega^{\text{Def}_{\text{An}}}$ and $\Omega^{\text{Def}_{\text{Cat}}}$ that yield the overall electrical connections (see Figure 1). The base area of the deflectors is denoted by $\Sigma$ . The intercalation electrodes $\Omega^j, j = {\text{An}},{\text{Cat}}$ consist themself of three phases, namely the active particle phase $\Omega_{\texttt{A}}^j$ , the conductive additive phase $\Omega_{\texttt{C}}^j$ and the electrolyte phase $\Omega_{\texttt{E}}^j$ , with $\Omega^j = \bigcup_{i \in \{{\texttt{A}},{\texttt{C}},{\texttt{E}} \}} \Omega_i^j$ . The union of the active phase and the conductive additive is frequently called solid phase $\Omega^j_{\texttt{S}} = \bigcup \Omega_{\texttt{A}}^j \cup \Omega_{\texttt{C}}^j$ . The porous separator consists of an electrolyte phase $\Omega_{\texttt{E}}^{\text{Sep}}$ and an polymeric additive $\Omega_{\texttt{P}}^{\text{Sep}}$ , with $\Omega^{\text{Sep}} = \Omega_{\texttt{E}}^{\text{Sep}} \cup \Omega_{\texttt{P}}^{\text{Sep}}$ . The whole electrolyte phase is further denoted by $\Omega_{\texttt{E}} = \bigcup_{j \in \{{\text{An}},{\text{Sep}},{\text{Cat}}\}} \Omega_{\texttt{E}}^j$ , the active phase as $\Omega_{\texttt{A}} = \bigcup_{j \in \{{\text{An}},{\text{Cat}}\}} \Omega_{\texttt{E}}^j$ and the solid phase $\Omega_{\texttt{S}} = \bigcup_{j \in \{{\text{An}},{\text{Cat}}\}} \Omega_{\texttt{S}}^j$ . The interface $\Sigma_{{\texttt{A}},{\texttt{E}}} = \Omega_{\texttt{A}} \cap \Omega_{\texttt{E}}$ captures the actual surface $\Sigma_{\texttt{A}}$ of the active particle as well as the electrochemical double layer ( $\Omega_{\texttt{A}}^{\text{SCL}}, \Omega_{\texttt{E}}^{\text{SCL}}$ ) forming at the interface, i.e. $\Sigma_{{\texttt{A}},{\texttt{E}}} = \Omega_{\texttt{A}}^{\text{SCL}} \cup \Sigma_{\texttt{A}} \cup \Omega_{\texttt{E}}^{\text{SCL}}$ . The domains $\Omega_{\texttt{E}}$ and $\Omega_{\texttt{A}}$ are thus electro-neutral, and we refer to [25, 48, 63] for details on the derivation. A similar argument holds for the interface $\Sigma_{{{\texttt{C}}{\texttt{E}}}} = \Omega_{\texttt{C}} \cap \Omega_{\texttt{E}}$ . The interface between the deflectors and a single phase of the porous electrode is denoted by $\Sigma^i_{j{\texttt{D}}} = \Omega^i_j \cup \Omega^{\text{Def}_i}, i = {\text{An}},{\text{Cat}}, j = {\texttt{A}},{\texttt{C}},{\texttt{E}}$ , for example, m denotes $\Sigma^{\text{Cat}}_{{\texttt{C}}{\texttt{D}}}$ the surface through which electrical current from the conductive phase of the electrode flows to the respective deflector.

Figure 1. Sketch of the porous electrochemical unit cell. During discharge, lithium ions flow from the anode to the cathode, while electrons drive an external electrical consumer.

In each phase of each domain, we have balance equations, which are coupled through respective boundary conditions modelling the intercalation reaction. We seek to employ periodic homogenisation to derive a homogenised PDE model for the electrochemical unit cell $\Omega = \Omega^{\text{Cat}} \cup \Omega^{\text{Sep}} \cup \Omega^{\text{An}}$ . In order to do so, we require a specific scaling with respect to the two essential length scales arising in the problem, (i) the length scale $\mathscr{L}$ of the macroscopic width W of the unit cell, i.e. $\Omega = \Sigma \times [0,W]$ with $\Sigma \subset {\mathbb{R}}^2$ being the area of the deflectors and (ii) the length scale $\ell$ of the intercalation particle radii $r_{\texttt{A}}$ . This yields a small parameter $\varepsilon = \frac{\mathscr{L}}{\ell}$ with respect to which we can perform a multi-scale asymptotic expansion and thus a periodic homogenisation.

2.1.1. Variables, chemical potentials and parameter (equilibrium)

The active particle $\Omega_{\texttt{A}}$ is a mixture of electrons $e^{-}$ , intercalated cations C and lattice ions ${M^+}$ and the electrolyte a mixture of solvated cations ${C^+}$ , solvated anions ${A^-}$ and solvent molecules S. The respective species densities are denoted with $n_\alpha({\textbf{x}},t), {\textbf{x}} \in \Omega_i$ . Further, we denote with

(2.1) \begin{align} \mu_\alpha = \frac{\partial \psi}{\partial n_\alpha} , \qquad \qquad i={\texttt{A}},{\texttt{E}}, \alpha = {{\texttt{E}}_A},{{\texttt{E}}_C},{{\texttt{E}}_S},\texttt{A}_{C},{{{\texttt{A}}_e}},{{{\texttt{A}}_M}},\end{align}

the chemical potential of the constituents, which are derived from a free energy density [16, 62] $\psi = \psi_{\texttt{A}} + \psi_{\texttt{E}}$ with $\psi_{\texttt{A}}=\hat \psi_{\texttt{A}}\!\left(n_{{{\texttt{A}}_e}},n_{{{\texttt{A}}_C}},n_{{{\texttt{A}}_M}}\right)$ of the active particle and $\psi_{\texttt{E}}=\hat \psi\!\left(n_{{\texttt{E}}_S},n_{{\texttt{E}}_A},n_{{\texttt{E}}_C}\right)$ of the electrolyte phase [27, 48, 52, 62].

Electrolyte

For the electrolyte, we consider exclusively the material model [23, 24, 52] of an incompressible liquid electrolyte accounting for solvation effects, i.e.

(2.2) \begin{align} \mu_\alpha = \psi_{\alpha,\tiny\text{ref}} + k_{\tiny\text{B}}T\, \text{ln}\!\left( {y_\alpha} \right) + v_\alpha \cdot p \quad \alpha = {{\texttt{E}}_S},{{\texttt{E}}_A},{{\texttt{E}}_C},\end{align}

with mole fraction

(2.3) \begin{align} y_\alpha = \frac{n_\alpha}{n_{{\texttt{E}},\tiny\text{tot}}}, \sum_\alpha y_\alpha = 1,\end{align}

molar concentration $n_\alpha$ , total molar concentration of the mixture (with respect to the number of mixing particles [52])

(2.4) \begin{align} n_{{\texttt{E}},\tiny\text{tot}} = n_{{\texttt{E}}_S} + n_{{\texttt{E}}_A} + n_{{\texttt{E}}_C},\end{align}

partial molar volume $v_\alpha$ of the specific constituent within the mixture, pressure p and reference molar free energy $\psi_{\alpha,\tiny\text{ref}}$ of the pure substance in the mixture, temperature T and Boltzmann constant $k_{\tiny\text{B}}\,$ .

A major aspect of the material model is the solvation effect, where each ion binds $\kappa_\alpha$ solvent molecules due to microscopic electrostatic interactions. These aspects are crucial for various aspects of the thermodynamic model and we refer to [21, 23, 49, 52] for its details as well as experimental validation. The bound solvent molecules do not contribute to the entropy of mixing [23], whereby $n_{{\texttt{E}}_S}$ in equation (2.4) denotes the number density of free ³ solvent molecules in the mixture. In turn, the partial molar volume $v_\alpha$ of the ionic constituents are increased roughly by $\kappa_\alpha \cdot v_{{\texttt{E}}_S}$ compared to the the volume of the bare central ion. Consider a single ion with a bare molar mass and molar volume of $m_{\alpha,\text{bare}}$ and $v_{\alpha,\text{bare}}$ , respectively. The corresponding solvated ion binds $\kappa_\alpha$ solvent molecules with mass and volume $m_{{{\texttt{E}}_S}}$ and $v_{{{\texttt{E}}_S}}$ , whereby the mass of the solvated ion is $m_\alpha = m_{\alpha,\text{bare}} + \kappa_\alpha m_{{\texttt{E}}_S}$ . A quite similar relation also holds for the volume of the ion, however, since the volume is not necessarily conserved upon solvation, for instance due to packing density aspects, such a relation is approximative, i.e. $v_\alpha \approx v_{\alpha,\text{bare}} + \kappa_\alpha v_{{\texttt{E}}_S}$ . From a thermodynamic point of view, it is convenient and meaningful to assume

(2.5) \begin{align} \frac{v_{{\texttt{E}}_C}}{v_{{\texttt{E}}_S}} = \frac{m_{{\texttt{E}}_C}}{m_{{\texttt{E}}_S}} \quad \text{and} \quad \frac{v_{{\texttt{E}}_A}}{v_{{\texttt{E}}_S}} = \frac{m_{{\texttt{E}}_A}}{m_{{\texttt{E}}_S}},\end{align}

whereby the incompressibility constraint [23, 24, 52] implies also a conservation of mass, i.e.

(2.6) \begin{align} \sum_{\alpha} v_\alpha n_\alpha = 1 \quad \Leftrightarrow \quad \sum_{\alpha} m_\alpha n_\alpha = \rho = \frac{m_{{\texttt{E}}_S}}{v_{{{\texttt{E}}_S},\tiny\text{ref}}} = \text{const.}\end{align}

If the molar volume (and mass) of the solvent molecules are far larger than the bare ion,⁴ i.e. $v_{{\texttt{E}}_S} \gg v_{\alpha,\text{bare}}$ , we may approximate further

(2.7) \begin{align} v_{{\texttt{E}}_C} = \kappa_{\texttt{E}} \cdot v_{{\texttt{E}}_S} \quad \text{and} \quad v_{{\texttt{E}}_A} = \kappa_{\texttt{E}} \cdot v_{{\texttt{E}}_S},\end{align}

whereby the volume of the solvated ion is mainly determined by the solvation shell and the solvation number is similar [52] for anion and cation ( $\kappa_{{\texttt{E}}_A} =\kappa_{{\texttt{E}}_C}=\!:\, \kappa_{\texttt{E}}$ ). The molar volume of the solvent is related to the mole density $n_{{{\texttt{E}}_S},\tiny\text{ref}}$ of the pure solvent as

(2.8) \begin{align} v_{{\texttt{E}}_S} = (n_{{{\texttt{E}}_S},\tiny\text{ref}})^{-1}\end{align}

and the solvation numbers $\kappa_\alpha$ can be determined from differential capacity measurements [52]. Note that the incompressibility constraint equation (2.6) is also used to determine the number densities $n_\alpha$ in terms of the mole fractions $y_\alpha$ , i.e.

(2.9) \begin{align}n_\alpha = \frac{y_\alpha}{n} \overset{\text{equation (2.6) }} = \frac{y_\alpha}{\sum_{\beta} v_\beta y_\beta}.\end{align}

Finally, the electrostatic potential is an additional variable in the electrolyte phase and denoted as $\varphi_{\texttt{E}}({\textbf{x}},t)$ .

Active Phase

As an example, we discuss an electrode active phase, which is applied to the anode and cathode phase in Section 2.1.4. For the active particle phase $\Omega_{\texttt{A}}$ , we consider exemplarily a classical lattice mixture model [4, 10, 11, 20, 22, 34, 51, 53], which we term in the following ideal electrode material. Note that the whole modelling as well as the numerical procedure can directly be adapted to other chemical potential functions $\mu_{{{\texttt{A}}_C}} = \mu_{{{\texttt{A}}_C}}(y_{\texttt{A}})$ . The chemical potential of intercalated lithium is stated as

(2.10) \begin{align} \mu_{{{\texttt{A}}_C}} = g_{{{\texttt{A}}_C}} + k_{\tiny\text{B}}T\, f_{\texttt{A}}(y_{\texttt{A}}) , \quad \text{with} \quad f_{\texttt{A}}(y_{\texttt{A}}) &\,:\!=\, \text{ln}\!\left( { \frac{y_{\texttt{A}}}{1\! +\! y_{\texttt{A}}} } \right) + \gamma_{\texttt{A}} \cdot (2 y_{\texttt{A}} \!-\! 1),\end{align}

with mole fraction

(2.11) \begin{align} y_{\texttt{A}} = \frac{n_{{{\texttt{A}}_C}}}{n_{{{\texttt{A}},\tiny\text{lat}}}} \end{align}

of intercalated cations in the active phase and partial molar Gibbs energy $g_{{{\texttt{A}}_C}} = \text{const.}$ The number density $n_{{{\texttt{A}},\tiny\text{lat}}}$ of lattice sites is constant, which corresponds to an incompressible lattice, and the enthalpy parameter $\gamma_{\texttt{A}} > -2.5$ . Note that $\gamma_{\texttt{A}} < - 2.5$ entails a phase separation [53]. The electrostatic potential in the solid phase is denoted as $\varphi_{\texttt{S}}({\textbf{x}},t)$ .

Electro-neutrality

The electro-neutrality condition of $\Omega_{\texttt{A}}$ , $\Omega_{\texttt{E}}$ and $\Omega_{\texttt{C}}$ can be obtained by an asymptotic expansion of the balance equations in the electrochemical double layer at the respective surface $\Sigma$ . We only briefly recapture the central conclusions and refer to [25, 26, 48, 50, 52, 63] for details on the modelling, validation and the asymptotic derivation. Most importantly, electro-neutrality implies (i) that the double layer is in thermodynamic equilibrium, (ii) that there exists a potential jump through the interface $\Sigma_{{{\texttt{A}}{\texttt{E}}}}$ and (iii) that for monovalent electrolytes the cation mole fraction (or number density) is equal to the anion mole fraction, i.e.

(2.12) \begin{align}y_{{\texttt{E}}_C} = y_{{\texttt{E}}_A} =\!:\, y_{\texttt{E}} .\end{align}

The total number density $n_{{\texttt{E}},\tiny\text{tot}}= n_{{\texttt{E}}_S} + n_{{\texttt{E}}_C} + n_{{\texttt{E}}_A}$ and the electrolyte concentration $n_{{\texttt{E}}_C}$ can be written as a function of $y_{\texttt{E}}$ using equation (2.6) as follows:

(2.13) \begin{align} n_{{\texttt{E}},\tiny\text{tot}} & = n_{{\texttt{E}}_S} \cdot \frac{1}{1+2(\kappa_{\texttt{E}}-1) y_{\texttt{E}}} = n_{{\texttt{E}},\tiny\text{tot}}(y_{\texttt{E}}) \end{align}

(2.14) \begin{align} n_{{\texttt{E}}_C} &= n_{{\texttt{E}}_S} \frac{y_{\texttt{E}}}{1+2(\kappa_{\texttt{E}}-1) y_{\texttt{E}}} = n_{{\texttt{E}}_C}(y_{\texttt{E}}).\end{align}

2.1.2. Transport equation relations

In the electrolyte $\Omega_{\texttt{E}}$ , we have two balance equations determining the concentration $n_{{\texttt{E}}_C}({\textbf{x}},t)$ (or mole fraction $y_{\texttt{E}}({\textbf{x}},t)$ ) and the electrostatic potential $\varphi_{\texttt{E}}({\textbf{x}},t)$ in the electrolyte [32, 33, 57, 64, 65, 66,], i.e.

(2.15) \begin{align} \frac{\partial n_{{\texttt{E}}_C}}{\partial t} &= -{\text{div}} {\textbf{J}}_{{\texttt{E}}_C} \quad \text{with} \quad & {\textbf{J}}_{{\texttt{E}}_C} &= - D_{\texttt{E}} \cdot n_{{\texttt{E}},\tiny\text{tot}}\, \Gamma_{\texttt{E}}^{\text{tf}} \cdot \nabla y_{\texttt{E}} + \frac{t_{{\texttt{E}}_C}}{e_0} {\textbf{J}}_{\texttt{E}_{q}} \end{align}

(2.16) \begin{align} \qquad\ 0 &= -{\text{div}} {\textbf{J}}_{{\texttt{E}}_q} \quad \text{with} \quad\ \ \ & {\textbf{J}}_{{\texttt{E}}_q} &= - S_{\texttt{E}} \cdot n_{{\texttt{E}},\tiny\text{tot}}\, \Gamma_{\texttt{E}}^{\text{tf}} \nabla y_{\texttt{E}} - \Lambda_{\texttt{E}} n_{{\texttt{E}}_C} \nabla \varphi_{\texttt{E}} \end{align}

effective electrolyte thermodynamic potential⁵

(2.17) \begin{align} f_{\texttt{E}}(y_{\texttt{E}}) \,:\!=\, \text{ln}\!\left( {y_{\texttt{E}}} \right) - \kappa_{\texttt{E}} \text{ln}\!\left( {(1-2\cdot y_{\texttt{E}})} \right) \end{align}

and (dimensionless) thermodynamic factor [24]

(2.18) \begin{align} \Gamma_{\texttt{E}}^{\text{tf}}(y_{\texttt{E}}) &= y_{\texttt{E}} \frac{\partial f_{\texttt{E}}}{\partial y_{\texttt{E}}} = 1 + 2 \kappa_{\texttt{E}} \frac{y_{\texttt{E}}}{1 - 2 y_{\texttt{E}}}, \end{align}

and electrolyte diffusion coefficient $D_{\texttt{E}}$ , transference number $t_{{\texttt{E}}_C}$ , molar electrolyte conductivity $\Lambda_{\texttt{E}}$ , electrolytic chemical conductivity $S_{\texttt{E}}$ , charge number $z_\alpha$ and elementary charge $e_0$ . In general, three of the transport parameters are independent and $S_{\texttt{E}}$ , $t_{{\texttt{E}}_C}$ and $\Lambda_{\texttt{E}}$ are connected via

(2.19) \begin{align} \frac{k_{\tiny\text{B}}T\,}{e_0} (2 t_C - 1) = \frac{S_{\texttt{E}}}{\Lambda_{\texttt{E}}} . \end{align}

The flux representations of (2.15) and (2.16) can be deduced on the basis of general Maxwell–Stefan type diffusional fluxes [32, 33, 45, 57, 64] or Nernst–Planck-type flux relations with cross-diffusional terms [16, 56]. For a three-component mixture $({{\texttt{E}}_S},{{\texttt{E}}_A},{{\texttt{E}}_C})$ , the Onsager reciprocal relations [16] state that two independent fluxes are present, i.e.

(2.20) \begin{align} {\textbf{J}}_{{\texttt{E}}_A} &= M_{AA} \nabla \tilde \mu_{{\texttt{E}}_A} + M_{AC} \nabla \tilde \mu_{{\texttt{E}}_C}\end{align}

(2.21) \begin{align} {\textbf{J}}_{{\texttt{E}}_C} &= M_{CA} \nabla \tilde \mu_{{\texttt{E}}_A} + M_{CC} \nabla \tilde \mu_{{\texttt{E}}_C}\end{align}

(2.22) \begin{align} & \text{ with } \tilde \mu_\alpha = \mu_\alpha - \frac{m_\alpha}{m_{{\texttt{E}}_S}} \mu_{{\texttt{E}}_S} + e_0 z_\alpha \varphi, \alpha = {{\texttt{E}}_A},{{\texttt{E}}_C},\end{align}

with mobilities $M_{\alpha\beta}$ , satisfying $M_{AC}= M_{CA}$ . Hence, three transport parameters are independent, i.e. $(M_{AA}, M_{AC}, M_{CC})$ . Rewriting the fluxes ${\textbf{J}}_{{\texttt{E}}_A}$ and ${\textbf{J}}_{{\texttt{E}}_C}$ , together with the charge flux

(2.23) \begin{align} {\textbf{J}}_{{\texttt{E}}_q} =e_0 z_{{\texttt{E}}_C} {\textbf{J}}_{{\texttt{E}}_C} + e_0 z_{{\texttt{E}}_A} {\textbf{J}}_{{\texttt{E}}_A} ,\end{align}

in the representations of (2.15) and (2.16), yields the definitions of $(D_E, t_E,S_E, \Lambda_E)$ in terms of the mobilities $(M_{AA}, M_{AC}, M_{CC})$ , which are given in Appendix A, as well as and the condition

(2.24) \begin{align} \frac{k_B T}{e_0} (2 t_{{\texttt{E}}_C} - 1) = \frac{S_{\texttt{E}}}{\Lambda_{\texttt{E}}} .\end{align}

We emphasise, however, that for simple Nernst–Planck-fluxes [24, 76], i.e.

(2.25) \begin{align} {\textbf{J}}_\alpha = D_\alpha^{\text{NP}} \frac{n_\alpha}{k_{\tiny\text{B}}T\,}\! \left( \nabla \mu_\alpha - \frac{m_\alpha}{m_{{\texttt{E}}_S}} \nabla \mu_{{\texttt{E}}_S} + e_0 z_\alpha n_\alpha \nabla \varphi_{\texttt{E}}\right) \qquad \alpha = {{\texttt{E}}_A},{{\texttt{E}}_C}, \end{align}

with constant diffusion coefficients $D_\alpha^{\text{NP}}$ , only two transport parameters are independent and the transport parameters of of (2.15) and (2.16) compute as

(2.26) \begin{align} D_{\texttt{E}} &= \frac{2 D_{{\texttt{E}}_C}^{\text{NP}} \cdot D_{{\texttt{E}}_A}^{\text{NP}}}{D_{{\texttt{E}}_A}^{\text{NP}} + D_{{\texttt{E}}_C}^{\text{NP}}} = \text{const.} & \quad t_{{\texttt{E}}_C} &= \frac{D_{{\texttt{E}}_C}^{\text{NP}}}{D_{{\texttt{E}}_A}^{\text{NP}} + D_{{\texttt{E}}_C}^{\text{NP}}} = \text{const.} \end{align}

(2.27) \begin{align} \Lambda_{\texttt{E}} &= \frac{e_0^2}{k_{\tiny\text{B}}T\,}\left(D_{{\texttt{E}}_A}^{\text{NP}} + D_{{\texttt{E}}_C}^{\text{NP}}\right) = \text{const.} & S_{\texttt{E}} &= e_0\!\left(D_{{\texttt{E}}_C}^{\text{NP}} - D_{{\texttt{E}}_A}^{\text{NP}}\right)= \text{const.} \end{align}

For the numerical simulations of Section 3, we assume constant values for the transport parameters $(t_{{\texttt{E}}_C}, S_{\texttt{E}}, D_{\texttt{E}},\Lambda_{\texttt{E}})$ , which is sufficient to show the general methodology of our framework approach. Note, however, that $(t_{{\texttt{E}}_C}, S_{\texttt{E}}, D_{\texttt{E}},\Lambda_{\texttt{E}})$ can in general be depend on the electrolyte concentration $n_{{\texttt{E}}_C}$ , which can be easily incorporated our framework.

In the active particle $\Omega_{\texttt{A}}$ , we have two balance equations determining the concentration $n_{{{\texttt{A}}_C}}({\textbf{x}},t)$ (or mole fraction $y_{\texttt{A}}$ ) of intercalated lithium and the electrostatic potential $\varphi_{\texttt{S}}({\textbf{x}},t)$ in the solid phase $\Omega_{\texttt{S}}$ , i.e. the balance equation for intercalated lithium ${\texttt{A}}_C$ and the balance equation for electrons⁶ ${\texttt{A}}_e$ in the entire solid phase,

(2.28) \begin{align} \frac{\partial n_{{{\texttt{A}}_C}}}{\partial t} &= -{\text{div}} {\textbf{J}}_{{{\texttt{A}}_C}} \quad \text{with} \quad & {\textbf{J}}_{{{\texttt{A}}_C}} &= - D_{\texttt{A}} \cdot n_{{{\texttt{A}},\tiny\text{lat}}}\, \Gamma_{\texttt{A}}^{\text{tf}} \cdot \nabla y_{\texttt{A}}, \quad {\textbf{x}} \in \Omega_{\texttt{A}}\end{align}

(2.29) \begin{align} \ \ 0 &= -{\text{div}} {\textbf{J}}_{\texttt{A}_{q}} \quad \text{with} & \ \ {\textbf{J}}_{{{\texttt{A}}_q}} &= - \sigma_{\texttt{S}}({\textbf{x}}) \nabla \varphi_{\texttt{S}}, \quad {\textbf{x}} \in \Omega_{\texttt{S}} \qquad \end{align}

and (dimensionless) thermodynamic factor

(2.30) \begin{align} \Gamma_{\texttt{A}}^{\text{tf}} = y_{\texttt{A}} \frac{\partial f_{\texttt{A}}}{\partial y_{\texttt{A}}} = 1 + \frac{y_{\texttt{A}}}{1- y_{\texttt{A}}} + 2 \gamma_{\texttt{A}} y_{\texttt{A}} = \Gamma_{\texttt{A}}^{\text{tf}}(y_{\texttt{A}}) . \end{align}

The (solid state) diffusion coefficient $D_{\texttt{A}}$ is further assumed to be

(2.31) \begin{align} D_{\texttt{A}} = D_{{\texttt{A}},\tiny\text{ref}} \cdot (1-y_{\texttt{A}}), \qquad D_{{\texttt{A}},\tiny\text{ref}} = \text{const.}, \end{align}

where the term $(1-y_{\texttt{A}})$ accounts for a reduced (solid state) diffusivity due to crowding [50]. Note that the electrical conductivity $\sigma_{\texttt{S}}$ is different in the active phase $\Omega_{\texttt{A}}$ and the conductive phase $\Omega_{\texttt{C}}$ , which form $\Omega_{\texttt{S}} = \Omega_{\texttt{A}} \cup \Omega_{\texttt{C}}$ . We account for this as

(2.32) \begin{align}\sigma_{\texttt{S}}({\textbf{x}}) = \begin{cases} \sigma_{{\texttt{A}}}, &\text{ if }{\textbf{x}} \in \Omega_{\texttt{A}}\\ \sigma_{{\texttt{C}}}, &\text{ if }{\textbf{x}} \in \Omega_{\texttt{C}}\end{cases}. \end{align}

In principle, $\sigma_{\texttt{A}}$ can be dependent on the amount of intercalated ions, i.e. $\sigma_{\texttt{A}} = \sigma_{\texttt{A}}(y_{\texttt{A}})$ , but for the sake of this work we assume $\sigma_{\texttt{A}} = \text{const.}$ and $\sigma_{\texttt{C}} = \text{const.}$

A remark on the diffusion equation for intercalated lithium

We emphasise that the balance equation (2.28) for intercalated lithium is naturally [50] a non-linear diffusion equation, which arises from the non-linear chemical potential function $\mu_{\texttt{A}}(y_{\texttt{A}})$ . For the sake of thermodynamic consistency of the entire model framework, it is necessary that the lithium flux in the solid phase is strictly represented as ${\textbf{J}}_{{{\texttt{A}}_C}} = M_{\texttt{A}}(y_{\texttt{A}}) \nabla \mu_{\texttt{A}}(y_{\texttt{A}})$ , where the mobility is $M_{\texttt{A}}$ is typically considered as $M_{\texttt{A}} = D_{\texttt{A}}(y_{\texttt{A}}) \cdot \frac{n_{\texttt{A}}}{k_{\tiny\text{B}}T\,}$ (Einstein relation [16]).

For the very special case $f_{\texttt{A}} = \text{ln}\!\left( {\frac{y_{\texttt{A}}}{1-y_{\texttt{A}}}} \right)$ , i.e. an ideal lattice gas, one obtains $\Gamma_{\texttt{A}}^{\text{tf}} = \frac{1}{1-y_{\texttt{A}}}$ and thus together with equation (2.31) a simple, linear diffusion equation

(2.33) \begin{align} \frac{\partial y_{\texttt{A}}}{\partial t} &= - D_{{\texttt{A}},\tiny\text{ref}} {\text{div}} \big(\nabla y_{\texttt{A}}\big) \quad \text{for} \ \ {\textbf{x}} \in \Omega_{\texttt{A}},\end{align}

as stated, for example, by Newman [19, 66]. We emphasise, however, that this entails for the OCP (c.f. the next Section 2.1.3 for a more detailed definition) of the active material, which is $E^{OCP} = \frac{1}{e_0}(\mu_{Li} - \mu_{\texttt{A}}(y_{\texttt{A}}))$ [50], also simple relation $E^{OCP} = E^0 - \frac{k_{\tiny\text{B}}T\,}{e_0}\text{ln}\!\left( {\frac{y_{\texttt{A}}}{1-y_{\texttt{A}}}} \right)$ with $E^0=\text{const.}$ This is, however, in contrast to the OCP relation of [66], in our notation $E^{OCP} = E^0 - \frac{k_{\tiny\text{B}}T\,}{e_0}\big(\text{ln}\!\left( {\frac{y_{\texttt{A}}}{1-y_{\texttt{A}}}} \right) + \beta y_{\texttt{A}} + \zeta \big)$ ( $\beta$ and $\zeta$ are considered as parameters in [66]).

2.1.3. Reaction rate and affinity

At the the interface $\Sigma_{{{\texttt{A}}{\texttt{E}}}}$ the intercalation reaction

(2.34) \begin{align} \textrm{Li}^+ \big|_{\texttt{E}} + \textrm{e}^{-}\big|_{\texttt{A}} \rightleftharpoons \textrm{Li}\big|_{\texttt{A}} + \kappa_{\texttt{E}} \cdot \textrm{S}\big|_{\texttt{E}}\end{align}

occurs, which is modelled on the basis of (surface) non-equilibrium thermodynamics [50]. Hence, the (surface) reaction rate ${\underset{s}{R}\,\!}$ can in general be written with $\alpha \in [0,1]$ as [17, 21, 25, 48, 77]

(2.35) \begin{align} {\underset{s}{R}\,\!} = {\underset{s}{L}\,\!} \cdot g\!\left( \frac{1}{k_{\tiny\text{B}}T\,} {\underset{s}{\lambda}\,\!}\right) \quad \text{with} \quad g(z) \,:\!=\, \left(\text{e}^{\alpha \cdot z}\, - \text{e}^{-(1-\alpha) \cdot z}\, \right),\end{align}

where

(2.36) \begin{align}{\underset{s}{\lambda}\,\!} {}_{{{\texttt{A}}{\texttt{E}}}} = e_0 \hat \varphi_{\texttt{S}}|_{{{\texttt{A}}{\texttt{E}}}} - e_0 \hat \varphi_{\texttt{E}}|_{{{\texttt{A}}{\texttt{E}}}} + k_{\tiny\text{B}}T\, \cdot \left( f_{\texttt{A}}^j\!\left(y_{\texttt{A}}|_{{{\texttt{A}}{\texttt{E}}}}\right) - f_{\texttt{E}}\!\left(y_{\texttt{E}}|_{{{\texttt{A}}{\texttt{E}}}}\right) \right)\end{align}

denotes the surface affinity of the reaction (2.34), which is pulled back through the double layer to the respective points (in an asymptotic sense) outside of the double layer. The quantity $\hat \varphi_{\texttt{E}} \,:\!=\, \varphi_{\texttt{E}} - E_{{Li^+},{\texttt{E}}}$ denotes the electrolyte potential vs. metallic Li and $\hat \varphi_{\texttt{S}}^j \,:\!=\, \varphi_{\texttt{S}}^j - E_{{Li^+},{\texttt{A}}}^j$ the electrode $^j$ potential vs. metallic Li [50]. Note again that $y_{\texttt{A}}|_{{{\texttt{A}}{\texttt{E}}}}$ denotes the evaluation of $y_{\texttt{A}}$ at the interface $\Sigma_{{\texttt{A}},{\texttt{E}}}$ and that the surface affinity (2.36) is dependent on the chemical potential (or the mole fraction) evaluated at the interface. The non-negative function ${\underset{s}{L}\,\!}$ in (2.35) ensures a non-negative entropy production ${\underset{s}{r}\,\!} {}_{\sigma,R}$ due to reactions on the surface, i.e. ${\underset{s}{r}\,\!} {}_{\sigma,R} = {\underset{s}{\lambda}\,\!} \cdot {\underset{s}{R}\,\!} > 0$ . For the sake of this work, we assume ${\underset{s}{L}\,\!} = \text{const.}$ and refer to [50] for a detailed discussion on concentration dependency.

A remark on the various formulations of the Butler–Volmer equation.

The reaction rate (2.35) is frequently called Butler–Volmer-type relation; however, various formulations in terms of the precise functional dependency on the thermodynamic state variables are present in the literature [19, 27, 38, 48, 50, 56, 66, 77, 81]. The experimentally derived Butler–Volmer equation [38] states a relation between the global current I and the deviation of the (macroscopic) cell potential E to the (macroscopic) equilibrium voltage $E^{OCP}$ , frequently called open circuit potential (OCP) [66], i.e. $ I = I_0 \cdot \left(\text{e}^{- \alpha \frac{e_0}{k_{\tiny\text{B}}T\,} \left(E - E^{OCP}\right)}\, - \text{e}^{- (1-\alpha) \frac{e_0}{k_{\tiny\text{B}}T\,} \left(E - E^{OCP}\right)}\right)$ . For intercalation materials, this OCP is concentration dependent, i.e. $E^{OCP}=E^{OCP}(y_{\texttt{A}})$ , which arises from the chemical potential function $\mu_{\texttt{A}}(y_{\texttt{A}})$ since in a half cell vs. metallic lithium $E^{OCP} = \frac{1}{e_0}(\mu_{Li} - \mu_{\texttt{A}}(y_{\texttt{A}}))$ [50]. The experimental Butler–Volmer equation can be deduced from the reaction rate model (2.35) by assuming that all processes, except the intercalation reaction rate, are in thermodynamic equilibrium [50]. However, in the general non-equilibrium setting of intercalation materials, the specific dependency of the reaction rate ${\underset{s}{R}\,\!}$ on (local) concentrations and electrostatic potentials is more difficile. Since the reaction rate couples the adjacent diffusion equations, i.e. ion transport in the electrolyte and intercalated lithium diffusion in the active phase, their concentration dependency also impacts the overall behaviour. From a thermodynamic point of view, it is most importantly that the very same chemical potential function $\mu_{\texttt{A}}$ is employed in the surface reaction rate and in the adjacent diffusion equations.

As shown in the remark of Section 2.1.2, the linear diffusion equation for intercalated lithium the active phase corresponds to an ideal lattice gas, which states $\mu_{\texttt{A}} = g_{{{\texttt{A}}_C}} + k_{\tiny\text{B}}T\, \text{ln}\!\left( {\frac{y_{\texttt{A}}}{1-y_{\texttt{A}}}} \right)$ , used for example in [19]. However, in the corresponding reaction boundary conditions of [19] (Butler–Volmer equation), the chemical potential function (in our notation) $\mu_{\texttt{A}} = g_{{{\texttt{A}}_C}} + frac{k_{\tiny\text{B}}T\,}{e_0}\!\left(\text{ln}\!\left( {\frac{y_{\texttt{A}}}{1-y_{\texttt{A}}}} \right) + \beta y_{\texttt{A}} + \zeta \right)$ , where $\beta$ and $\zeta$ are fit parameters [66], is used and thus reflects the non-ideal OCP behaviour on a solid lattice. This is accompanied with a concentration-dependent exchange current density, which is ${\underset{s}{L}\,\!}$ in our formulation (2.35). Albeit good results on experimental validation, this model is not within our framework as different chemical potential functions $\mu_{\texttt{A}}(y_{\texttt{A}})$ for lithium diffusion in the solid phase and the intercalation reaction at the boundary are used. Reference [50] addresses this aspect more in detail and we emphasise again that the goal of the presented framework is overall thermodynamic consistency for rather general chemical potential functions $\mu_{\texttt{A}}(y_{\texttt{A}})$ .

2.1.4. Balance equations

Applying the above modelling procedure for $j = {\text{An}},{\text{Cat}}$ yields the following equation system,

(2.37) \begin{align} \frac{\partial n_{{{\texttt{A}}_C}}^j}{\partial t} & = -{\text{div}} {\textbf{J}}_{{{\texttt{A}}_C}}^j & \text{with} \quad & {\textbf{J}}_{{{\texttt{A}}_C}}^j = - D_{\texttt{A}}^j \cdot n_{{\texttt{A}},\tiny\text{lat}}^j \, \Gamma_{\texttt{A}}^j \cdot \nabla y_{\texttt{A}}^j && {\textbf{x}} \in \Omega_{\texttt{A}}^j, \end{align}

(2.38) \begin{align} 0 &= -{\text{div}} {\textbf{J}}_{\texttt{S}_{q}}^j & \text{with} \quad & {\textbf{J}}_{{\texttt{S}}_q}^j = - \sigma_{\texttt{S}}^j({\textbf{x}}) \nabla \varphi_{\texttt{S}}^j && {\textbf{x}} \in \Omega_{\texttt{S}}^j, \end{align}

(2.39) \begin{align} \frac{\partial n_{{\texttt{E}}_C}}{\partial t} &= -{\text{div}} {\textbf{J}}_{{\texttt{E}}_C} & \text{with} \quad & {\textbf{J}}_{{\texttt{E}}_C} = - D_{\texttt{E}} \cdot n_{{\texttt{E}},\tiny\text{tot}} \, \Gamma_{\texttt{E}} \cdot \nabla y_{\texttt{E}} + \frac{t_{{\texttt{E}}_C}}{e_0} {\textbf{J}}_{{\texttt{E}}_q} && {\textbf{x}} \in \Omega_{\texttt{E}},\end{align}

(2.40) \begin{align} 0 &= -{\text{div}} {\textbf{J}}_{{\texttt{E}}_q} & \text{with} \quad & {\textbf{J}}_{{\texttt{E}}_q} = - S_{\texttt{E}} \cdot n_{{\texttt{E}},\tiny\text{tot}}\, \Gamma_{\texttt{E}} \nabla y_{\texttt{E}} - \Lambda_{\texttt{E}} n_{{\texttt{E}}_C} \nabla \varphi_{\texttt{E}} && {\textbf{x}} \in \Omega_{\texttt{E}},\end{align}

where (2.37) is the balance of intercalated lithium within the active phase, (2.38) the charge balance in the solid phase and (2.38) $_2$ the electron flux, (2.39) the balance of lithium ions in the electrolyte phase and (2.40) the charge balance in the electrolyte, where (2.40) $_2$ is the flux of the charge $q_{\texttt{E}}$ in the electrolyte. Note that $\sigma_{\texttt{S}}^j({\textbf{x}})$ incorporates the fact that the conductivity is far larger in the conductive additive phase $\Omega_{\texttt{C}}^j$ then in the active particle phase $\Omega_{\texttt{A}}^j$ . The index ${}^{j}$ is necessary because anode and cathode are in general different materials, hence having different parameters and material functions, but (2.37) yields a compact typeface.

The boundary conditions at the interface $\Sigma_{{{\texttt{A}}{\texttt{E}}}}^j = \Omega^j_{\texttt{A}} \cap \Omega^j_{\texttt{E}}$ , where the intercalation reaction ${\text{Li}^+ + \text{e}^- \rightleftharpoons \text{Li}}$ occurs, read

(2.41) \begin{align} {\textbf{J}}_{{{\texttt{A}}_C}}^j \cdot {\textbf{n}} & = -{\underset{s}{L}\,\!}^j \cdot g\left(\frac{1}{k_{\tiny\text{B}}T\,} {\underset{s}{\lambda}\,\!}_{{{\texttt{A}}{\texttt{E}}}}^j \right) && \text{on } \Sigma_{{{\texttt{A}}{\texttt{E}}}}^j , \end{align}

(2.42) \begin{align} {\textbf{J}}_{{\texttt{S}}_q}^j \cdot {\textbf{n}} &= - e_0 {\underset{s}{L}\,\!}^j \cdot g\left(\frac{1}{k_{\tiny\text{B}}T\,} {\underset{s}{\lambda}\,\!}_{{{\texttt{A}}{\texttt{E}}}}^j \right) && \text{on } \Sigma_{{{\texttt{A}}{\texttt{E}}}}^j , \end{align}

(2.43) \begin{align} {\textbf{J}}_{{\texttt{E}}_C} \cdot {\textbf{n}} & = {\underset{s}{L}\,\!}^j\cdot g\left(\frac{1}{k_{\tiny\text{B}}T\,} {\underset{s}{\lambda}\,\!}_{{{\texttt{A}}{\texttt{E}}}}^j \right) && \text{on } \Sigma_{{{\texttt{A}}{\texttt{E}}}}^j , \end{align}

(2.44) \begin{align} {\textbf{J}}_{{\texttt{E}}_q} \cdot {\textbf{n}} &= e_0 {\underset{s}{L}\,\!}^j \cdot g\left(\frac{1}{k_{\tiny\text{B}}T\,} {\underset{s}{\lambda}\,\!}_{{{\texttt{A}}{\texttt{E}}}}^j \right) && \text{on } \Sigma_{{{\texttt{A}}{\texttt{E}}}}^j ,\end{align}

where by convention ${\textbf{n}}$ is pointing from $\Omega_{\texttt{A}}^j$ into $\Omega_{\texttt{E}}^j$ . Note this formulation of the boundary conditions neglects double-layer charging, i.e. currents arising from temporal variations of the boundary layer charge [43, 48, 50, 63], which are rather small compared to the intercalation current [50]. At the interface $\Sigma_{{{\texttt{C}}{\texttt{E}}}}^j = \Omega^j_{\texttt{C}} \cap \Omega^j_{\texttt{E}}$ we have homogenous Neumann boundary conditions, i.e.

(2.45) \begin{align} {\textbf{J}}_{{{\texttt{A}}_C}}^j \cdot {\textbf{n}} = {\textbf{J}}_{{\texttt{S}}_q}^j \cdot {\textbf{n}} = {\textbf{J}}_{{\texttt{E}}_C}\cdot {\textbf{n}} = {\textbf{J}}_{{\texttt{E}}_q} \cdot {\textbf{n}} = 0 \qquad\qquad \text{on } \Sigma_{{{\texttt{C}}{\texttt{E}}}}^j.\end{align}

Further we have global boundary conditions for the electric current density i leaving the anode (discharge, $i>0$ ) or entering the anode (charge, $i<0$ ), i.e.

(2.46) \begin{align} {\textbf{J}}_{{\texttt{S}}_q}^{\text{Cat}} \cdot {\textbf{n}} = - i \qquad\qquad \text{on } \Sigma_{{{\texttt{C}}{\texttt{D}}}}^{\text{Cat}},\end{align}

as well as a reference value for the electrostatic potential, which reads

(2.47) \begin{align} \varphi_S^{\text{An}} = 0 \qquad\qquad \text{on } \Sigma_{{{\texttt{C}}{\texttt{D}}}}^{\text{An}}.\end{align}

2.2. Homogenisation

In this section, the coupled transport equation system (2.37)–(2.47) is spatially homogenised with asymptotic expansion methods [1, 15, 40, 43, 59, 74]. This is exemplarily done for a single porous electrode consisting of spherical particles. Based on an induced scaling the equation system is (i) non-dimensionalised in Section 2.2.2, (ii) homogenised Section 2.2.3, and finally re-dimensionalised in Section 2.2.4. This yields a one to one assignment of the equation system (2.37)–(2.47) to the derived homogenised equations (2.101)–(2.104) and the global boundary conditions. The intermediate part Sections 2.2.1–2.2.3 can be skipped for a better reading flow.

2.2.1. Introduction

An important feature of the coupled transport equation system (2.37)–(2.40) is the circumstance that the solid-state diffusion $D_{\texttt{A}}^j$ is far smaller than the electrolyte diffusivity $D_{\texttt{E}}$ . This accompanies, however, with smaller diffusion pathways on the length scale $\ell$ of the intercalation particles. The diffusivity of Li in LiCoO_2 is, for instance, about $D_{\texttt{A}} \approx 10^{-12} \,/\,{\frac{\text{cm}}{\,\text{s}}}\,$ [83], while the diffusivity of ${Li^+}$ in DMC is in the order of $D_{\texttt{E}} \approx 10^{-5}\,/\,{\frac{\text{cm}}{\,\text{s}}}\,$ [47]. The macro-length scale is $\mathscr{L} \approx 1/{\unicode{x03BC}} {\rm m}$ [19], while the micro-scale is $\ell \approx 1 \,/\,{\text{nm}}\,$ (see Figure 2). Hence, the length scale ration $\frac{\ell}{\mathscr{L}} = \varepsilon \approx 10^{-3}$ and $D_{\texttt{A}}^j \approx \varepsilon^2 D_{\texttt{E}}$ . This motivates the re-scaling

(2.48) \begin{align} D_{\texttt{A}}^j = \varepsilon^2 \cdot D_{\texttt{A}}^{j,\varepsilon} && j = {\text{An}},{\text{Cat}},\end{align}

which is subsequently exploited in the homogenisation procedure.

Figure 2. Sketch of the homogenisation procedure. The porous electrode is simplified as a network of interconnected active phase spheres, yielding a unit cell $\omega$ containing one electrode particle.

2.2.2. Non-dimensionalisation

For the sake of the homogenisation as well as parameter studies and numerical implementations, it is necessary to non-dimensionalise the overall equation system. This is performed in several steps, starting from an initial non-dimensionalisation for the homogenisation, briefly sketched here. Consider the dimensional equation system

(2.49) \begin{align} \frac{\partial w}{\partial t} &= {\text{div}} D \nabla w && {\textbf{x}} \in \Omega_{\texttt{E}} \end{align}

(2.50) \begin{align} D \nabla w \cdot {\textbf{n}} = R \qquad\qquad\qquad\qquad\quad\qquad\qquad {\textbf{x}} \in \Sigma_{{{\texttt{A}}{\texttt{E}}}}.\quad\end{align}

For the sake of the homogenisation, it is convenient to introduce the scaling

(2.51) \begin{align} w &= u \cdot n & \tilde D & = \frac{T}{\mathscr{L}^{\,2}} \cdot D \\[-5pt]\nonumber \end{align}

(2.52) \begin{align} t &= \tau \cdot T & \tilde R &= \frac{T}{n} \frac{1}{\ell} \cdot R \\[-7pt]\nonumber \end{align}

(2.53) \begin{align} x &= \xi \cdot \mathscr{L} & \varepsilon &= \frac{\ell}{\mathscr{L}},\end{align}

which yields

(2.54) \begin{align} \frac{\partial u}{\partial \tau} &= {\text{div}}_\xi \tilde D \nabla_\xi u && \xi \in \tilde \Omega_{\texttt{E}},\\[-5pt]\nonumber \end{align}

(2.55) \begin{align} \tilde D \nabla_\xi u \cdot {\textbf{n}} &= \varepsilon \tilde R && \xi \in \tilde \Sigma_{{{\texttt{A}}{\texttt{E}}}}.\end{align}

2.2.3. General homogenisation framework

For a single electrode, dropping the super-index $^j$ and denoting in this subsection x as non-dimensional space variable, we have a coupled equation system of the form⁷

(2.56) \begin{align} \frac{\partial u_{\texttt{E}}}{\partial \tau} &= {\text{div}} \big(\tilde D_{\texttt{E}} h_{{\texttt{E}}}(u_{\texttt{E}}) \nabla u_{\texttt{E}}\big) && \in \Omega_{\texttt{E}}, \\[-6pt]\nonumber \end{align}

(2.57) \begin{align} \frac{\partial u_{\texttt{A}}}{\partial \tau} &= {\text{div}} \big(\varepsilon^2 \tilde D_{\texttt{A}} h_{{\texttt{A}}}(u_{\texttt{A}}) \nabla u_{\texttt{A}} \big) && \in \Omega_{\texttt{A}}, \\[-6pt]\nonumber \end{align}

(2.58) \begin{align} \frac{\partial u_{\texttt{S}}}{\partial \tau} &= {\text{div}} \big(\tilde D_{\texttt{S}}^\varepsilon h_{{\texttt{S}}}(u_{\texttt{S}},{\textbf{x}}) \nabla u_{\texttt{S}} \big) && \in \Omega_{\texttt{S}}, \\[-6pt]\nonumber \end{align}

(2.59) \begin{align} h_{{\texttt{E}}}(u_{\texttt{E}}) \tilde D_{\texttt{E}} \nabla u_{\texttt{E}} \cdot {\textbf{n}} &= \varepsilon \tilde R_{\texttt{E}} \!\left(u_{\texttt{E}} \big|_{\Sigma_{{{\texttt{A}}{\texttt{E}}}}}^+ ,u_{\texttt{A}} \big|_{\Sigma_{{{\texttt{A}}{\texttt{E}}}}}^-,u_{\texttt{S}} \big|_{\Sigma_{{{\texttt{A}}{\texttt{E}}}}}^-\right) && \text{on } \Sigma_{{{\texttt{A}}{\texttt{E}}}}, \\[-6pt]\nonumber \end{align}

(2.60) \begin{align} h_{{\texttt{A}}}(u_{\texttt{A}}) \varepsilon^2 \tilde D_{\texttt{A}} \nabla u_{\texttt{A}} \cdot {\textbf{n}} &= \varepsilon \tilde R_{\texttt{A}} \!\left(u_{\texttt{E}} \big|_{\Sigma_{{{\texttt{A}}{\texttt{E}}}}}^+ ,u_{\texttt{A}} \big|_{\Sigma_{{{\texttt{A}}{\texttt{E}}}}}^-,u_{\texttt{S}} \big|_{\Sigma_{{{\texttt{A}}{\texttt{E}}}}}^-\right) && \text{on } \Sigma_{{{\texttt{A}}{\texttt{E}}}}, \\[-6pt]\nonumber \end{align}

(2.61) \begin{align} h_{{\texttt{S}}}(u_{\texttt{S}},{\textbf{x}}) \tilde D_{\texttt{S}} \nabla u_{\texttt{S}} \cdot {\textbf{n}} &= \varepsilon \tilde R_{\texttt{S}} \!\left(u_{\texttt{E}} \big|_{\Sigma_{{{\texttt{A}}{\texttt{E}}}}}^+ ,u_{\texttt{A}} \big|_{\Sigma_{{{\texttt{A}}{\texttt{E}}}}}^-,u_{\texttt{S}} \big|_{\Sigma_{{{\texttt{A}}{\texttt{E}}}}}^-\right) && \text{on } \Sigma_{{{\texttt{A}}{\texttt{E}}}}, \\[-6pt]\nonumber \end{align}

(2.62) \begin{align} h_{{\texttt{E}}}(u_{\texttt{E}}) \tilde D_{\texttt{E}} \nabla u_{\texttt{E}} \cdot {\textbf{n}} &= h_{{\texttt{A}}}(u_{\texttt{A}}) \varepsilon^2 \tilde D_{\texttt{A}} \nabla u_{\texttt{A}} \cdot {\textbf{n}} = h_{{\texttt{S}}}(u_{\texttt{S}},{\textbf{x}}) \tilde D_{\texttt{S}} \nabla u_{\texttt{S}} \cdot {\textbf{n}} = 0 && \text{on } \Sigma_{{{\texttt{C}}{\texttt{E}}}},\end{align}

where $u_i, i = {\texttt{A}},{\texttt{E}},{\texttt{S}}$ , denotes the respective phase variable, $(\tilde D_i^\varepsilon, \tilde R_i^\varepsilon)$ the $\varepsilon$ -dependent non-equilibrium parameter and $h_{i}$ capture the non-linearities w.r.t. $u_i$ . For the solid-phase $\Omega_{\texttt{S}}$ , the ${\textbf{x}}$ dependency of $h_{{\texttt{S}}}\!\left(u_{\texttt{S}},\frac{{\textbf{x}}}{\varepsilon}\right)$ accounts for different conductivities in $\Omega_{\texttt{C}} \subset \Omega_{\texttt{S}}$ and $\Omega_{\texttt{A}} \subset \Omega_{\texttt{S}}$ , e.g.

(2.63) \begin{align} h_{\texttt{S}}(u_{\texttt{S}},{\textbf{x}}) = h_{\texttt{S}}^u(u_{\texttt{S}}) \cdot h_{\texttt{S}}^x\!\left(\frac{{\textbf{x}}}{\varepsilon}\right) \quad \text{with} \quad h_{\texttt{S}}^x = \begin{cases} h_{{\texttt{S}}}^{x,{\texttt{A}}}, &\text{ if } {\textbf{x}} \in \Omega_{\texttt{A}}\\[4pt] h_{{\texttt{S}}}^{x,{\texttt{C}}}, &\text{ if } {\textbf{x}} \in \Omega_{\texttt{C}} \end{cases}.\end{align}

Note that we abbreviate

(2.64) \begin{align} u_i \big|_{\Sigma_{{{\texttt{A}}{\texttt{E}}}}}^+ \,=\!:\, u_{\texttt{E}}|^+ \quad \text{and} \quad u_i \big|_{\Sigma_{{{\texttt{A}}{\texttt{E}}}}}^- \,=\!:\, u_i|^-.\end{align}

We consider a multi-scale expansion ( $i={\texttt{A}},{\texttt{E}},{\texttt{S}}$ )

(2.65) \begin{align} u_i({\textbf{x}},t) = \sum_{k=0}^{\infty} \varepsilon^k u_i^k(x,y,t) \quad \text{with} \quad y = \frac{x}{\varepsilon} ,\end{align}

whereby

(2.66) \begin{align} \nabla &= \nabla_x + \varepsilon^{-1} \nabla_y \end{align}

(2.67) \begin{align} {\text{div}} &= {\text{div}}_x + \varepsilon^{-1} {\text{div}}_y. \end{align}

For non-linear functions $h=h(u)$ , we consider the $\varepsilon-$ Taylor expansion

(2.68) \begin{align} h(u) = \sum_{k=0}^{\infty} \frac{1}{k!} \frac{d^k h}{d u^k}(u-u^0) &= h(u^0) + \varepsilon u^1 h'(u^0) + {\mathcal{O}}\!\left(\varepsilon^2\right)\end{align}

and for $g=g\!\left(u_{\texttt{E}}^+,u_{\texttt{A}}^-,u_{\texttt{S}}^-\right)$ a multi-dimensional Taylor expansions, i.e.

(2.69) \begin{align} g\!\left(u_{\texttt{E}}|^+,u_{\texttt{A}}|^-,u_{\texttt{S}}^-\right) = g\!\left( u_{\texttt{E}}^0|^+,u_{\texttt{A}}^0|^-,u_{\texttt{S}}^0|^-\right) + \sum_{i={\texttt{E}},{\texttt{A}},{\texttt{S}}} \varepsilon \cdot \partial_{u_i} g \big|_{u_{\texttt{E}}^0|^+,u_{\texttt{A}}^0|^-,u_{\texttt{S}}^0|^-} \cdot u_i^1|^\pm + {\mathcal{O}}\!\left(\varepsilon^2\right).\end{align}

We abbreviate

(2.70) \begin{align} &h^0 \,:\!=\, h(u^0)\, , h^1 \,:\!=\, u_1 \cdot h'(u^0) , \quad \text{and} \quad g^0 \,:\!=\, g\!\left( u_{\texttt{E}}^0|^+,u_{\texttt{A}}^0|^-,u_{\texttt{S}}^0|^-\right) \end{align}

(2.71) \begin{align} &g^1 \,:\!=\, u_{\texttt{E}}^1 \partial_{u_{\texttt{E}}} g \!\left(u_{\texttt{E}}^0|^+,u_{\texttt{A}}^0|^-,u_{\texttt{S}}^0|^-\right) + u_{\texttt{A}}^1 \partial_{u_{\texttt{A}}} g \!\left(u_{\texttt{E}}^0|^+,u_{\texttt{A}}^0|^-,u_{\texttt{S}}^0|^-\right) + u_{\texttt{S}}^1 \partial_{u_{\texttt{S}}} g \!\left(u_{\texttt{E}}^0|^+,u_{\texttt{A}}^0|^-,u_{\texttt{S}}^0|^-\right)\end{align}

and expand thus

(2.72) \begin{align} h_i &= h_i^0 + \varepsilon h_i^1 \qquad\qquad g = g^0 + \varepsilon g^1.\end{align}

Insertion of the multi-scale expansion yields a sequence of PDEs in the orders of $\varepsilon$ .

Order $\varepsilon^{-2}$ :

For $i={\texttt{E}},{\texttt{S}}$ we have

(2.73) \begin{align} {\text{div}}_y \tilde D_i h_{i}^0 \nabla_y u_i^0 = 0 \quad \text{in} \quad \omega_i\end{align}

with

(2.74) \begin{align} \tilde D_i h_{i}^0 \nabla_y u_i^0 \cdot {\textbf{n}} = 0 \quad \text{on} \quad \sigma_{{{\texttt{A}}{\texttt{E}}}}.\end{align}

This yields $u_i^0 = u_i^0({\textbf{x}},t)$ by the maximum principle.

Order $\varepsilon^{-1}$ :

For $i={\texttt{E}},{\texttt{S}}$ , we have due to $u_i^0 = u_i^0({\textbf{x}},t)$ we have

(2.75) \begin{align}{\text{div}}_y \tilde D_i h_{i}^0 \nabla_y u_i^1 = 0 \quad \text{in} \quad \omega_i\end{align}

with

(2.76) \begin{align} \left(\tilde D_i h_{i}^0 \nabla_x u_i^0 + \tilde D_i h_{i}^0 \nabla_y u_i^1\right) \cdot {\textbf{n}} = 0 \quad \text{on} \quad \sigma_{{{\texttt{A}}{\texttt{E}}}}.\end{align}

This yields

(2.77) \begin{align} u_i^1(x,y,t) = - \nabla_x u_i^0 \cdot \chi_i(y),\end{align}

where ${\boldsymbol{\chi}}_{\texttt{E}} = \left(\chi_{\texttt{E}}^1,\chi_{\texttt{E}}^2,\chi_{\texttt{E}}^3\right)$ satisfies the cell problem ( $k=1,2,3$ )

(2.78) \begin{align}(\textbf{CP}_{\texttt{E}}) \qquad \begin{cases} {\text{div}}_y \nabla_y \chi_{\texttt{E}}^k & = 0 \qquad {\textbf{y}} \in \omega_{\texttt{E}}, \, i ={\texttt{E}},{\texttt{S}} \\[4pt] \nabla \chi_{\texttt{E}}^k \cdot {\textbf{n}} & = n_k \qquad\ \ \text{on } \sigma_{{{\texttt{A}}{\texttt{E}}}} \\[4pt] \quad \chi_{\texttt{E}}^k & \text{ is periodic }\end{cases}\end{align}

and ${\boldsymbol{\chi}}_{\texttt{S}} = (\chi_{\texttt{S}}^1,\chi_{\texttt{S}}^2,\chi_{\texttt{S}}^3)$ satisfies the cell problem ( $k=1,2,3$ )

(2.79) \begin{align}(\textbf{CP}_{\texttt{S}}) \qquad \begin{cases} {\text{div}}_y \!\left( h_{{\texttt{S}}}^x({\textbf{y}}) \cdot \left( {\textbf{e}}^k - \nabla_y \chi_{\texttt{S}}^k \right) \right) & = 0 \qquad {\textbf{y}} \in \omega_{\texttt{S}} \\[4pt] h_{{\texttt{S}}}^x({\textbf{y}}) \cdot \left( {\textbf{e}}^k - \nabla_y \chi_{\texttt{S}}^k \right) \cdot {\textbf{n}} & = 0 \qquad \ \text{on } \sigma_{{{\texttt{A}}{\texttt{E}}}} \\[4pt] \quad \chi_{\texttt{S}}^k & \text{ is periodic}.\end{cases}\end{align}

Note that $n_k$ denotes the kth component of the corresponding normal vector ${\textbf{n}}$ on the surface $\sigma_{{{\texttt{A}}{\texttt{E}}}}$ .

Order $\varepsilon^{0}$

Since $u_{\texttt{E}}^0 = u_{\texttt{E}}^0({\textbf{x}},t)$ we have for $i={\texttt{E}},{\texttt{S}}$

(2.80) \begin{align} {\text{div}}_x \tilde D_i h_{i}^0 \nabla_x u_i^0 + {\text{div}}_x \tilde D_i h_{i}^0 \nabla_y u_i^1 + {\text{div}}_y \tilde D_i h_{i}^0 \nabla_x u_i^1 + {\text{div}}_y \tilde D_i h_{i}^1 \nabla_y u_i^1 = \frac{\partial u_i^0}{\partial \tau}\end{align}

and

(2.81) \begin{align} {\text{div}}_y \tilde D_{\texttt{A}} h_{{\texttt{A}}}^0 \nabla_y u_{\texttt{A}}^0 = \frac{\partial u_{\texttt{A}}^0}{\partial \tau}\end{align}

with

(2.82) \begin{align}\Big(\tilde D_{\texttt{E}} h_{{\texttt{E}}}^0 \nabla_x u_{\texttt{E}}^1 + \tilde D_{\texttt{E}} h_{{\texttt{E}}}^1 \nabla_x u_{\texttt{E}}^0 + \tilde D_{\texttt{E}} h_{{\texttt{E}}}^1 \nabla_y u_{\texttt{E}}^1 \Big) \cdot {\textbf{n}}& = \tilde R_{\texttt{E}}^0\end{align}

(2.83) \begin{align}\Big(\tilde D_{\texttt{S}} h_{{\texttt{S}}}^0 \nabla_x u_{\texttt{S}}^1 + \tilde D_{\texttt{S}} h_{{\texttt{S}}}^1 \nabla_x u_{\texttt{S}}^0 + \tilde D_{\texttt{S}} h_{{\texttt{S}}}^1 \nabla_y u_{\texttt{S}}^1 \Big) \cdot {\textbf{n}}& = \tilde R_{\texttt{S}}^0\end{align}

(2.84) \begin{align}\Big(\tilde D_{\texttt{A}} h_{{\texttt{A}}}^0 \nabla_y u_{\texttt{A}}^0 \Big) \cdot {\textbf{n}} & = \tilde R_{\texttt{A}}^0.\end{align}

Equation (2.80) leads (by an integration over $\omega_{\texttt{E}}$ and integration by parts) for $i = {\texttt{E}}$ to

(2.85) \begin{align} \frac{\partial u_{\texttt{E}}^0}{\partial \tau} & = {\text{div}}_x \big( \tilde D_{\texttt{E}} h_{{\texttt{E}}}^0 \pi_{\texttt{E}} \cdot \nabla_x u_{\texttt{E}}^0 \big) + \frac{a_{{{\texttt{A}}{\texttt{E}}}}}{\psi_{\texttt{E}} } \cdot \frac{1}{{\text{area}\{{\sigma_{{{\texttt{A}}{\texttt{E}}}}}\}} } \int_{\sigma_{{{\texttt{A}}{\texttt{E}}}}} \tilde R_{\texttt{E}}^0\!\left(u_{\texttt{E}}^0, u_{\texttt{S}}^0,u_{\texttt{A}}^0|_{\partial \omega_{\texttt{A}}}\right) dA({\textbf{y}})\end{align}

with

(2.86) \begin{align} {\boldsymbol{\pi}}_{\texttt{E}} &\,:\!=\, \frac{1}{{\text{vol}\{{\omega_{\texttt{E}}}\}} } \int_{\omega_{\texttt{E}}} ({ \textbf{Id}} - \boldsymbol{\nabla}_y {\boldsymbol{\chi}}_{\texttt{E}} ) dV(y) = { \textbf{Id}} - \frac{1}{{\text{vol}\{{\omega_{\texttt{E}}}\}} } \int_{\omega_{\texttt{E}}}\boldsymbol{\nabla}_y {\boldsymbol{\chi}}_{\texttt{E}} dV(y)\\[-4pt]\nonumber \end{align}

(2.87) \begin{align} \psi_{\texttt{E}} &= \frac{\int_{\omega_{\texttt{E}}} 1 \,dV }{\int_{\omega} 1 \,dV } = \frac{{\text{vol}\{{\omega_{\texttt{E}}}\}}}{{\text{vol}\{{\omega}\}} } \\[-4pt]\nonumber \end{align}

(2.88) \begin{align} a_{{{\texttt{A}}{\texttt{E}}}} &= \frac{\int_{\sigma_{{{\texttt{A}}{\texttt{E}}}}} 1\,dA}{{\int_{\omega} 1 \,dV } } = \frac{{\text{area}\{{\sigma_{{{\texttt{A}}{\texttt{E}}}}}\}}}{{\text{vol}\{{\omega}\}}}.\end{align}

Note that $u_{\texttt{A}}^0=u_{\texttt{A}}^0({\textbf{x}},{\textbf{y}},t)$ and that $u_{\texttt{A}}^0|_{\partial \omega_{\texttt{A}}}$ denotes an evaluation of $u_{\texttt{A}}^0$ at the boundary of the micro-domain $\omega_{\texttt{A}}$ . Equation (2.80) leads (by an integration over $\omega_{\texttt{S}}$ ) for $i = {\texttt{S}}$ to

(2.89) \begin{align}\frac{\partial u_{\texttt{S}}^0}{\partial \tau} & = {\text{div}}_x \big( \tilde D_{\texttt{S}}^0 h_{{\texttt{S}}}^0 {\boldsymbol{\pi}}_{\texttt{S}} \cdot \nabla_x u_{\texttt{S}}^0 \big) + \frac{a_{{{\texttt{A}}{\texttt{E}}}}}{\psi_{\texttt{S}} } \frac{1}{{\text{area}\{{\sigma_{{{\texttt{A}}{\texttt{E}}}}}\}} } \int_{\sigma_{{{\texttt{A}}{\texttt{E}}}}} \tilde R_{\texttt{S}}^0\!\left(u_{\texttt{E}}^0, u_{\texttt{S}}^0,u_{\texttt{A}}^0|_{\partial \omega_{\texttt{A}}}\right) dA({\textbf{y}})\end{align}

with

(2.90) \begin{align} {\boldsymbol{\pi}}_{\texttt{S}} &\,:\!=\, \frac{1}{{\text{vol}\{{\omega_{\texttt{S}}}\}} } \int_{\omega_{\texttt{S}}} h_{{\texttt{S}},3} (y) ({ \textbf{Id}} - \boldsymbol{\nabla}_y {\boldsymbol{\chi}}_{\texttt{S}} ) dV(y) \\[-4pt]\nonumber \end{align}

(2.91) \begin{align} \psi_{\texttt{S}} &= \frac{\int\limits_{\omega_{\texttt{S}}} 1 \,dV }{\int\limits_{\omega} 1 \,dV } = \frac{{\text{vol}\{{\omega_{\texttt{S}}}\}}}{{\text{vol}\{{\omega}\}} }\end{align}

Hence we obtain for the equation system (2.56)–(2.62) via periodic homogenisation the coupled micro-macro balance equation system

(2.92) \begin{align}\frac{\partial u_{\texttt{E}}^0({\textbf{x}},t)}{\partial \tau} & = {\text{div}}_x \big( \tilde D_{\texttt{E}} h_{{\texttt{E}}}^0 {\boldsymbol{\pi}}_{\texttt{E}} \cdot \nabla_x u_{\texttt{E}}^0 \big) + \frac{a_{{{\texttt{A}}{\texttt{E}}}}}{\psi_{\texttt{E}}} \frac{1}{{\text{area}\{{\sigma_{{{\texttt{A}}{\texttt{E}}}}}\}} } \int_{\sigma_{{{\texttt{A}}{\texttt{E}}}}} \tilde R_{\texttt{E}}^0(u_{\texttt{E}}^0, u_{\texttt{S}}^0,u_{\texttt{A}}^0|_{\partial \omega_{\texttt{A}}}) dA({\textbf{y}}) \end{align}

(2.93) \begin{align}\frac{\partial u_{\texttt{S}}^0({\textbf{x}},t)}{\partial \tau} & = {\text{div}}_x \big( \tilde D_{\texttt{S}}^0 h_{{\texttt{S}}}^0 {\boldsymbol{\pi}}_{\texttt{S}} \cdot \nabla_x u_{\texttt{S}}^0 \big) + \frac{a_{{{\texttt{A}}{\texttt{E}}}}}{\psi_{\texttt{S}}} \frac1{{\text{area}\{{\sigma_{{{\texttt{A}}{\texttt{E}}}}}\}} } \int_{\sigma_{{{\texttt{A}}{\texttt{E}}}}} \tilde R_{\texttt{S}}^0(u_{\texttt{E}}^0, u_{\texttt{S}}^0,u_{\texttt{A}}|_{\partial \omega_{\texttt{A}}}) dA({\textbf{y}}) \end{align}

(2.94) \begin{align}\frac{\partial u_{\texttt{A}}^0(x,y,t)}{\partial \tau} & = {\text{div}}_y \tilde D_{\texttt{A}} h_{{\texttt{A}}}^0 \nabla_y u_{\texttt{A}}^0\end{align}

with

(2.95) \begin{align} \Big(\tilde D_{\texttt{A}} h_{{\texttt{A}}}^0 \nabla_y u_{\texttt{A}}^0 \Big) \cdot {\textbf{n}} = \tilde R_{\texttt{A}}^0.\end{align}

Spherical particles

In the following, we assume on the macro-scale a 1D approximation $x \in [0,1]$ as well as spherical particles $\omega_{\texttt{A}}$ , i.e. $\omega = [{-}0.5,0.5]^3$ ( ${\text{vol}\{{\omega}\}} = 1$ ), $\omega_{\texttt{A}} = B_{\tilde r_{\texttt{A}}}(\textbf{0})$ , $ \tilde r_{\texttt{A}} < 0.5$ , $\omega_{\texttt{E}} = \omega \backslash \omega_{\texttt{A}}$ , yielding

(2.96) \begin{align}\psi_{\texttt{E}} \frac{\partial u_{\texttt{E}}^0(x,t)}{\partial \tau} & = \partial_x \big( \psi_{\texttt{E}} \tilde D_{\texttt{E}} h_{{\texttt{E}}}^0 {\boldsymbol{\pi}}_{\texttt{E}} \cdot \partial_x u_{\texttt{E}}^0 \big)\!+\!\theta_{{{\texttt{A}}{\texttt{E}}}} \tilde R_{\texttt{E}}^0\!\left(u_{\texttt{E}}^0,u_{\texttt{S}}^0,u_{\texttt{A}}^0\!\left(x,\tilde r_{\texttt{A}},t\right) \right) , && x \in [0,1]\\[-4pt]\nonumber\end{align}

(2.97) \begin{align}\psi_{\texttt{S}} \frac{\partial u_{\texttt{S}}^0(x,t)}{\partial \tau} & = \partial_x \big( \psi_{\texttt{S}} \tilde D_{\texttt{S}}^0 h_{{\texttt{S}}}^0 {\boldsymbol{\pi}}_{\texttt{S}} \cdot \partial_x u_{\texttt{S}}^0 \big)\!+\!\theta_{{{\texttt{A}}{\texttt{E}}}} \tilde R_{\texttt{S}}^0\!\left(u_{\texttt{E}}^0,u_{\texttt{S}}^0,u_{\texttt{A}}^0\!\left(x,\tilde r_{\texttt{A}},t\right) \right) , && x \in [0,1]\\[-4pt]\nonumber\end{align}

(2.98) \begin{align}\frac{\partial u_{\texttt{A}}^0(x,r,t)}{\partial \tau} & = \frac{1}{r^2} \partial_r \big( r^2 \tilde D_{\texttt{A}} h_{{\texttt{A}}}^0 \partial_r u_{\texttt{A}}^0 \big) , && x \in [0,1], r \in (0,\tilde r_{\texttt{A}})\\[-4pt]\nonumber\end{align}

with

(2.99) \begin{align} \Big( \tilde D_{\texttt{A}} h_{{\texttt{A}}}^0 \partial_r u_{\texttt{A}}^0 \Big) \Big|_{r=r_{\texttt{A}}} & = \tilde R_{\texttt{A}}^0 \!\left(u_{\texttt{E}}^0(x,t),u_{\texttt{S}}^0(x,t),u_{\texttt{A}}^0\!\left(x,\tilde r_{\texttt{A}},t\right)\right)\end{align}

and

(2.100) \begin{align} \theta_{{{\texttt{A}}{\texttt{E}}}} = \frac{{\text{area}\{{\sigma_{{{\texttt{A}}{\texttt{E}}}}}\}}}{{\text{vol}\{{\omega}\}}} = 4 \pi \tilde r_{\texttt{A}}^2 - {\text{area}\{{\sigma_{{\texttt{A}},{\texttt{C}}}}\}}.\end{align}

For spherical particles, various possibilities regarding their micro-structural packing arise, most prominent (i) simple cubic, (ii) body-centered cubic and (iii) face-centered cubic, see Figure 3. For a given micro-structure (or periodic unit cell $\omega$ ), the porous media parameters $(\psi_{\texttt{E}},{\boldsymbol{\pi}}_{\texttt{E}},\theta_{{{\texttt{A}}{\texttt{E}}}})$ and $({\boldsymbol{\pi}}_{\texttt{S}},\psi_{\texttt{S}})$ can (numerically) be computed by solving the cell problems (2.78) and (2.79), respectively. Treating the particle radii $r_{\texttt{A}}$ as a parameter yields then the diffusion corrector ${\boldsymbol{\pi}}_{\texttt{E}}$ as function of the porosity $\psi_{\texttt{E}}$ [54] (solid lines in Figure 3). Quite commonly, the (scalar) tortuosity corrector $\tau_{\texttt{E}}$ is introduced via an effective diffusion coefficient [14, 29, 66, 80], which is in our notation $\tau_{\texttt{E}} = (\pi_{\texttt{E}})^{-1}$ [54]. The Bruggeman approximation states that $\tau_{\texttt{E}} = \psi_{\texttt{E}}^{-\alpha}$ [54], where $\alpha$ is to be determined (Bruggeman fit, dashed lines in Figure 3).

Figure 3. Porous media parameters for simple cubic, body-centered cubic and face centered cubic micro-structures (Figure  15 of [54], reprinted with permission of Elsevier).

Note that for equally sized spherical particles the actual micro-structure is exclusively encoded in the porous media parameters $(\psi_{\texttt{E}},{\boldsymbol{\pi}}_{\texttt{E}},\theta_{{{\texttt{A}}{\texttt{E}}}},{\boldsymbol{\pi}}_{\texttt{S}},\psi_{\texttt{S}})$ of the homogenised equation system (2.96)–(2.98). We refer to [54] for more complex micro-structure geometries as well as their meshing and numerical solution of (2.78) and proceed for the sake of this work with a simple cubic micro-structure.

2.2.4. Homogenised battery model

Applying this homogenisation scheme to the equation system (2.37)–(2.40) and its boundary conditions (2.41)–(2.45), dropping the leading order index $^0$ and reinserting the scalings of Sections 2.2.1 and 2.2.2 yields the following equation system ( $j = {\text{An}},{\text{Cat}}$ ):

(2.101) \begin{align} \frac{\partial n_{{{\texttt{A}}_C}}^j}{\partial t} & = - \frac{1}{r^2} \partial_r J_{{{\texttt{A}}_C}}^j && (r,x) \in [0,r_{\texttt{A}}^j] \times I^j, \nonumber \\ & \text{with} \quad J_{{{\texttt{A}}_C}}^j = - D_{\texttt{A}}^j \cdot n_{{\texttt{A}},\tiny\text{lat}}^j \, \Gamma_{\texttt{A}}^j \cdot r^2 \, \partial_r y_{\texttt{A}}^j \end{align}

(2.102) \begin{align} 0 &= -\partial_x J_{{\texttt{S}}_q}^j - e_0 \theta_{{{\texttt{A}}{\texttt{E}}}} \frac{1}{\ell} {\underset{s}{L}\,\!}^j \cdot g && x \in I^j, \nonumber \\ & \text{with} \quad J_{{\texttt{S}}_q}^j = - \psi_{\texttt{S}} {\boldsymbol{\pi}}_{{\texttt{S}}}^\sigma \sigma_{\texttt{C}}^j \partial_x \varphi_{\texttt{S}}^j \end{align}

(2.103) \begin{align} \psi_{\texttt{E}} \frac{\partial n_{{\texttt{E}}_C}}{\partial t} &= -\partial_x J_{\texttt{E}} + \theta_{{{\texttt{A}}{\texttt{E}}}} \frac{1}{\ell} {\underset{s}{L}\,\!}^j \cdot g && x \in I,\nonumber \\ & \text{with} \quad J_{{\texttt{E}}_C} = - \psi_{\texttt{E}} {\boldsymbol{\pi}}_{\texttt{E}} D_{\texttt{E}} \cdot n_{{\texttt{E}},\tiny\text{tot}} \, \Gamma_{\texttt{E}} \cdot \partial_x y_{\texttt{E}} + \frac{t_{{\texttt{E}}_C}}{e_0} J_{{\texttt{E}}_q} \end{align}

(2.104) \begin{align} 0 &= -\partial_x J_{{\texttt{E}}_q} + e_0 \theta_{{{\texttt{A}}{\texttt{E}}}} \frac{1}{\ell} {\underset{s}{L}\,\!}^j\cdot g && x \in I, \nonumber \\ & \text{with} \quad J_{{\texttt{E}}_q} = - \psi_{\texttt{E}} {\boldsymbol{\pi}}_{\texttt{E}} S_{\texttt{E}} \cdot n_{{\texttt{E}},\tiny\text{tot}}\, \Gamma_{\texttt{E}} \nabla y_{\texttt{E}} - \psi_{\texttt{E}} {\boldsymbol{\pi}}_{\texttt{E}} \Lambda_{\texttt{E}} n_{{\texttt{E}}_C} \nabla \varphi_{\texttt{E}},\end{align}

where

(2.105) \begin{align} I^{\text{An}} &= \left[0,W^{\text{An}}\right], I^{\text{Sep}} = \left[W^{\text{An}},W^{\text{An}}\!+\!W^{\text{Sep}}\right], I^{\text{Cat}} = \left[W^{\text{An}}\!+\!W^{\text{Sep}},W\right],\\[-4pt]\nonumber \end{align}

(2.106) \begin{align} W &= W^{\text{An}} + W^{\text{Sep}} + W^{\text{Cat}}, I = I^{\text{An}} \cup I^{\text{Sep}} \cup I^{\text{Cat}} = [0,W],\\[-4pt]\nonumber \end{align}

(2.107) \begin{align} g &= g\!\left(y_{\texttt{E}}(x,t),\varphi_{\texttt{E}}(x,t),\varphi_S(x,t),y_{\texttt{A}}(x,r_{\texttt{A}}^j,t)\right),\end{align}

and boundary conditions

(2.108) \begin{align} J_{\texttt{A}}^j \Big|_{r=r_{\texttt{A}}} \!\!= - r_{\texttt{A}}^2 {\underset{s}{L}\,\!}^j \cdot g, \quad J_{\texttt{A}}^j \Big|_{r=0} \!\!= 0, \quad J_{{\texttt{S}}_q}^{\text{Cat}}\Big|_{x = W} \!\!= - i, \quad \varphi_S^{\text{An}} \Big|_{x = 0} = 0 .\end{align}

Note that $t_{{\texttt{E}}_C} = \text{const.}$ allows us to rewrite (2.103) with (2.104) as

(2.109) \begin{align}\psi_{\texttt{E}} \frac{\partial n_{{\texttt{E}}_C}}{\partial t} &= -\partial_x J_{\texttt{E}} + \theta_{{{\texttt{A}}{\texttt{E}}}} \frac{1}{\ell} (1-t_{{\texttt{E}}_C}) {\underset{s}{L}\,\!}^j \cdot g && x \in I,\nonumber \\ & \text{with} \quad J_{{\texttt{E}}_C} = - \psi_{\texttt{E}} {\boldsymbol{\pi}}_{\texttt{E}} D_{\texttt{E}} \cdot n_{{\texttt{E}},\tiny\text{tot}} \, \Gamma_{\texttt{E}} \cdot \partial_x y_{\texttt{E}} ,\end{align}

which is further used.

Cell Voltage, Capacity, C-Rate

To introduce proper scalings an non-dimensionalisations, some definitions of important (global) quantities are required. For a LIB, the following quantities are of most importance:

1. (1). the cell voltage

   (2.110) \begin{align} \qquad\qquad\qquad\qquad\qquad\qquad E &\,:\!=\, \varphi_{\texttt{S}}^{\text{An}}|_{x=0} - \varphi_{\texttt{S}}^{\text{Cat}}|_{x=W} \qquad\qquad\qquad\qquad\qquad\left[ {\text{V}} \right], \end{align}

2. (2). the basic (electrode) capacity ( $j={\text{An}},{\text{Cat}}$ )

   (2.111) \begin{align} \qquad\qquad Q^{j,0} &\,:\!=\, \int_{I^j} \frac{4 \pi}{\ell^3} \int_0^{r_{\texttt{A}}^j} e_0 n_{{\texttt{A}},\tiny\text{lat}}^j r^2 \, dr dx = W^j \psi_{\texttt{A}}^j q_{\texttt{A}}^j = \text{const.} && \qquad\qquad\left[ {\text{Ah} \, \text{m}{^2}} \right], \end{align}
   with $q_{\texttt{A}}^j \,:\!=\, e_0 n_{{\texttt{A}},\tiny\text{lat}}$ ,

3. (3). the present (electrode) capacity ( $j={\text{An}},{\text{Cat}}$ )

   (2.112) \begin{align} \qquad\qquad\qquad\qquad Q^{j}(t) &\,:\!=\, \int_{I^j} \frac{4 \pi}{\ell^3} \int_0^{r_{\texttt{A}}^j} e_0 n_{{{\texttt{A}}_C}}^j r^2 \, dr dx && \qquad\qquad\qquad\quad\left[ {\text{Ah} \, \text{m}{^2}} \right], \end{align}

4. (4). the (electrode) status of charge

   (2.113) \begin{align} \qquad\qquad\qquad\qquad \bar y_{\texttt{A}}^j(t) & \,:\!=\, \frac{Q^{j}}{Q^{j,0}} = \frac{1}{W^j} \int_{I^j} \frac{3}{r_{\texttt{A}}^3} \int_0^{r_{\texttt{A}}^j} y_{\texttt{A}}^j \cdot r^2 \, dr dx && \qquad\qquad\quad\left[ {1} \right], \end{align}

5. (5). and the 1-C current density

   (2.114) \begin{align} \qquad\qquad\qquad\qquad\qquad\qquad\qquad i_C & \,:\!=\, \frac{Q^{{\text{Cat}},0}}{1 \!\left[ {\text{h}} \right]} && \qquad\qquad\qquad\qquad\qquad\left[ {\text{A} \text{m}{^2}} \right], \end{align}
   which yields for the current density the scaling
   (2.115) \begin{align} i = C_h \cdot i_C, \end{align}
   where $C_h \in {\mathbb{R}}$ is the C-Rate.

Note that

(2.116) \begin{align} Q^{\text{Cat}}(t)\!-\!Q^{{\text{Cat}},0} & = \int_{I^{\text{Cat}}} \frac{e_0 4 \pi}{(\ell^{\text{Cat}})^3} \int_0^{r_{\texttt{A}}^{\text{Cat}}} \int_0^t r^2 \partial_ t n_{{{\texttt{A}}_C}}^{\text{Cat}} \, dt \, dr \, dx = \int_0^t i \, dt .\end{align}

For a constant current discharge ( $i = \text{const.}$ ), we obtain thus

(2.117) \begin{align} \bar y_{\texttt{A}}^{\text{Cat}}(t) - \bar y_{\texttt{A}}^{\text{Cat}}(t=0) = \frac{C_h}{1 \!\left[ {\text{h}} \right]} \cdot t,\end{align}

which introduces the time re-scaling $\tau = \frac{C_h}{1 \!\left[ {\text{h}} \right]} \cdot t$ .

Scalings

Summarised, we use the scalings

(2.118) \begin{align} t &= \frac{1 \!\left[ {\text{h}} \right]}{C_h} \tau \quad,\; \tau \in [0,1] & i_C & = \frac{Q^{{\text{Cat}},0}}{1 \!\left[ {\text{h}} \right]} = \frac{W^{\text{Cat}} \psi_{\texttt{A}}^{\text{Cat}} q_{\texttt{A}}^{\text{Cat}}}{1 \!\left[ {\text{h}} \right]}\\[-6pt]\nonumber \end{align}

(2.119) \begin{align} n_{{\texttt{E}},\tiny\text{tot}} &= \tilde n_{{\texttt{E}},\tiny\text{tot}} n_{{\texttt{E}},\tiny\text{ref}} & x &= \xi W \quad,\;\xi \in [0,1] \\[-6pt]\nonumber \end{align}

(2.120) \begin{align}n_{{{\texttt{A}}_C}}^j &= n_{{\texttt{A}},\tiny\text{lat}}^j y_{\texttt{A}}^j ,\; y_{\texttt{A}}^j\in[0,1] & \varphi &= \frac{k_{\tiny\text{B}}T\,}{e_0} \tilde \varphi \\[-6pt]\nonumber \end{align}

(2.121) \begin{align} r &= r_{\texttt{A}} \nu \quad,\; \nu \in [0,1] & c_{\texttt{E}} &= \frac{\partial n_{{\texttt{E}}_C}(y_{\texttt{E}})}{\partial y_{\texttt{E}}} \frac{1}{n_{{\texttt{E}},\tiny\text{ref}}} \\[-6pt]\nonumber \end{align}

(2.122) \begin{align}n_{{\texttt{E}}_C} &= \tilde n_{{\texttt{E}},\tiny\text{tot}} \, n_{{\texttt{E}},\tiny\text{ref}} \, y_{\texttt{E}} \quad,\; y_{\texttt{E}} \in [0,1] & \eta_n^{{\texttt{E}},j} &\,:\!=\, \frac{n_{{\texttt{A}},\tiny\text{lat}}^j}{n_{{\texttt{E}},\tiny\text{ref}}^j} \\[-6pt]\nonumber \end{align}

(2.123) \begin{align}\eta_x^j & = \frac{W^j}{W} & & \\[-6pt]\nonumber \end{align}

(2.124) \begin{align}\tilde r_{\texttt{A}}^j &\,:\!=\, \frac{\ell^j}{r_{\texttt{A}}^j} & \eta_W^{\text{Cat}} &\,:\!=\, \frac{\psi_{\texttt{A}}^{\text{Cat}} \cdot W^{\text{Cat}}}{W} \\[-6pt]\nonumber \end{align}

(2.125) \begin{align} \tilde \sigma_{\texttt{S}}^j &\,:\!=\, \frac{1 \!\left[ {\text{h}} \right]}{W^2} \frac{\frac{k_{\tiny\text{B}}T\,}{(e_0)^2}}{n_{{\texttt{A}},\tiny\text{lat}}^j } \cdot \sigma_{\texttt{S}}^j& \tilde \Lambda_{\texttt{E}} &\,=\!:\, \frac{1 \!\left[ {\text{h}} \right]}{W^2} \frac{k_{\tiny\text{B}}T\,}{e_0^2} \Lambda_{\texttt{E}} \\[-6pt]\nonumber \end{align}

(2.126) \begin{align} \tilde D_{\texttt{E}} &\,:\!=\, \frac{1\!\left[ {\text{h}} \right] }{W^2} D_{\texttt{E}} & \tilde S_{\texttt{E}} &\,:\!=\, (2 t_C - 1) \tilde \Lambda_{\texttt{E}} \\[-6pt]\nonumber \end{align}

(2.127) \begin{align} \underset{s}{\tilde{L}} {}^j &\,:\!=\, \frac{1}{n_{{\texttt{A}},\tiny\text{lat}}^j} \big( {1 \!\left[ {\text{h}} \right]} \big) \frac{1}{\ell^j} {\underset{s}{L}\,\!}^j & \tilde D_{{\texttt{A}},\tiny\text{ref}}^j & \,:\!=\, \frac{1 \!\left[ {\text{h}} \right]}{(\ell^j)^2} \cdot D_{{\texttt{A}},\tiny\text{ref}}^j\end{align}

and definitions

(2.128) \begin{align} \Omega^{\text{Cat}} &= \left[0,\eta_x^{\text{Cat}}\right], \Omega^{\text{Sep}} = \left[\eta_x^{\text{Cat}},\eta_x^{\text{Cat}}\!+\!\eta_x^{\text{Sep}}\right], \Omega^{\text{An}} = \left[\eta_x^{\text{Cat}}\!+\!\eta_x^{\text{Sep}},1\right], \Omega^{\text{El}} = \Omega^{\text{Cat}} \cup \Omega^{\text{An}} \end{align}

(2.129) \begin{align} \Omega_{\texttt{A}}^{\text{Cat}} & = \left[0,\eta_x^{\text{Cat}}\right] \times [0,1], \Omega_{\texttt{A}}^{\text{An}} = \left[\eta_x^{\text{Cat}}+\eta_x^{\text{Sep}},1\right] \times [0,1] \quad \Omega_{\texttt{A}} = \Omega_{\texttt{A}}^{\text{Cat}} \cup \Omega_{\texttt{A}}^{\text{An}}\end{align}

as well as

(2.130) \begin{align} y_{\texttt{A}} &\,:\!=\, \begin{cases} y_{\texttt{A}}^{\text{Cat}} & (\xi,r) \in \Omega_{\texttt{A}}^{\text{Cat}} \\[4pt] y_{\texttt{A}}^{\text{An}} & (\xi,r) \in \Omega_{\texttt{A}}^{\text{An}} \\[4pt] \end{cases} &&& \hat \sigma_{\texttt{S}} &\,:\!=\, \begin{cases} \psi_{\texttt{S}}^{\text{Cat}} {\boldsymbol{\pi}}_{{\texttt{S}},\sigma}^{\text{Cat}} \tilde \sigma_{\texttt{C}}^{\text{Cat}} & \xi \in \Omega^{\text{Cat}} \\[4pt] \psi_{\texttt{S}}^{\text{An}} {\boldsymbol{\pi}}_{{\texttt{S}},\sigma}^{\text{An}} \tilde \sigma_{\texttt{C}}^{\text{An}} & \xi \in \Omega^{\text{An}} \\[4pt] \end{cases}\\[-7pt]\nonumber\end{align}

(2.131) \begin{align} \varphi_{\texttt{S}} &\,:\!=\, \begin{cases} \varphi_{\texttt{S}}^{\text{Cat}} & \xi \in \Omega^{\text{Cat}} \\[4pt] \varphi_{\texttt{S}}^{\text{An}} & \xi \in \Omega^{\text{An}} \\[4pt] \end{cases} &&& \hat D_{\texttt{A}} &\,:\!=\, \begin{cases} \tilde D_{{\texttt{A}},\tiny\text{ref}}^{\text{Cat}} \, (1-y_{\texttt{A}}) \Gamma_{\texttt{A}}^{\text{Cat}}(y_{\texttt{A}}) \cdot \nu^2 & (\xi,r) \in \Omega_{\texttt{A}}^{\text{Cat}} \\[4pt] \tilde D_{{\texttt{A}},\tiny\text{ref}}^{\text{An}} \, (1-y_{\texttt{A}}) \Gamma_{\texttt{A}}^{\text{An}}(y_{\texttt{A}}) \cdot \nu^2 & (\xi,r) \in \Omega_{\texttt{A}}^{\text{An}} \\[4pt] \end{cases}\\[-7pt]\nonumber\end{align}

(2.132) \begin{align} \psi_{\texttt{E}} &\,:\!=\, \begin{cases} \psi_{\texttt{E}}^{\text{Cat}} & \xi \in \Omega^{\text{Cat}} \\[4pt] \psi_{\texttt{E}}^{\text{Sep}} & \xi \in \Omega^{\text{Sep}} \\[4pt] \psi_{\texttt{E}}^{\text{An}} & \xi \in \Omega^{\text{An}} \\[4pt] \end{cases} &&& \hat D_{\texttt{E}} &\,:\!=\, \begin{cases} \psi_{\texttt{E}}^{\text{Cat}} {\boldsymbol{\pi}}_{\texttt{E}}^{\text{Cat}} \tilde D_{\texttt{E}} \cdot \tilde n_{{\texttt{E}},\tiny\text{tot}}(y_{\texttt{E}}) \, \Gamma_{\texttt{E}}(y_{\texttt{E}}) & \xi \in \Omega^{\text{Cat}} \\[4pt] \psi_{\texttt{E}}^{\text{Sep}} {\boldsymbol{\pi}}_{\texttt{E}}^{\text{Sep}} \tilde D_{\texttt{E}} \cdot \tilde n_{{\texttt{E}},\tiny\text{tot}}(y_{\texttt{E}}) \, \Gamma_{\texttt{E}}(y_{\texttt{E}}) & \xi \in \Omega^{\text{Sep}} \\[4pt] \psi_{\texttt{E}}^{\text{An}} {\boldsymbol{\pi}}_{\texttt{E}}^{\text{An}} \tilde D_{\texttt{E}} \cdot \tilde n_{{\texttt{E}},\tiny\text{tot}}(y_{\texttt{E}}) \, \Gamma_{\texttt{E}}(y_{\texttt{E}}) & \xi \in \Omega^{\text{An}} \\[4pt] \end{cases} \end{align}

(2.133) \begin{align} \psi_{\texttt{S}} &\,:\!=\, \begin{cases} \psi_{\texttt{S}}^{\text{Cat}} & \xi \in \Omega^{\text{Cat}} \\[4pt] \psi_{\texttt{S}}^{\text{An}} & \xi \in \Omega^{\text{An}} \\[4pt] \end{cases} &&& \hat S_{\texttt{E}} &\,:\!=\, \begin{cases} \psi_{\texttt{E}}^{\text{Cat}} {\boldsymbol{\pi}}_{\texttt{E}}^{\text{Cat}} \cdot \tilde S_{\texttt{E}} \cdot \tilde n_{{\texttt{E}},\tiny\text{tot}}(y_{\texttt{E}}) \, \Gamma_{\texttt{E}}(y_{\texttt{E}}) & \xi \in \Omega^{\text{Cat}} \\[4pt] \psi_{\texttt{E}}^{\text{Sep}} {\boldsymbol{\pi}}_{\texttt{E}}^{\text{Sep}} \cdot \tilde S_{\texttt{E}} \cdot \tilde n_{{\texttt{E}},\tiny\text{tot}}(y_{\texttt{E}}) \, \Gamma_{\texttt{E}}(y_{\texttt{E}}) & \xi \in \Omega^{\text{Sep}} \\[4pt] \psi_{\texttt{E}}^{\text{An}} {\boldsymbol{\pi}}_{\texttt{E}}^{\text{An}} \cdot \tilde S_{\texttt{E}} \cdot \tilde n_{{\texttt{E}},\tiny\text{tot}}(y_{\texttt{E}}) \, \Gamma_{\texttt{E}}(y_{\texttt{E}}) & \xi \in \Omega^{\text{An}} \\[4pt] \end{cases}\\[-7pt]\nonumber \end{align}

(2.134) \begin{align} \eta_n^{\texttt{E}} &\,:\!=\, \begin{cases} \eta_n^{{\texttt{E}},{\text{Cat}}} & \xi \in \Omega^{\text{Cat}} \\[4pt] 0 & \xi \in \Omega^{\text{Sep}} \\[4pt] \eta_n^{{\texttt{E}},{\text{An}}} & \xi \in \Omega^{\text{An}} \\[4pt] \end{cases} &&& \hat \sigma_{\texttt{E}} &\,:\!=\, \begin{cases} \psi_{\texttt{E}}^{\text{Cat}} {\boldsymbol{\pi}}_{\texttt{E}}^{\text{Cat}} \tilde \Lambda_{\texttt{E}} \cdot \tilde n_{{\texttt{E}},\tiny\text{tot}}(y_{\texttt{E}}) \, y_{\texttt{E}} & \xi \in \Omega^{\text{Cat}} \\[4pt] \psi_{\texttt{E}}^{\text{Sep}} {\boldsymbol{\pi}}_{\texttt{E}}^{\text{Sep}} \tilde \Lambda_{\texttt{E}} \cdot \tilde n_{{\texttt{E}},\tiny\text{tot}}(y_{\texttt{E}}) \, y_{\texttt{E}} & \xi \in \Omega^{\text{Sep}} \\[4pt] \psi_{\texttt{E}}^{\text{An}} {\boldsymbol{\pi}}_{\texttt{E}}^{\text{An}} \tilde \Lambda_{\texttt{E}} \cdot \tilde n_{{\texttt{E}},\tiny\text{tot}}(y_{\texttt{E}}) \, y_{\texttt{E}} & \xi \in \Omega^{\text{An}} \\[4pt] \end{cases}\\[-7pt]\nonumber\end{align}

(2.135) \begin{align} \theta &\,:\!=\, \begin{cases} \theta_{{{\texttt{A}}{\texttt{E}}}}^{\text{Cat}} & \xi \in \Omega^{\text{Cat}} \\[4pt] 0 & \xi \in \Omega^{\text{Sep}} \\[4pt] \theta_{{{\texttt{A}}{\texttt{E}}}}^{\text{An}} & \xi \in \Omega^{\text{An}} \\[4pt] \end{cases} &&& R &\,:\!=\, \begin{cases} \underset{s}{\tilde L}^{\text{Cat}} g\!\left(\underset{s}{\tilde \lambda} {}^{\text{Cat}}\right) & \xi \in \Omega^{\text{Cat}} \\[4pt] 0 & \xi \in \Omega^{\text{Sep}} \\[4pt] \tilde L^{\text{An}} g\!\left({{\underset{s}{\tilde{\lambda}}}\,\!} {}^{\text{An}}\right) & \xi \in \Omega^{\text{An}} \\[4pt]\end{cases}\\[-7pt]\nonumber\end{align}

(2.136) \begin{align} \tilde r_{\texttt{A}} &\,:\!=\, \begin{cases} \tilde r_{\texttt{A}}^{\text{Cat}} & \xi \in \Omega^{\text{Cat}} \\[4pt] \tilde r_{\texttt{A}}^{\text{An}} & \xi \in \Omega^{\text{An}} \\[4pt] \end{cases}\end{align}

which yields

(2.137) \begin{align} \nu^2 \tilde r_{\texttt{A}}^2 C_h \frac{\partial y_{\texttt{A}}}{\partial \tau} & = - \partial_\nu \tilde J_{{{\texttt{A}}_C}} & \text{with} \quad & \tilde J_{{{\texttt{A}}_C}} \!=\! - \hat D_{\texttt{A}} \, \partial_\nu y_{\texttt{A}} && (\nu,\xi) \in \Omega_{\texttt{A}} \end{align}

(2.138) \begin{align} 0 &= -\partial_\xi \tilde J_{{\texttt{S}}_q} - \theta \cdot R & \text{with} \quad & \tilde J_{{\texttt{S}}_q} \!=\! - \hat \sigma_{\texttt{S}} \partial_\xi \tilde \varphi_{\texttt{S}} && \xi \in \Omega^{\text{El}} , \end{align}

(2.139) \begin{align} \psi_{\texttt{E}} C_h c_{\texttt{E}} \frac{\partial y_{\texttt{E}}}{\partial \tau} &= -\partial_\xi \tilde J_{{\texttt{E}}_C} + \eta_n^{\texttt{E}} \!\left(1\!-\!t_{{\texttt{E}}_C}\right) \!\theta \cdot R & \text{with} \quad & \tilde J_{{\texttt{E}}_C} \!=\! - \hat D_{\texttt{E}}(y_{\texttt{E}}) \partial_\xi y_{\texttt{E}} && \xi \in [0,1] , \end{align}

(2.140) \begin{align} 0 &= -\partial_\xi \tilde J_{{\texttt{E}}_q} + \eta_n^{\texttt{E}} \theta \cdot R & \text{with} \quad & \tilde J_{{\texttt{E}}_q} \!=\! - \hat S_{\texttt{E}}(y_{\texttt{E}}) \partial_\xi y_{\texttt{E}} - \hat \sigma_{\texttt{E}} \partial_\xi \tilde \varphi_{\texttt{E}} && \xi \in [0,1] .\end{align}

The boundary conditions read

(2.141) \begin{align} \tilde J_{{{\texttt{A}}_C}} \Big|_{\nu=1} &= - \tilde r_{\texttt{A}} R , & \tilde J_{{{\texttt{A}}_C}} \Big|_{\nu=0} &= 0 , \end{align}

(2.142) \begin{align} \hat \sigma_{\texttt{C}}^{\text{Cat}} \partial_\xi \tilde \varphi_{\texttt{S}}^{\text{Cat}} \Big|_{\xi = 1} &= C_h \eta_W^{\text{Cat}} , & \tilde \varphi_S^{\text{An}} \Big|_{\xi = 0} & = 0 ,\end{align}

with additional homogenous Neumann boundary conditions for all unspecified boundaries.

A remark on the sign of the current density

For the reaction Li⁺ + e^– ⇌ Li + κS, we have $\lambda \,:\!=\, \mu_{Li}|^{\texttt{A}} + \kappa \mu_S|^{\texttt{E}} - \mu_{Li^+}|^{\texttt{E}} - \mu_{e^-}|^{\texttt{A}}$ whereby $\lambda > 0$ entails $ r = L \cdot g(\lambda) > 0$ . Since $J_{\texttt{A}} = - \hat D_{\texttt{A}} \nabla y_{\texttt{A}}$ we have at the boundary

(2.143) \begin{align} J_{\texttt{A}} \cdot {\textbf{n}} = ({+}1) r,\end{align}

where (+1) is the stoichiometric coefficient of the product Li.

2.3. Initial values and potential

The initial values, also used for the Newton solver, are defined as

(2.144) \begin{align} y_{\texttt{A}}(\xi,r,\tau=0) &= y_{\texttt{A}}^0 & y_{\texttt{E}}(\xi,\tau=0) &=y_{\texttt{E}}^0 \end{align}

(2.145) \begin{align} \tilde\varphi_{\texttt{E}}(\xi,\tau=0) &=\tilde \varphi_{\texttt{E}}^0 & \tilde\varphi_{\texttt{S}}(\xi,\tau=0) &=\tilde \varphi_{\texttt{S}}^0 ,\end{align}

with

(2.146) \begin{align} y_{\texttt{A}}^0 &\,:\!=\, \begin{cases} y_{\texttt{A}}^{{\text{Cat}},0} & \qquad (\xi,r) \in \Omega_{\texttt{A}}^{\text{Cat}} \\[4pt] y_{\texttt{A}}^{{\text{An}},0} &\qquad (\xi,r) \in \Omega_{\texttt{A}}^{\text{An}} \\[4pt] \end{cases} ,&&& y_{\texttt{E}}^0 &\,:\!=\, \frac{n_{{\texttt{E}}_C}}{n_{{\texttt{E}}_S} - 2 \cdot \kappa_{\texttt{E}} \cdot n_{{\texttt{E}}_C}} ,\\[-4pt]\nonumber \end{align}

(2.147) \begin{align} \tilde \varphi_{\texttt{S}}^0 &\,:\!=\, \begin{cases} \left. f_{\texttt{A}}^{\text{An}}\left(y_{\texttt{A}}^{{\text{An}},0}\right) - f_{\texttt{A}}^{\text{Cat}}\left(y_{\texttt{A}}^{{\text{Cat}},0}\right)\right) & \qquad \xi \in \Omega^{\text{Cat}} \\[4pt] 0 & \qquad \xi \in \Omega^{\text{An}} \\[4pt] \end{cases} , &&& \tilde \varphi_{\texttt{E}}^0 &\,:\!=\,f_{\texttt{A}}^{\text{An}}\left(y_{\texttt{A}}^{{\text{An}},0}\right) - f_{\texttt{E}}\left(y_{\texttt{E}}^0\right).\end{align}

2.3.1. Parameters

All parameters and their values for the subsequent numerical calculations are summarised in Appendix B.

3. Discretisation and model order reduction of the battery model

The modelling approach discussed in the previous section formulates the multi-scale battery model as a homogeneous medium in which electrolyte and the solid electrode materials exist at every point. This homogenisation approach results in a so-called macroscopic model where the details of the micro-structure are taken into account. In this macroscopic model, the intercalation of Li-ions in the electrode particles is incorporated through a coupled diffusion equation in radial direction of the particles in each macroscopic quadrature point. In this way, we get a pseudo-2D model (i.e. 1D + 1D model) for the full battery cell. The model is given by a system of nonlinear PDEs for the homogenised electric potentials $(\tilde{\varphi}_E , \tilde{\varphi}_S )$ and the mole fraction of $\text{Li}^+$ -ions $(y_E, y_A)$ in the electrolyte and in the positive and negative electrode materials, respectively. The equation for the mole fraction in the electrode materials is calculated in the additional pseudo dimension, namely in the particle radius $\nu$ . The overall PDE system can be written in the following abstract form:

\begin{align*}- \alpha(\cdot, y_A) \, \frac{\partial y_A}{\partial \tau}& - \nabla_{\nu} \cdot [{-} \beta(\cdot,y_A) \nabla_{\nu} y_A] & =& \,\, 0, \\&- \nabla_{\xi} \cdot [{-} \gamma(\cdot) \nabla_{\xi} \tilde{\varphi}_S] &- R_1\!\left(\cdot ,\tilde{\varphi}_E, \tilde{\varphi}_S, y_E, y_A \vert_ {\nu=1}\right) =& \,\, 0 , \\-\delta(\cdot, y_E) \, \frac{\partial y_E}{\partial \tau} &- \nabla_{\xi} \cdot [{-} \kappa(\cdot,y_E) \, \nabla_{\xi} y_E ] & + R_2\!\left(\cdot ,\tilde{\varphi}_E, \tilde{\varphi}_S, y_E, y_A \vert_ {\nu=1}\right) =& \,\, 0, \\&- \nabla_{\xi} \cdot [{-} \omega (\cdot, y_E) \, \nabla_{\xi} y_E - \rho(\cdot, y_E) \, \nabla_{\xi} \tilde{\varphi}_E ]\!\!\!\!\! &\!\!\!\!\!+ R_3\!\left(\cdot ,\tilde{\varphi}_E, \tilde{\varphi}_S, y_E, y_A \vert_ {\nu=1}\right) =& \,\, 0.\end{align*}

The linear and nonlinear coefficient functions $\alpha, \beta, \gamma, \delta, \kappa, \omega$ and $ \rho$ correspond to the representation from the equations (2.137)–(2.140) and depend on the domain for which the system is defined (anode, cathode and separator). $R_i , i = \{1, 2, 3, 4\}$ represents the reaction rate functions (2.135) in addition to the previous constants. The system is completed by the boundary conditions as well as interface conditions (2.141)–(2.142). Note that there are corresponding Neumann boundary conditions $\beta(\cdot,y_A) \nabla_{\nu} y_A \cdot n = - R_4(\cdot,\tilde{\varphi}_E, \tilde{\varphi}_S, y_E, y_A)$ at the boundary of the electrode particles, i.e. at $\nu = 1$ , which couples the microscopic equation with the macroscopic equations.

In Section 3.1, we briefly introduce the discretisation of the above battery PDE system by the finite element method in space and the backward Euler method in time. Numerical simulations such as high-fidelity approximations techniques like the finite element method can become prohibitively expensive when we expect them to deal quickly and efficient with the repetitive solution of PDEs. In Section 3.2,we present the reduced basis method that replaces the high-fidelity problem by a reduced order model (ROM) of much lower numerical complexity. To achieve a offline–online decomposition, the discrete empirical interpolation method (DEIM) is performed on the full finite element model operator in Section 3.3.

3.1. Discretisation

For the battery model in the abstract form above, let $\Omega_{1D}$ denote a computational one-dimensional domain in the macroscopic direction and $\Omega_{\nu} = \Omega^\xi_{\nu}$ the transverse/radial directions in the electrode particles associated with each $\xi \in \Omega_{1D}$ . Let $\mathscr{P} \subset \mathbb{R}^P, P \geqslant 1$ denote the parameter space. We define the solution space $V = V_1 \oplus V_2$ with $ V_1 = H^1(\Omega_{\nu}), V_2 = (H^1(\Omega_{1D}))^3$ and $V^{\prime}$ , the dual space of V. For a corresponding variational weak formulation, we obtain, after a semi-discretisation in time t by the implicit Euler method, that the battery model can be formulated as the following nonlinear system:

(3.1) \begin{align}\text{Find}\ u^{t+1}=\left[u^{t+1}_1, u^{t+1}_2, u^{t+1}_3, u^{t+1}_4\right]^\top \in V: \quad G_{\mu}(u^{t+1},\psi)=f_{\mu}(u^{t},\psi) \ \ \forall \ \psi \in V,\end{align}

where the operator $G_{\mu}(\cdot,\cdot)\,:\, V \times V \rightarrow \mathbb{R}$ represents the non-linear time-discrete PDE system. The index $\mu \in \mathscr{P}$ indicates the dependence of the problem on certain parameters, such as the charge rate, the diffusion coefficient or the reaction rate. $f_{\mu}(u^{t}, \cdot) \in V^{\prime}$ contains the solid potential Neumann boundary conditions.

For the discretisation in space, the finite element method is used [79], i.e. we project (3.1) to a finite dimensional, continuous and piecewise polynomial space $V_h \subset V$ . We hence obtain for each time step a fully discrete non-linear system of the form:

(3.2) \begin{align}\text{Find}\ u_h^{t+1}=\left[u^{t+1}_{1_h}, u^{t+1}_{2_h}, u^{t+1}_{3_h}, u^{t+1}_{4_h}\right]^\top \in V_h\,:\, \quad G_{\mu}\!\left(u_h^{t+1},\psi_i\right)=f_{\mu}\!\left(u_h^{t},\psi_i\right) \quad \forall \, i=1,\dots n,\end{align}

where $\psi_i, i=1, \dots, n$ denotes the standard Lagrange basis of the finite element space $V_h = \bigoplus_{i=1}^4 V_{i_h}$ . Henceforth, the operator $G_{\mu}$ can be called the finite element operator which operates on $V_h \times V_h$ . The developed discretisation does not depend on the specific choice of the parameters to be varied.

3.2. Reduced basis method

The computational studies of the battery model require many evaluations of the finite element system (3.2) with different parameter settings and therefore involve a large amount of time and experimental effort. Hence, we apply the reduced basis method to further reduce the parameterised discretised battery model to obtain a reduced order model that is cheap to evaluate, as e.g. detailed in [73]. Due to the parametrisation of the battery model, the finite element solution space is restricted to a solution manifold that contains the solutions relevant for the considered application. The reduced basis method is based on the idea of performing a Galerkin projection of the discrete equations onto low-dimensional subspaces $\tilde{V} \subset V_h$ that approximates this solution manifold in order to accelerate the repeated solution of (3.2) for varying parameters $\mu$ . Under this projection, the reduced problem is given by

(3.3) \begin{align}\text{Find}\ \tilde{u}^{t+1}=\left[\tilde{u}^{t+1}_1, \tilde{u}^{t+1}_2, \tilde{u}^{t+1}_3, \tilde{u}^{t+1}_4\right]^\top \in \tilde{V}\,:\, \quad G_{\mu}\!\left(\tilde{u}^{t+1},\tilde{\psi}_i\right)=f_{\mu}(\tilde{u}^{t},\tilde{\psi}_i)\ \forall \, i=1,\dots m,\end{align}

where $m \ll n$ and $\tilde{\psi}_i, i=1,\dots,m$ , represents the basis of the reduced (basis) space $\tilde{V}$ . The basic idea of the reduced basis method is to perform a so-called offline/online decomposition. In the preceding offline phase, high-dimensional computations are performed for a chosen set of parameters to construct the reduced space such that a small error between the high-dimensional and the reduced approximation can be certified for each valid parameter value. After projecting onto this reduced space, the resulting low-dimensional problem (ROM) can be solved for any suitable parameter value in a following online phase. Various methods for constructing the reduced space for time-dependent problems have been considered in the literature, such as the POD-Greedy [37]. The POD-Greedy method produces approximation spaces with quasi-optimal $l^{\infty}$ -in- $\mu$ , $l^2$ -in-time reduction error [36].

As the POD-Greedy method relies on the usage of rigorous and efficient to evaluate a posteriori error estimators, which are not available for the non-linear battery model at hand, in our numerical experiments, the reduced space is generated by the POD method. The basic idea of POD method is that a pre-selected set of solutions trajectories of the problem (3.2), called snapshots, are transformed into a new set of uncorrelated variables, called POD modes. The first few modes contain the most relevant features from the battery system and form the reduced basis. In our case, we apply the POD method separately to the respective solution components, cf. [69]. More precisely, let $\mathcal{P} = \{ \mu^1, \dots, \mu^{n_s} \}$ be a set of $n_s$ parameter samples and $\{ u_h(\mu^1), \dots, u_h(\mu^{n_s}) \}$ the corresponding snapshot set. At each time step, the snapshot of the set is calculated using the following prescription of Newton’s method:

\begin{align*}D_uF_{\mu}\!\left(u_h^{t+1,k},\psi_i\right) \ \delta u_h &= - F_{\mu}\!\left(u_h^{t+1,k},\psi_i\right) \qquad \forall \, i=1,\dots,n,\\[4pt]u_h^{t+1,k+1}(\mu) &= u_h^{t+1,k}(\mu) + \delta u_h,\end{align*}

where $D_u F_{\mu}(z,\psi_i)$ is the Fréchet derivative of $F_{\mu} (u,\cdot) = G_{\mu}(u,\cdot) - f_{\mu}(u_h^t,\cdot)$ with respect to u at $z \in V_h$ . We separate the snapshots into the respective components, $u_{1_h}, \ldots, u_{4_h}$ and generate the reduced basis separately. We define the corresponding snapshot matrices $S_i \in \mathbb{R}^{N_h^{\iota} \cdot \vert t \vert \times n_s}$ with $i=1,2,3,4$ and $\iota \in \{ \Omega_{1D}, \Omega_{\nu}\}$ as

\begin{align*}S_i &= \left[ u_i^1, \dots, u_i^{n_s}\right],\end{align*}

where the vectors $u_i^j \in \mathbb{R}^{N_h^{\iota} \cdot \vert t \vert},\, 1 \leqslant j \leqslant n_s,\,$ denote the degrees of freedom of the functions $u_{i_h}(\mu^j)\vert_{t} \in V_{i_h}$ . The singular value decomposition of $S_i$ is given through $S_i = U_i \Sigma_i Z_i^T,$ where $U_i \in \mathbb{R}^{N_h^{\iota} \cdot \vert t \vert \times N_h^{\iota} \cdot \vert t \vert}$ and $Z_i \in \mathbb{R}^{ns \times ns}$ represent orthogonal matrices, and $\Sigma_i = \text{diag}\!\left(\sigma_i^1, \dots, \sigma_i^{z_i}\right) \in \mathbb{R}^{N_h^{\iota} \cdot \vert t \vert \times ns}$ with $\sigma_i^1 \geqslant \sigma_i^2 \geqslant \dots \geqslant \sigma_i^{z_i}, \, z_i \leqslant \min\!\left(N_h^{\iota} \cdot \vert t \vert, ns\right)$ contain the singular values. The left singular vectors

\begin{align*}U_i = \left[ \zeta_i^1 |\dots | \zeta_i^{N_h^{\iota}\cdot \vert t \vert} \right]\end{align*}

span the reduced space $\tilde{V_{i}}$ using only the singular vectors whose singular values are above a fixed threshold value $\varepsilon$ . Due to the fundamental properties of POD, the projection error consist of the $l^2$ -sum of the corresponding truncated singular values, and $\tilde{V}_i, i=1,2,3,4,$ are $l^2$ -in-space, $l^2$ -in-time best-approximation spaces for the considered training set of snapshots. The overall reduced space is defined as $\tilde{V} \,:\!=\, \bigoplus\limits_{i=1}^4 \tilde{V_{i}}$ .

3.3. Empirical interpolation and hierarchical approximate POD

The model reduction approach described in Section 3.2 still suffers from a huge offline cost, due to the global POD construction of the reduced basis and from an inefficient online computational cost, as the reduced model (3.3) cannot be decomposed efficiently into high- and low-dimensional computations due to the presence of non-linearities in the system. In order to obtain a more efficient ROM with reduced computational cost in the offline phase and a full so called offline–online decomposition, we replace the global POD reduced basis construction by an incremental hierarchical approximate POD (HAPOD) [42] and construct an affine approximation of the non-linear differential operator $G_{\mu}$ by an empirical operator interpolation [13, 28].

In detail, the offline–online decomposition consists of precomputing parameter- and solution-independent terms, a so called collateral basis that allows to interpolate evaluations of the operator $G_{\mu}$ . The construction of the collateral basis and associated interpolation functionals is done once in the offline phase to allow a fast evaluation of the interpolated operator $I_M(G_{\mu})$ during the online phase. In detail, the empirical operator interpolation generates a separable approximation by interpolating at M selected degrees of freedom of $V_h$ . The approximation of the operator is of the following form:

(3.4) \begin{align}G_{\mu}(\tilde{u},\cdot) \approx I_M(G_{\mu}(\tilde{u},\cdot)) = \sum_{q=1}^M \theta_{\mu}^q(\tilde{u}) G^q,\end{align}

where $\{ G^q\}_{q=1}^M$ is the collateral basis, i.e. a basis for a subspace of

\begin{equation*}\mathscr{M}_G = \left\{ G_{\mu}(u_h,\cdot) | \mu \in \mathscr{P} \right\},\end{equation*}

and $\theta_{\mu}^q$ are interpolation coefficients recalculated for each $\mu$ and $\tilde{u}$ during the online phase. The collateral basis is obtained by applying the POD method to the set of snapshots obtained as images under the operator, i.e. $G_{\mu}(u_h,\cdot)$ . The set of snapshots thereby includes the Newton stages in addition to the corresponding solution trajectory of $G_{\mu}(u_h,\cdot)$ . Moreover, in analogy to the construction of the reduced basis described above, also the collateral basis is constructed separately for each solution component. Hence, we define for $i=1,2,3,4$ :

\begin{align*}\left[ G_i^1, \dots, G_i^{M_i} \right] = \text{POD}\!\left( \left[G_{\mu_1}^1(u_i), \dots, G_{\mu_1}^{s_1}(u_i), G_{\mu_2}^1(u_i), \dots, G_{\mu_{n_s}}^{s_{n_s}}(u_i)\right], \epsilon_{POD} \right),\end{align*}

where the vectors $G_{\mu_j}^k(u_i) \in \mathbb{R}^{N_h^\iota}, 1 \leqslant j \leqslant n_s,\, 1 \leqslant k \leqslant s_j,$ represent the degrees of freedom of the functions $G_{\mu_j}\!\left(u_{i_h}^k(\mu^j)\right)$ and $\epsilon_{POD}$ the error tolerance for the POD. We define $\left[ G^1, \dots, G^M \right] = \left[ G_1^1, \dots, G_4^{M_4} \right]$ with $M = \sum_{i=1}^4 M_i$ .

Figure 4. Sketch of the tree structure of the incremental HAPOD, introduced in [42].

To calculate the interpolation coefficients $\theta_{\mu}^q(\tilde{u}), q=1, \dots,M$ for given $\mu$ and $\tilde{u}$ , the interpolation constraints are imposed at M interpolation points. The interpolation points are selected iteratively from the indices of basis $\{ G^q\}_{q=1}^M$ using a greedy procedure. This procedure determines each new interpolation point by the minimisation of the interpolation error over the snapshots set measured in the maximum norm. For more details, we refer to [2, 13].

By replacing the operator $G_{\mu}$ in (3.3) by the fast-to-evaluate interpolated operator $I_M (G_{\mu})$ , we obtain the completely offline–online decomposable reduced problem

(3.5) \begin{align}\text{Find}\ \tilde{u}^{t+1} \in \tilde{V}\,: \quad I_M \!\left(G_{\mu}\!\left(\tilde{u}^{t+1},\tilde{\psi}_i\right)\right)=f_{\mu}(\tilde{u}^t,\tilde{\psi}_i)\ \forall \, i=1,\dots m,\end{align}

which can be solved efficiently for varying parameters $\mu$ .

For large-scale time-dependent applications such as our battery model, computing the POD algorithm can be expensive. Especially if we include the evaluations of $G_{\mu}$ at all Newton levels of the selected solution trajectories in the operator snapshot set. The hierarchical approximate POD (HAPOD) algorithm is an efficient approach, which approximates the standard POD algorithm based on tree hierarchies, where the task of computing a POD for a given large snapshot set S is replaced by multiple small PODs [42]. More specifically, we use the special case of incremental HAPOD. In this case, the tree structure is such that each node of this tree represents either a leaf or has exactly one leaf and one non-leaf as children, cf. Figure 4. In detail, first the vectors of the given snapshot set are assigned to the leaves of the tree $\beta_1, \cdots , \beta_B$ . Starting with two leaves $\beta_1, \beta_2$ , a POD of each local snapshot data is computed. The resulting modes scaled by their singular values are the input to the parent node $\alpha_1$ , which is again used to calculate a POD. This newly generated input and the local snapshot data assigned to leaf $\beta_3$ are the input of the parent node $\alpha_2$ . The final HAPOD modes $\rho$ are reached, when the last leaf $\beta_B$ has entered. For the calculation of the collateral basis, we e.g. define for $i = 1,2,3,4$ :

\begin{align*}\left[ G_i^1, \dots, G_i^{M_i} \right] = \text{HAPOD}\!\left( \left[G_{\mu_1}^1(u_i), \dots, G_{\mu_1}^{s_1}(u_i), G_{\mu_2}^1(u_i), \dots, G_{\mu_{n_s}}^{s_{n_s}}(u_i)\right], \epsilon_{POD}, \omega \right),\end{align*}

where $\epsilon_{POD}$ is the desired approximation error tolerance for the resulting HAPOD space. Depending on $\omega$ , one might get more modes than needed for a POD with the same tolerance. The local tolerances in the HAPOD algorithm are computed from $\epsilon_{POD}$ and $\omega \in (0,1)$ . More details can be found in [42].

4. Numerical results

In order to create a test environment for our modelling framework, we develop an experimental implementation of the ageing effects of the battery model from Section 2. We investigate the efficiency of the reduced order simulations by evaluating electrochemical characteristics over the cyclisation $n = 1, \dots, 1000$ for different ageing models. The electrochemical characteristics are the voltage-capacity spectrum $E(\bar y_A , n;\, C_h)$ and the status of charge $\bar y_A^{Cat}$ at a specific voltage value $E_{\min}$ (see equations (2.102)–(2.109)). We assume that the ageing effects are modeled by given functions dependent on the cycle number n for the reaction rate $\hat L(n)$ and diffusion coefficient $\hat D_A(n)$ . These functions are used to investigate the qualitative behaviour of the ageing effects. In this context, we assume that the parameter dependence can be interpreted as follows; The parameter dependence of the reaction rate $\hat L$ should investigate the degradation of the solid electrolyte interphase. This might illustrates the increase in reaction resistance due to cyclisation. Furthermore, we want to investigate the degradation of the porous electrodes, which can be interpreted by a decrease in the diffusion coefficient $\hat D_A$ . This effect can be caused by micro-cracks within a particle. To efficiently analyze this forward modelling of ageing effects, we consider three scenarios. In the first scenario, we consider the unaged battery and calculate the voltage against the state of charge for varying charge rates $C_h$ . In the next case, we set $C_h = 1$ and examine the ageing effects of $\hat D_A$ and $\hat L$ by alternately choosing one of the parameters fixed. In the last case, we vary all parameters $\hat D_A, \hat L$ and $C_h$ .

It should be noted that the degradation experiments are intended to serve as a proof of concept. The next step would be to compare the battery model with real data. However, that is not the focus. The focus of this work is to present the battery modelling framework and provide a proof of concept through forward modelling of the ageing effects. Therefore, in a very general setting, our work provides the methodology to systematically investigate degradation mechanisms on the basis of cycle number dependent parameters.

4.1. Implementation aspects

Let us first introduce the settings used in the numerical experiments. We implement a (pseudo) 2D grid, where the bottom of the grid corresponds to the 1D grid on which the macroscopic equations are computed (see Figure 5). The length of the bottom $L_{cat} + L_{sep} + L_{ano}$ is divided into $N_h^{\Omega_{1D}}=300$ grid points. Here, each of the cathode, separator and anode is discretised into 100 grid points. The pseudo-dimension for the microscopic equation, which goes vertically at each bottom grid point, consists of $N_h^{\Omega_{\nu}}$ grid points. In the simulations $N_h^{\Omega_{\nu}} = 100$ is chosen. All in all, this results in 30,900 degrees of freedom (dofs) for the overall system. The Neumann boundary condition of the microscopic equation couples the equation to the macroscopic equations and is defined for $\nu = 1$ . To ensure that the bottom grid corresponds to $\nu=1$ a transformation of the microscopic equation with $\tilde{\nu}=(1-\nu)$ is performed. In addition, as an essential step for the stability of the model, a variable transformation of the microscopic variable $y_A$ of the following form is applied:

\begin{align*}y_A(g) = \frac{e^g}{1+ e^g}, \quad g(y_A) = \text{ln} \!\left( \frac{y_A}{1-y_A} \right).\end{align*}

The time discretisation is performed by an implicit Euler method on a $T = 1$ time interval with a time step size of $\Delta t = 10^{-2}$ . The nonlinear battery system is solved by a Newton method to a relative error accuracy of $10^{-5}$ and a termination condition of $\min_{\xi \in \Omega^{Cat}} \tilde{\varphi}_S(\xi) \leqslant E_{\min}$ with the voltage value $E_{\min}=-0.2 $ . To generate the reduced space $\tilde{V}$ , we compute a snapshot set $S_{train}$ on training sets of equidistant parameters. For the experiments 1 and 3, we choose $\vert \mathcal{P}_{train} \vert = 15$ and for experiment 2 we choose $\vert\mathcal{P}_{train} \vert= 10$ .

Figure 5. Sketch of the pseudo 2D grid. The blue line shows the computational domain of the macroscopic equations, while the red lines illustrate the computational domain for the microscopic equation.

The analytical solution of a partial differential equation is called an exact solution. In many applications, an analytical solution is not available, especially for the battery problem at hand. The finite element method is an approximation approach. We call the high-dimensional finite element solution the truth (full order) solution. The reduced basis method reduces the number of unknowns compared to the finite element model by Galerkin projection. The solution of this reduced model is called reduced solution. Note that the focus is on approximating the truth solution, assuming that the discretisation error between the truth and the respective exact solution is negligible. In this work, we do not consider the discretisation error. Our focus is on the error between the reduced and the truth solution and the corresponding speed up due to the model order reduction. As a measure for the model reduction error, we determine the relative $l^2$ -in-space, $l^2$ - in-time error averaged over a set of random test parameters $\mathcal{P}_{test}$ given by

(4.1) \begin{align} \frac{1}{\vert \mathcal{P}_{test} \vert}\sum_{\mu \in \mathcal{P}_{test}} \frac{\vert\vert u_h(\mu) - \tilde{u}(\mu) \vert\vert_2}{\vert \vert \tilde{u}(\mu) \vert \vert_2}.\end{align}

For the truth solution, we want to ensure that we have achieved grid convergence up to a certain error $\varepsilon$ . Thus, we assume that the high-dimensional discretisation adopts a resolution of $\vert\vert u_n - u_m \vert\vert_2 < \varepsilon \leqslant 10^{-5}$ with 30,900 grid points for $u_n$ and 31,512 grid points $\left(N_h^{\Omega_{1D}}=303, N_h^{\Omega_{\nu}} = 101\right)$ for $u_m$ .

The empirical operator interpolation approximates the full discrete and nonlinear model operator. If we do not interpolate the nonlinearities of the model operator sufficiently well, then the system is unstable. We detect the instabilities by the model reduction error and then increase the dimension of the empirical interpolation accordingly. More precisely, we determine the basis using the HAPOD method. As a result, we get a set of basis vectors. From this set we take the first basis vectors, check the corresponding system for instabilities and increase the number of basis vectors until we achieve stability and the desired model reduction error. In this way we determined the dimension of the empirical interpolation in the experiments. For a principal algorithmic approach, we refer to [5, 28]. We proceed in the same way for the reduced basis.

The programming language is Python. All simulations of the high-dimensional model are computed with the finite element sofware NGSolve [67]. NGSolve provides the ability to construct the complex grid structure and define the battery model for each subdomain. For the implementation of the reduced basis method, the NGSolve code has integrated into the model order reduction library pyMOR [72]. We include the evaluations of $G_{\mu}$ on all Newton stages of the selected solution trajectories in the operator snapshot set to compute the collateral basis. This leads to a stabilization of the reduced model. In order to speed up the computation of the collateral basis for the empirical interpolation data via POD, the HAPOD algorithm is used instead of the standard POD algorithm. For the generation of the reduced basis, we use the HAPOD algorithm as well. We choose $\epsilon = 4e-8$ and $\omega = 0.9$ in both cases. For illustration purposes, the reduced space and empirical interpolation are calculated for a training set $\vert\mathcal{P}_{train} \vert= 5$ for the variation of $C_h$ . In this case, the calculation without using HAPOD takes $107.32$ min, while only $7.31$ min CPU-time were needed using HAPOD. This corresponds to a speedup of $14.68$ in the offline phase. All tests were performed on the same computer and software basis.

4.2. Experiment 1

In this subsection, we consider the variation of the charge rate $C_h$ with $C_h \in [0.01,4]$ for an unaged battery. In this case, we choose $\hat L=0.5$ and $\hat D_A=0.5$ . For the reduced order model the number of basis functions for the four variables are set to $ \# u_1 = 3, \#u_2 = 3, \#u_3 = 5$ and $\#u_4 = 4$ . For the empirical operator interpolation, the number of interpolation points is $ \# G(u_1) = 19, \#G(u_2) = 15, \#G(u_3) = 60$ and $\#G(u_4) = 8$ . These numbers are obtained by calculating the relative model reduction error (4.1) with successive increase of the basis size up to an accuracy of order $10^{-5}$ (see Table 1). To ensure the stability of the reduced model for empirical operator interpolation, the number of interpolation points (or collateral basis) must be chosen large enough. Especially the number of interpolation points for $G(u_3)$ are crucial here, since to achieve stability we need a much higher number of interpolation points than for $G(u_1), G(u_2)$ and $G(u_4)$ .

Figure 6. (a) Voltage-capacity spectrum compared to the open circuit potential. Solution plots of the four components, (b)-(c) with $C_h = 0.1$ and (d)-(e) with $C_h = 4$ for $t=0.2$ .

The voltage-capacity spectrum is shown in Figure 6(a)), where we achieve a model reduction error of less than $10^{-4}$ for a simulation time of 2.58 min. Therefore, we obtain a speed up of 15.41. In addition, during the computation of the voltage-capacitance spectrum, the four components of the solution of the battery model are determined and stored. As an example, two solution plots at a fixed time $t=0.2$ for the charge rate $C_h \in \{ 0.1, 4 \}$ are illustrated in Figure 6(b)–(e). At a low charge rate, we almost reach the OCP, which can be observed by the fact that nearly constant functions are obtained. For larger charge rates, we observe higher gradients in macroscopic and microscopic directions due to transport limitations. Table 1 presents the average simulation times and the corresponding relative model reduction errors for the computation of a single battery solution trajectory for different reduced basis sizes.

Table 1. Relative model reduction error (4.1) and reduced simulation times for a battery simulation trajectory with $\hat L=0.5, \hat D_A=0.5$ and $\mathcal{P}_{test} = 10$ . The number of the reduced basis consists of the four variables, e.g. $11 = \# u_1 + \#u_2 + \#u_3 +\#u_4 = 2+2+4+3$ . When the reduced basis is increased, each variable is added one basis. The number of interpolation point is 102. The average time for the solution trajectory of the high-dimensional model is 363.48 s

4.3. Experiment 2

In the following, the evaluation of the degradation simulations of the solid electrolyte interphase and the porous electrodes are presented. To investigate the qualitative behavior of these ageing effects, we consider the decrease in value of the parameters $\hat L$ and $\hat D_A$ equally in anode and cathode over n with $n = 0 ,\dots, N$ cycles and set $C_h =1$ . A cycle consists of a discharge process, assuming that the battery is charged in such a way that the chosen initial conditions apply at the beginning of each cycle. In addition, we assume that the evolution of the parameters $\hat L(n)$ and $\hat D_A(n)$ as a function of the number of cycles n satisfies an ordinary differential equation

(4.2) \begin{align}\frac{\textrm{d} F(n)}{\textrm{d} n} &= a_F \, F(n),\end{align}

with the unknown parameter $a_F(\beta)$ such that:

(4.3) \begin{align}F(0) &= F_0, \\F(N) &= \beta F_0, \quad \beta < 1. \notag\end{align}

It follows that $F(n)= F_0 \, e^{\log (\beta) n}/{N}$ and $a_F(\beta) = \frac{\log(\beta)}{N}$ (see Figure 7). $F_0$ is the corresponding initial parameter value. In our case, $\hat D_{A,0}= \hat L_0 = 0.5$ . Under this assumption, we impose that the ageing mechanism in one cycle depends on the value of the degradation of the previous cycle. (This assumption is based on the fact that under laboratory conditions, the cell is always discharged in the same way.)

The characteristic spectrum of cell voltage E for the degradation of reaction rate $\hat L$ for $\hat D_A=0.5$ is shown in Figure 8(a), (c) for two different choices of $\beta$ . Furthermore, the capacity $\bar y_A^{Cat}$ at the specific voltage value $E_{\min}=-0.2$ over the number of cycles $N = 1000$ with a variation of $\beta$ is illustrated in Figure 8(e). The same scenario is shown for the degradation of the diffusion coefficient $\hat D_A$ for $\hat L=0.5$ in Figure 8(b),(d),(f). As expected, the graphs show when the parameters degrade faster and more severely, the cell voltage and capacity decrease more rapidly. However, the variation of $\hat L$ has a larger effect in terms of absolute capacity loss.

Figure 7. Various evolutions of the parameter functions satisfying the ordinary differential equation (4.2) with $F_0 =0.5$ and $N = 1000$ .

Figure 8. (a)-(d) Evolution of the capacity-dependent voltage E of $\hat D_A(n)$ and $\hat L(n)$ compared to the open circuit potential for $\beta = 0.1$ or $\beta = 0.4$ . (e)-(f) Effect of the different degradation models of $\hat D_A(n)$ and $\hat L(n)$ on the capacity at voltage $E_{\min}$ over the number of cycles n at $C_h=1$ .

For the implementation of the reduced order model (ROM) with dependence on the parameters $\mu=[ \hat D_A$ , $\hat L ]$ , the number of basis functions for the four variables is set to $ \# u_1 = 3, \#u_2 = 3, \#u_3 = 5$ and $\#u_4 = 3$ . For the empirical operator interpolation, the number of interpolation points is $ \# G(u_1) = 9, \#G(u_2) = 9, \#G(u_3) = 15$ and $\#G(u_4) = 9$ . As in Experiment 1, these numbers are obtained by calculating the relative model reduction error (4.1), taking into account an accuracy of order $10^{-6}$ with successive increase of the basis size. To ensure the stability of the reduced model for empirical operator interpolation, the number of interpolation points must be increased for larger reduced basis space dimensions.

In Table 2, we observe a rapid decay of the model reduction error, which stagnates already for relatively small reduced basis space dimensions. In this manner, we obtain relative reduction errors as small as $10^{-5}$ with battery simulation times of less than 10 s. Note that the different conditions between experiment 1 and experiment 2 that leads to a significant reduction in time is that we have 102 interpolation points in experiment 1 and only 42 interpolation points in experiment 2. Thus, the approximation of the battery operator via empirical operator interpolation in experiment 2 has significantly fewer computational operations, which is significant in terms of time. When calculating the voltage spectra, we achieve an average relative reduction error of about $10^{-3}$ and an average speed up of 24.37. Calculating the capacity at the voltage value $E_{\min}$ over the number of cycles requires about 118.92 h for the full model. By using the reduced model, we obtain approximations in about 2.53 h with a relative error of $10^{-5}$ . It is a speed up of 46.83.

Table 2. Relative model reduction error (4.1) and reduced simulation times for a battery simulation trajectory with $C_h=1$ and $\mathcal{P}_{test} = 10$ . The number of the reduced basis consists of the four variables, e.g. $10 = \# u_1 + \#u_2 + \#u_3 +\#u_4 = 2+2+4+2$ . In each column, a basis is added to each variable. The number ob interpolation points amounts to 42. The average time for the solution trajectory of the high-dimensional model is 356.49 s

4.4. Experiment 3

In the previous section, degradation simulations of the solid electrolyte interphase and the porous electrodes were considered by decreasing the value of the parameters $\hat L$ and $\hat D_A$ over N cycles. We assumed that the value of the parameter depends on the value of the parameter in the previous cycle. Thereby, the charge rate $C_h$ was set constant to 1 for each simulation. In this experiment, a variation of the charge rate is a matter of interest as well. This implies that the ROM depends on the following parameter vector $\mu=[C_h, \hat D_A, \hat L ]$ . Furthermore, as in Experiment 2, we assume that the evolution of the parameters $\hat D_A$ and $ \hat L$ as a function of the number of cycles n satisfies the ordinary differential equation that additionally depends on the charge rate $C_h$

(4.4) \begin{align}\frac{\textrm{d} F(n)}{\textrm{d} n} &= a_F \, F(n) \, C_h,\end{align}

with the unknown parameter $a_F(\beta)$ such that:

(4.5) \begin{align}F(0) &= F_0, \nonumber\\F(N) &= \beta F_0, \quad \text{for} \, \beta < 1 \, \text{at} \ C_h=1. \end{align}

It follows that $F(n)= F_0 \, e^{{ C_h \log\! (\beta) n}/{N}}$ and $a_F(\beta) = \frac{\log(\beta)}{N}$ (see Figure 9). Here we set $\beta = 0.6$ and $F_0$ again represents the corresponding initial parameter value. Note that for $C_h=1 $ the situation is the same as in Experiment 2.

In this experiment, we set the number of the basis functions for the four variables to $ \# u_1 = 4,$ $ \#u_2 = 4, \#u_3 = 5, \#u_4 = 4$ and the number of interpolation points to $ \# G(u_1) = 22, \#G(u_2) = 20, \#G(u_3) = 110, \#G(u_4) = 20$ , compare Table 3, taking into account an accuracy of order $10^{-5}$ . Note, as in the first experiment, that the number of interpolation points, in particular the number of interpolation points for $G(u_3)$ , must be chosen large enough to ensure the stability of the reduced model.

Table 3. Relative model reduction error (4.1) and reduced simulation times for a battery simulation trajectory with $\mathcal{P}_{test} = 10$ . The number of the reduced basis consists of the four variables, e.g. $13 = \# u_1 + \#u_2 + \#u_3 +\#u_4 = 3+3+4+3$ . When the reduced basis is increased, each variable is added one basis. The number ob interpolation points amounts is 172. The average time for the solution trajectory of the high-dimensional model is 357.31 s

Figure 9. Various degradation models that satisfy the ordinary differential equation (4.4) with $F_0 =0.5, \beta=0.6$ and $N = 1000$ .

Figure 10. (a)-(b) The spectrum of cell voltage E for the degradation of $\hat D_A(n)$ and $\hat L(n)$ compared to the open circuit potential for $C_h = 2$ . (c)-(d) The effect of the different degradation models of $\hat D_A(n)$ and $\hat L(n)$ on the capacity at voltage $E_{\min}$ over the number of cycles n with $\beta=0.6$ .

Comparing the run times from the calculation of the spectrum of the cell voltage (see Figure 10(a)–(b)) for the full and reduced models, we obtain a speed up of about $7.97$ with a relative reduction error of about $10^{-5}$ . In addition, the degradation of the capacity at the voltage $E_{\min}$ over N cycles is shown in Figure 10(c)–(d). This resulted in a speed up of about $23.43$ with a relative reduction error of order $10^{-3}$ .

5. Conclusion

We developed a mathematical model framework for an intercalation battery, consisting of multi-phase porous electrodes, on the basis of non-equilibrium thermodynamics. The framework is very flexible and applicable to a wide range of materials, either for the active intercalation phases or the (liquid) electrolyte, by a stringent formulation of the entire PDE problem in terms of general chemical potential functions. Special emphasis is put on thermodynamic consistency of the transport equations and their respective reaction boundary conditions by employing the very same chemical potential function entirely throughout the model. Periodic homogenisation theory is applied to derive a general set of PDEs for the porous battery cell, where a special scaling of the micro-scale diffusion coefficient leads to a coupled micro-macro scale problem. Spherical symmetry of the intercalation particles is further employed, as well as a 1D approximation of the macro-scale yields an effective 1D + 1D non-linear PDE system. The (dis-)charge current, effectively characterised by the C-rate $C_h$ , enters the PDE system as boundary condition for the electron flux. This allows, on the basis of numerical simulations, the computation of the time- and space-dependent thermodynamic state variables, i.e. the electrolyte potential $\varphi_{\texttt{E}}(x,t)$ , solid potential $\varphi_{\texttt{S}}(x,t)$ , electrolyte concentration $y_{\texttt{E}}(x,t)$ and active phase concentration $y_{\texttt{A}}(x,r,t)$ . Subsequently, this yields important characteristics of an intercalation battery, i.e. the cell voltage E as function of the status of charge $\bar y_{\texttt{A}}$ , parametrically dependent on the C-rate $C_h$ . Further, battery degradation is considered in terms of cycle number n dependent parameters, where exemplarily some degradation models in terms of simple evolution equations were stated. In order to simulate degradation effects, repeated numerical computations of the PDE system are required. For efficient numerical simulations, model reduction techniques were applied to the electrochemical battery model, i.e. the reduced basis method combined with an empirical operator interpolation. We demonstrated the efficient applicability of these method with numerical studies on several ageing scenarios. For degradation effects that impact the diffusion coefficient in the active phase or the intercalation reaction rate, we obtained capacity curves over the number of cycles with a speedup of about 46, compared to full numerical simulations of the same implementation. A speedup factor of about 23 was achieved by additionally investigating the effect of different choices of the charge rate. Numerical relative accuracy of order $10^{-3}$ (at least) was ensured within our simulations.

6. List of symbols ( $j = {\text{An}},{\text{Cat}}$ )