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.


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 . 1 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). 2 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.

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.

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.

Domains, definitions and initial scaling
We seek to model a porous electrochemical unit cell, consisting of a porous anode An ⊂ R 3 , a porous separator Sep ⊂ R 3 and a porous cathode Cat ⊂ R 3 , which are packed between two metallic deflectors Def An and Def Cat that yield the overall electrical connections (see Figure 1). The base area of the deflectors is denoted by . The intercalation electrodes j , j = An, Cat consist themself of three phases, namely the active particle phase j A , the conductive additive phase j C and the electrolyte phase j E , with j = i∈{A,C,E} j i . The union of the active phase and the conductive additive is frequently The domains E and A are thus electro-neutral, and we refer to [25,48,63] for details on the derivation. A similar argument holds for the interface CE = C ∩ E . The interface between the deflectors and a single phase of the porous electrode is denoted by i jD = i j ∪ Def i , i = An, Cat, j = A, C, E, for example, m denotes Cat CD the surface through which electrical current from the conductive phase of the electrode flows to the respective deflector.
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 = Cat ∪ Sep ∪ 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 L of the macroscopic width W of the unit cell, i.e. = × [0, W] with ⊂ R 2 being the area of the deflectors and (ii) the length scale of the intercalation particle radii r A . This yields a small parameter ε = L with respect to which we can perform a multi-scale asymptotic expansion and thus a periodic homogenisation.

Variables, chemical potentials and parameter (equilibrium)
The active particle 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 α (x, t), x ∈ i . Further, we denote with the chemical potential of the constituents, which are derived from a free energy density [16,62] [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.
with mole fraction 3) molar concentration n α , total molar concentration of the mixture (with respect to the number of mixing particles [52]) partial molar volume v α of the specific constituent within the mixture, pressure p and reference molar free energy ψ α,ref of the pure substance in the mixture, temperature T and Boltzmann constant k B . A major aspect of the material model is the solvation effect, where each ion binds κ α 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 E S in equation (2.4) denotes the number density of free 3 solvent molecules in the mixture. In turn, the partial molar volume v α of the ionic constituents are increased roughly by κ α · v 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 α,bare and v α,bare , respectively. The corresponding solvated ion binds κ α solvent molecules with mass and volume m E S and v E S , whereby the mass of the solvated ion is m α = m α,bare + κ α m 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 α ≈ v α,bare + κ α v E S . From a thermodynamic point of view, it is convenient and meaningful to assume whereby the incompressibility constraint [23,24,52] implies also a conservation of mass, i.e.
If the molar volume (and mass) of the solvent molecules are far larger than the bare ion, 4 i.e. v E S v α,bare , we may approximate further 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 (κ E A = κ E C =: κ E ). The molar volume of the solvent is related to the mole density n E S ,ref of the pure solvent as 3 Consider a total number of N E S ,tot solvent molecules, N E A anions and N E C cations with respective solvation numbers. Within the mixture, each ion binds κ α solvent molecules whereby the number of free solvent molecules in the mixture is 4 For instance with dimethyl carbonate (DMC) [18] as solvent with molar mass m E S = 90.08 g/mol and lithium as electrolyte cation (m E C ,bare = 6.9 g/mol), we have Note, however, that this assumption only simplifies the equation system a bit and is not a necessary assumption. and the solvation numbers κ α 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 α in terms of the mole fractions y α , i.e.
Finally, the electrostatic potential is an additional variable in the electrolyte phase and denoted as ϕ E (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 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 μ A C = μ A C (y A ). The chemical potential of intercalated lithium is stated as with mole fraction of intercalated cations in the active phase and partial molar Gibbs energy g A C = const. The number density n A,lat of lattice sites is constant, which corresponds to an incompressible lattice, and the enthalpy parameter γ A > −2.5. Note that γ A < −2.5 entails a phase separation [53]. The electrostatic potential in the solid phase is denoted as ϕ S (x, t).

Electro-neutrality
The electro-neutrality condition of A , E and C can be obtained by an asymptotic expansion of the balance equations in the electrochemical double layer at the respective surface . 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 AE and (iii) that for monovalent electrolytes the cation mole fraction (or number density) is equal to the anion mole fraction, i.e.
The total number density n E,tot = n E S + n E C + n E A and the electrolyte concentration n E C can be written as a function of y E using equation (2.6) as follows: (2.14)

Transport equation relations
In the electrolyte E , we have two balance equations determining the concentration n E C (x, t) (or mole fraction y E (x, t)) and the electrostatic potential ϕ E (x, t) in the electrolyte [32, 33, 57, 64, 65, 66,], i.e.
For the numerical simulations of Section 3, we assume constant values for the transport parameters (t E C , S E , D E , E ), which is sufficient to show the general methodology of our framework approach. Note, 5 Note that this arises from the fact that the diffusional flux of cations is [24,27,50] and since however, that (t E C , S E , D E , E ) can in general be depend on the electrolyte concentration n E C , which can be easily incorporated our framework.
In the active particle A , we have two balance equations determining the concentration n A C (x, t) (or mole fraction y A ) of intercalated lithium and the electrostatic potential ϕ S (x, t) in the solid phase S , i.e. the balance equation for intercalated lithium A C and the balance equation for electrons 6 A e in the entire solid phase, and (dimensionless) thermodynamic factor The (solid state) diffusion coefficient D A is further assumed to be where the term (1 − y A ) accounts for a reduced (solid state) diffusivity due to crowding [50]. Note that the electrical conductivity σ S is different in the active phase A and the conductive phase C , which form S = A ∪ C . We account for this as In principle, σ A can be dependent on the amount of intercalated ions, i.e. σ A = σ A (y A ), but for the sake of this work we assume σ A = const. and σ C = 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 μ A (y 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 [16]). For the very special case f A = ln yA 1−yA , i.e. an ideal lattice gas, one obtains tf A = 1 1−yA and thus together with equation (2.31) a simple, linear diffusion equation 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 This is, however, in contrast to the OCP relation of [66], in our notation E OCP = E 0 − k B T e 0 ln yA 1−yA + βy A + ζ (β and ζ are considered as parameters in [66]).

Reaction rate and affinity
At the the interface AE the intercalation reaction occurs, which is modelled on the basis of (surface) non-equilibrium thermodynamics [50]. Hence, the (surface) reaction rate R s can in general be written with α ∈ [0, 1] as [17,21,25,48,77] 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φ E := ϕ E − E Li + ,E denotes the electrolyte potential vs. metallic Li andφ j S := ϕ j S − E j Li + ,A the electrode j potential vs. metallic Li [50]. Note again that y A | AE denotes the evaluation of y A at the interface A,E and that the surface affinity (2.36) is dependent on the chemical potential (or the mole fraction) evaluated at the interface. The nonnegative function L s in (2.35) ensures a non-negative entropy production r s σ ,R due to reactions on the surface, i.e. r For the sake of this work, we assume L s = 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.
. For intercalation materials, this OCP is concentration dependent, i.e. E OCP = E OCP (y A ), which arises from the chemical potential function μ A (y A ) since in a half cell vs. metallic lithium E OCP = 1 e 0 (μ Li − μ A (y 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 R s 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 μ 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 μ A = g A C + k B T ln yA 1−yA , used for example in [19]. However, in the corresponding reaction boundary conditions of [19] (Butler-Volmer equation), the chemical potential function (in our notation) μ A = g A C + frack B T e 0 ln yA and ζ 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 L s in our formulation (2.35). Albeit good results on experimental validation, this model is not within our framework as different chemical potential functions μ A (y 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 μ A (y A ).

Balance equations
Applying the above modelling procedure for j = An, Cat yields the following equation system, Note that σ j S (x) incorporates the fact that the conductivity is far larger in the conductive additive phase j C then in the active particle phase j A . 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 j where by convention n is pointing from j A into j E . 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 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.

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

Introduction
An important feature of the coupled transport equation system (2.37)-(2.40) is the circumstance that the solid-state diffusion D j A is far smaller than the electrolyte diffusivity D E . This accompanies, however, with smaller diffusion pathways on the length scale of the intercalation particles. The diffusivity of Li in LiCoO_2 is, for instance, about D A ≈ 10 −12 / cm s [83], while the diffusivity of Li + in DMC is in the order of D E ≈ 10 −5 / cm s [47]. The macro-length scale is L ≈ 1/μm [19], while the micro-scale is ≈ 1 / nm (see Figure 2). Hence, the length scale ration L = ε ≈ 10 −3 and D j A ≈ ε 2 D E . This motivates the re-scaling which is subsequently exploited in the homogenisation procedure.

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 For the sake of the homogenisation, it is convenient to introduce the scaling

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 7 where u i , i = A, E, S, denotes the respective phase variable, (D ε i ,R ε i ) the ε-dependent non-equilibrium parameter and h i capture the non-linearities w.r.t. u i . For the solid-phase S , the x dependency of h S u S , x ε accounts for different conductivities in C ⊂ S and A ⊂ S , e.g. (2.63) Note that we abbreviate We consider a multi-scale expansion (i = A, E, S) For non-linear functions h = h(u), we consider the ε−Taylor expansion We abbreviate and expand thus Insertion of the multi-scale expansion yields a sequence of PDEs in the orders of ε.
Order ε −2 : Order ε −1 : This yields where and Note that n k denotes the kth component of the corresponding normal vector n on the surface σ AE .

Equation (2.80) leads (by an integration over ω E and integration by parts) for
and that u 0 A | ∂ωA denotes an evaluation of u 0 A at the boundary of the micro-domain ω A . Equation (2.80) leads (by an integration over ω S ) for i = S to Hence we obtain for the equation system (2.56)-(2.62) via periodic homogenisation the coupled micromacro balance equation system

Spherical particles
In the following, we assume on the macro-scale a 1D approximation x ∈ [0, 1] as well as spherical 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 ω), the porous media parameters (ψ E , π E , θ AE ) and (π S , ψ S ) can (numerically) be computed by solving the cell problems (2.78) and (2.79), respectively. Treating the particle radii r A as a parameter yields then the diffusion corrector π E as function of the porosity ψ E [54] (solid lines in Figure 3). Quite commonly, the (scalar) tortuosity corrector τ E is introduced via an effective diffusion coefficient [14,29,66,80], which is in our notation τ E = (π E ) −1 [54]. The Bruggeman approximation states that τ E = ψ −α E [54], where α is to be determined (Bruggeman fit, dashed lines in Figure 3).
Note that for equally sized spherical particles the actual micro-structure is exclusively encoded in the porous media parameters (ψ E , π E , θ AE , π S , ψ 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.  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).

Homogenised battery model
where and boundary conditions Note that t E C = const. allows us to rewrite (2.103) with (2.104) as 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). the cell voltage (2). the basic (electrode) capacity (j = An, Cat) . the present (electrode) capacity (j = An, Cat) 113) (5). and the 1-C current density which yields for the current density the scaling where C h ∈ R is the C-Rate.
Note that For a constant current discharge (i = const.), we obtain thus

Scalings
Summarised, we use the scalings as well as The boundary conditions read 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 λ : where (+1) is the stoichiometric coefficient of the product Li.

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

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

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 microstructure 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 (φ E ,φ S ) and the mole fraction of 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 ν. The overall PDE system can be written in the following abstract form: The linear and nonlinear coefficient functions α, β, γ , δ, κ, ω and ρ 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 β(·, y A )∇ ν y A · n = −R 4 (·,φ E ,φ S , y E , y A ) at the boundary of the electrode particles, i.e. at ν = 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.

Discretisation
For the battery model in the abstract form above, let 1D denote a computational one-dimensional domain in the macroscopic direction and ν = ξ ν the transverse/radial directions in the electrode particles associated with each ξ ∈ 1D . Let P ⊂ R P , P 1 denote the parameter space. We define the solution 3 and V , 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: where the operator G μ (·, ·) : V × V → R represents the non-linear time-discrete PDE system. The index μ ∈ P indicates the dependence of the problem on certain parameters, such as the charge rate, the diffusion coefficient or the reaction rate. f μ (u t , ·) ∈ V 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 ⊂ V. We hence obtain for each time step a fully discrete non-linear system of the form: where ψ i , i = 1, . . . , n denotes the standard Lagrange basis of the finite element space V h = 4 i=1 V i h . Henceforth, the operator G μ can be called the finite element operator which operates on V h × V h . The developed discretisation does not depend on the specific choice of the parameters to be varied.

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 V ⊂ V h that approximates this solution manifold in order to accelerate the repeated solution of (3.2) for varying parameters μ. Under this projection, the reduced problem is given by where m n andψ i , i = 1, . . . , m, represents the basis of the reduced (basis) spaceṼ. 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 lowdimensional 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 ∞ -in-μ, 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 P = {μ 1 , . . . , μ ns } be a set of n s parameter samples and {u h (μ 1 ), . . . , u h (μ ns )} the corresponding snapshot set. At each time step, the snapshot of the set is calculated using the following prescription of Newton's method: We separate the snapshots into the respective components, u 1 h , . . . , u 4 h and generate the reduced basis separately. We define the corresponding snapshot matrices S i ∈ R N ι h ·|t|×ns with i = 1, 2, 3, 4 and ι ∈ { 1D , ν } as where the vectors u j i ∈ R N ι h ·|t| , 1 j n s , denote the degrees of freedom of the functions span the reduced spaceṼ i using only the singular vectors whose singular values are above a fixed threshold value ε. Due to the fundamental properties of POD, the projection error consist of the l 2sum of the corresponding truncated singular values, andṼ 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Ṽ := 4 i=1Ṽ i .

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 μ by an empirical operator interpolation [13,28].
In detail, the offline-online decomposition consists of precomputing parameter-and solutionindependent terms, a so called collateral basis that allows to interpolate evaluations of the operator G μ . 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 μ ) 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:

4)
where {G q } M q=1 is the collateral basis, i.e. a basis for a subspace of and θ q μ are interpolation coefficients recalculated for each μ andũ 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 μ (u h , ·). The set of snapshots thereby includes the Newton stages in addition to the corresponding solution trajectory of G μ (u h , ·). 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: where the vectors G k μ j (u i ) ∈ R N ι h , 1 j n s , 1 k s j , represent the degrees of freedom of the functions G μ j u k i h (μ j ) and POD the error tolerance for the POD. We define G 1 , . . . , To calculate the interpolation coefficients θ q μ (ũ), q = 1, . . . , M for given μ andũ, the interpolation constraints are imposed at M interpolation points. The interpolation points are selected iteratively from the indices of basis {G q } M q=1 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 μ in (3.3) by the fast-to-evaluate interpolated operator I M (G μ ), we obtain the completely offline-online decomposable reduced problem which can be solved efficiently for varying parameters μ. 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 μ 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 β 1 , · · · , β B . Starting with two leaves β 1 , β 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 α 1 , which is again used to calculate a POD. This newly generated input and the local snapshot data assigned to leaf β 3 are the input of the parent node α 2 . The final HAPOD modes ρ are reached, when the last leaf β B has entered. For the calculation of the collateral basis, we e.g. define for i = 1, 2, 3, 4: where POD is the desired approximation error tolerance for the resulting HAPOD space. Depending on ω, one might get more modes than needed for a POD with the same tolerance. The local tolerances in the HAPOD algorithm are computed from POD and ω ∈ (0, 1). More details can be found in [42].

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, . . . , 1000 for different ageing models. The electrochemical characteristics are the voltage-capacity spectrum E(ȳ A , n; C h ) and the status of chargeȳ Cat A 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 rateL(n) and diffusion coefficientD 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 rateL 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 coefficientD 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 ofD A andL by alternately choosing one of the parameters fixed. In the last case, we vary all parametersD A ,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.

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 1D h = 300 grid points. Here, each of the cathode, separator and anode is discretised into 100 grid points. The pseudodimension for the microscopic equation, which goes vertically at each bottom grid point, consists of N ν h grid points. In the simulations N ν h = 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 ν = 1. To ensure that the bottom grid corresponds to ν = 1 a transformation of the microscopic equation withν = (1 − ν) 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: y A (g) = e g 1 + e g , g(y A ) = ln The time discretisation is performed by an implicit Euler method on a T = 1 time interval with a time step size of 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 ξ ∈ Catφ S (ξ ) E min with the voltage value E min = −0.2. To generate the reduced spaceṼ, we compute a snapshot set S train on training sets of equidistant parameters.
For the experiments 1 and 3, we choose |P train | = 15 and for experiment 2 we choose |P train | = 10.
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 P test given by For the truth solution, we want to ensure that we have achieved grid convergence up to a certain error ε. Thus, we assume that the high-dimensional discretisation adopts a resolution of ||u n − u m || 2 < ε 10 −5 with 30,900 grid points for u n and 31,512 grid points N 1D 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 μ 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 = 4e − 8 and ω = 0.9 in both cases. For illustration purposes, the reduced space and empirical interpolation are calculated for a training set |P train | = 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.

Experiment 1
In this subsection, we consider the variation of the charge rate C h with C h ∈ [0.01, 4] for an unaged battery. In this case, we chooseL = 0.5 andD 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 ).
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 ∈ {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.

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 parametersL andD A equally in anode and cathode over n with n = 0, . . . , 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 parametersL(n) andD A (n) as a function of the number of cycles n satisfies an ordinary differential equation with the unknown parameter a F (β) such that: It follows that F(n) = F 0 e log (β)n /N and a F (β) = log (β) N (see Figure 7). F 0 is the corresponding initial parameter value. In our case,D A,0 =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 rateL forD A = 0.5 is shown in Figure 8(a), (c) for two different choices of β. Furthermore, the capacityȳ Cat A at the specific voltage value E min = −0.2 over the number of cycles N = 1000 with a variation of β is illustrated in Figure 8(e). The same scenario is shown for the degradation of the diffusion coefficientD A forL = 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 ofL has a larger effect in terms of absolute capacity loss.
For the implementation of the reduced order model (ROM) with dependence on the parameters μ = [D A ,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

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 parametersL andD 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 μ = [C h ,D A ,L]. Furthermore, as in Experiment 2, we assume that the evolution of the parametersD A andL as a function of the number of cycles n satisfies the ordinary differential equation that additionally depends on the charge rate C h with the unknown parameter a F (β) such that: It follows that F(n) = F 0 e C h log(β)n/N and a F (β) = log (β) N (see Figure 9). Here we set β = 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.
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 Figure 9. Various degradation models that satisfy the ordinary differential equation (4.4) with F 0 = 0.5, β = 0.6 and N = 1000. 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 .

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 ϕ E (x, t), solid potential ϕ S (x, t), electrolyte concentration y E (x, t) and active phase concentration y 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ȳ 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.

List of symbols (j =
whereby only two of the transport parameters (t E C , S E , D E , E ) are independent and all of them are constant.

Appendix B: Parameters
The following