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 [Reference Dreyer, Guhlke and Muller25, Reference Landstorfer48, Reference Newman63] 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 [Reference de Groot and Mazur16, Reference Müller62] $\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 [Reference Dreyer, Guhlke and Müller27, Reference Landstorfer48, Reference Landstorfer, Guhlke and Dreyer52, Reference Müller62].

Electrolyte

For the electrolyte, we consider exclusively the material model [Reference Dreyer, Guhlke and Landstorfer23, Reference Dreyer, Guhlke and Muller24, Reference Landstorfer, Guhlke and Dreyer52] 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 [Reference Landstorfer, Guhlke and Dreyer52])

(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 [Reference Dreyer, Guhke, Landstorfer and Müller21, Reference Dreyer, Guhlke and Landstorfer23, Reference Landstorfer49, Reference Landstorfer, Guhlke and Dreyer52] for its details as well as experimental validation. The bound solvent molecules do not contribute to the entropy of mixing [Reference Dreyer, Guhlke and Landstorfer23], whereby $n_{{\texttt{E}}_S}$ in equation (2.4) denotes the number density of free Footnote ³ 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 [Reference Dreyer, Guhlke and Landstorfer23, Reference Dreyer, Guhlke and Muller24, Reference Landstorfer, Guhlke and Dreyer52] 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,Footnote ⁴ 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 [Reference Landstorfer, Guhlke and Dreyer52] 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 [Reference Landstorfer, Guhlke and Dreyer52]. 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 [Reference Bazant4, Reference Cahn10, Reference Cahn and Hilliard11, Reference Dreyer, Gaberescek, Guhlke, Huth and Jamnik20, Reference Dreyer, Guhlke and Huth22, Reference Garcia, Bishop and Carter34, Reference Landstorfer, Funken and Jacob51, Reference Landstorfer and Jacob53], 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 [Reference Landstorfer and Jacob53]. 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 [Reference Dreyer, Guhlke and Muller25, Reference Dreyer, Guhlke and Muller26, Reference Landstorfer48, Reference Landstorfer50, Reference Landstorfer, Guhlke and Dreyer52, Reference Newman63] 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 [Reference Fuoss32, Reference Fuoss and Onsager33, Reference Liu and Monroe57, Reference Newman, Bennion and Tobias64, Reference Newman and Chapman65, Reference Newman and Thomas66,], 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 potentialFootnote ⁵

(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 [Reference Dreyer, Guhlke and Muller24]

(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 [Reference Fuoss32, Reference Fuoss and Onsager33, Reference Kim and Srinivasan45, Reference Liu and Monroe57, Reference Newman, Bennion and Tobias64] or Nernst–Planck-type flux relations with cross-diffusional terms [Reference de Groot and Mazur16, Reference Latz and Zausch56]. For a three-component mixture $({{\texttt{E}}_S},{{\texttt{E}}_A},{{\texttt{E}}_C})$ , the Onsager reciprocal relations [Reference de Groot and Mazur16] 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 [Reference Dreyer, Guhlke and Muller24, Reference Sanfeld76], 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 electronsFootnote ⁶ ${\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 [Reference Landstorfer50]. 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 [Reference Landstorfer50] 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 [Reference de Groot and Mazur16]).

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 [Reference Doyle, Fuller and Newman19, Reference Newman and Thomas66]. 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}}))$ [Reference Landstorfer50], 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 [Reference Newman and Thomas66], 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 [Reference Newman and Thomas66]).

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 [Reference Landstorfer50]. Hence, the (surface) reaction rate ${\underset{s}{R}\,\!}$ can in general be written with $\alpha \in [0,1]$ as [Reference Defay, Prigogine and Sanfeld17, Reference Dreyer, Guhke, Landstorfer and Müller21, Reference Dreyer, Guhlke and Muller25, Reference Landstorfer48, Reference Sanfeld, Passerone, Ricci and Joud77]

(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 [Reference Landstorfer50]. 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 [Reference Landstorfer50] 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 [Reference Doyle, Fuller and Newman19, Reference Dreyer, Guhlke and Müller27, Reference Hamann and Vielstich38, Reference Landstorfer48, Reference Landstorfer50, Reference Latz and Zausch56, Reference Newman and Thomas66, Reference Sanfeld, Passerone, Ricci and Joud77, Reference Vetter, Bruckenstein, Howard and Technica81]. The experimentally derived Butler–Volmer equation [Reference Hamann and Vielstich38] 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) [Reference Newman and Thomas66], 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}}))$ [Reference Landstorfer50]. 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 [Reference Landstorfer50]. 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 [Reference Doyle, Fuller and Newman19]. However, in the corresponding reaction boundary conditions of [Reference Doyle, Fuller and Newman19] (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 [Reference Newman and Thomas66], 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 [Reference Landstorfer50] 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 [Reference Hunt, Brosa Planella, Theil and Widanage43, Reference Landstorfer48, Reference Landstorfer50, Reference Newman63], which are rather small compared to the intercalation current [Reference Landstorfer50]. 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 [Reference Allaire1, Reference Ciucci and Lai15, Reference Hennessy and Moyles40, Reference Hunt, Brosa Planella, Theil and Widanage43, Reference Marquis, Sulzer, Timms, Please and Chapman59, Reference Richardson, Denuault and Please74]. 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}}}\,$ [Reference Xia, Lu and Ceder83], while the diffusivity of ${Li^+}$ in DMC is in the order of $D_{\texttt{E}} \approx 10^{-5}\,/\,{\frac{\text{cm}}{\,\text{s}}}\,$ [Reference Landesfeind and Gasteiger47]. The macro-length scale is $\mathscr{L} \approx 1/{\unicode{x03BC}} {\rm m}$ [Reference Doyle, Fuller and Newman19], 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 formFootnote ⁷

(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}}$ [Reference Landstorfer, Prifling and Schmidt54] (solid lines in Figure 3). Quite commonly, the (scalar) tortuosity corrector $\tau_{\texttt{E}}$ is introduced via an effective diffusion coefficient [Reference Chung, Ebner, Ely, Wood and García14, Reference Ebner and Wood29, Reference Newman and Thomas66, Reference Tjaden, Cooper, Brett, Kramer and Shearing80], which is in our notation $\tau_{\texttt{E}} = (\pi_{\texttt{E}})^{-1}$ [Reference Landstorfer, Prifling and Schmidt54]. The Bruggeman approximation states that $\tau_{\texttt{E}} = \psi_{\texttt{E}}^{-\alpha}$ [Reference Landstorfer, Prifling and Schmidt54], 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 [Reference Landstorfer, Prifling and Schmidt54], 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 [Reference Landstorfer, Prifling and Schmidt54] 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}

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 [Reference Thomée79], 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 [Reference Quarteroni, Manzoni and Negri73]. 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 [Reference Haasdonk and Ohlberger37]. The POD-Greedy method produces approximation spaces with quasi-optimal $l^{\infty}$ -in- $\mu$ , $l^2$ -in-time reduction error [Reference Haasdonk36].

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. [Reference Ohlberger, Rave and Schindler69]. 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) [Reference Himpe, Leibner and Rave42] and construct an affine approximation of the non-linear differential operator $G_{\mu}$ by an empirical operator interpolation [Reference Chaturantabut and Sorensen13, Reference Drohmann, Haasdonk and Ohlberger28].

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 [Reference Himpe, Leibner and Rave42].

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 [Reference Barrault, Maday, Nguyen and Patera2, Reference Chaturantabut and Sorensen13].

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 [Reference Himpe, Leibner and Rave42]. 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 [Reference Himpe, Leibner and Rave42].

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 [Reference Benaceur, Ehrlacher, Ern and Meunier5, Reference Drohmann, Haasdonk and Ohlberger28]. 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}$ .