MULTIBAT: Unified workflow for fast electrochemical 3D simulations of lithium-ion cells combining virtual stochastic microstructures, electrochemical degradation models and model order reduction

TL;DR

This paper introduces a comprehensive workflow combining stochastic microstructure generation, detailed electrochemical modeling, and model order reduction to enable efficient 3D simulations of lithium-ion batteries, focusing on microstructure-performance interactions.

Contribution

It presents a unified simulation workflow that integrates microstructure modeling, detailed electrochemical physics, and model order reduction for fast 3D battery simulations.

Findings

Microstructure generation reduces reliance on tomographic imaging.

The integrated model captures lithium plating and stripping effects.

Model order reduction significantly decreases computation time.

Abstract

We present a simulation workflow for efficient investigations of the interplay between 3D lithium-ion electrode microstructures and electrochemical performance, with emphasis on lithium plating. Our approach addresses several challenges. First, the 3D microstructures of porous electrodes are generated by a parametric stochastic model, in order to significantly reduce the necessity of tomographic imaging. Secondly, we integrate a consistent microscopic, 3D spatially-resolved physical model for the electrochemical behavior of the lithium-ion cells taking lithium plating and stripping into account. This highly non-linear mathematical model is solved numerically on the complex 3D microstructures to compute the transient cell behavior. Due to the complexity of the model and the considerable size of realistic microstructures even a single charging cycle of the battery requires several hours computing time. This renders large scale parameter studies extremely time consuming. Hence, we develop a mathematical model order reduction scheme. We demonstrate how these aspects are integrated into one unified workflow, which is a step towards computer aided engineering for the development of more efficient lithium-ion cells.

Tables

Table 1. Table 1: Overview of interface conditions between the different material domains for the ion fluxes N 𝑁 N and current densities j 𝑗 j . “cont.” mathematically means no interface but a continuous flux according to the transport equation. Between the graphite electrode and the lithium foil as counter electrode there is obviously no interface.

	graphite	electrolyte	lithium foil	current collector	plated lithium
graphite	cont.	(2.6),(2.7)	no interface	$N=0,j=\text{cont.}$	$N=0,j=\text{cont.}$
electrolyte	(2.6),(2.7)	cont.	(2.6),(2.8)	$N=-1,j=0$	(2.6),(2.14)
lithium foil	no interface	(2.6),(2.8)	cont.	$N=0,j=\text{cont.}$	no interface
​​ current col.	$N=0,j=\text{cont.}$​​​	$N=0,j=0$​​​	$N=0,j=\text{cont.}$ ​​​​​	cont.	$N=0,j=\text{cont.}$​​​
​​ plated lithium	$N=0,j=\text{cont.}$	(2.6),(2.14)	no interface	$N=0,j=\text{cont.}$​​​	cont.

Table 2. Table 2: Extrapolated timings for the model reduction experiment. Time for single full model simulation: 15h 38m 35s (median), time for single reduced simulation: 7m 48s (median), time for generation of reduced model from snapshot data: 13h 43m 12s. ‘without MOR’ is the required time if all simulations are performed with the full order model ( 2.4 ), ‘with MOR’ is the required time if the reduced order model ( 2.5 ) is used for all simulations after the first 3 snapshot computations (including reduced order model construction). All computations have been performed on a single core of an Intel Xeon E5-2698 v3 CPU.

simulations	without MOR			with MOR			speedup
1		15h	38m	–			–
2	1d	7h	17m	–			–
3	1d	22h	55m	–			–
4	2d	14h	34m	2d	12h	46m	1.0
5	3d	6h	12m	2d	12h	54m	1.3
10	6d	12h	25m	2d	13h	33m	2.5
50	32d	14h	9m	2d	18h	45m	11.7
100	65d	4h	18m	3d	1h	15m	21.4
limit (= full model vs. reduced order model)							120.3

Equations

$${\partial }_t{c}_{\mathrm{Gr}}$$

$$0$$

$${\partial }_t{c}_{\mathrm{El}}$$

$$0$$

$$\frac{\partial \mu }{\partial c_{\mathrm{El}}}$$

$$\begin{array}{rl}{\mathbf{j}}_{\mathrm{El}}\cdot {\mathbf{n}}_{\mathrm{So}-\mathrm{El}}& =i_{\mathrm{interface}},\\ {\mathbf{j}}_{\mathrm{So}}\cdot {\mathbf{n}}_{\mathrm{So}-\mathrm{El}}& =i_{\mathrm{interface}},\\ {\mathbf{N}}_{\mathrm{El}}\cdot {\mathbf{n}}_{\mathrm{So}-\mathrm{El}}& =\frac{i_{\mathrm{interface}}}{F},\\ {\mathbf{N}}_{\mathrm{So}}\cdot {\mathbf{n}}_{\mathrm{So}-\mathrm{El}}& =\frac{i_{\mathrm{interface}}}{F},\end{array}$$

$$i_{\mathrm{interface}}$$

$$=2\cdot i_{\mathrm{Gr}-\mathrm{El}}^{00}\cdot \sqrt{c_{\mathrm{Gr}}\cdot c_{\mathrm{El}}}\cdot \sinh \left(\frac{F}{2\cdot R\cdot T}\cdot {\eta }_{\mathrm{Gr}-\mathrm{El}}\right),$$

$$i_{\mathrm{interface}}=i_{\mathrm{CE}-\mathrm{El}}$$

$${\mathrm{Li}}_{\mathrm{Electrolyte}}^{+}+{\mathrm{e}}_{\mathrm{Electrode}}^{-}\rightleftharpoons {\mathrm{Li}}_{\mathrm{Plated}}^0.$$

$$F\cdot {\eta }_{Pl/St}={\tilde{\mu }}_{\mathrm{L}{\mathrm{i}}^{+}}^{\mathrm{Pl}}-{\tilde{\mu }}_{\mathrm{L}{\mathrm{i}}^{+}}^{\mathrm{El}}.$$

$${\eta }_{Pl/St}={\Phi }_{\mathrm{Pl}}-{\phi }_{\mathrm{L}{\mathrm{i}}^{+}}^{\mathrm{El}},$$

$$i_{\mathrm{Pl}-\mathrm{El}}=$$

$$\cdot \left(f_{\mathrm{pre}}\left(n_{\mathrm{Li}}\right)\cdot \exp \left(\frac{F{\eta }_{\mathrm{Pl}-\mathrm{El}}}{2\cdot R\cdot T}\right)-\exp \left(-\frac{F{\eta }_{\mathrm{Pl}-\mathrm{El}}}{2\cdot R\cdot T}\right)\right).$$

$$f_{\mathrm{pre}}\left(n_{\mathrm{Li}}\right)$$

$$i_{\mathrm{interface}}$$

$$\left[\begin{array}{c}\frac{1}{\Delta t^{(n+1)}}(c_{\mu }^{(n+1)}-c_{\mu }^{(n)})\\ 0\end{array}\right]+A_{\mu }\left(\left[\begin{array}{c}{c}_{\mu }^{(n+1)}\\ {\phi }_{\mu }^{(n+1)}\end{array}\right]\right)=0,$$

$$(c_{\mu }^{(n+1)},{\phi }_{\mu }^{(n+1)})\in V_h\oplus V_h.$$

$$P_{\tilde{V}}\left\{\left[\begin{array}{c}\frac{1}{\Delta t^{(n+1)}}({\tilde{c}}_{\mu }^{(n+1)}-{\tilde{c}}_{\mu }^{(n)})\\ 0\end{array}\right]+A_{\mu }\left(\left[\begin{array}{c}{\tilde{c}}_{\mu }^{(n+1)}\\ {\tilde{\phi }}_{\mu }^{(n+1)}\end{array}\right]\right)\right\}=0,$$

$$({\tilde{c}}_{\mu }^{(n+1)},{\tilde{\phi }}_{\mu }^{(n+1)})\in \tilde{V},$$

$${\mathcal{I}}_M[O](x)\in \tilde{Y}\text{and}{\pi }_m({\mathcal{I}}_M[O](x))={\pi }_m(O(x)),1\le m\le M.$$

$$A_{\mu }=A_{\mu }^{(aff)}+A^{(bv)}+A^{(1/c)},$$

$$A_{\mu }^{(aff)}=A^{(const)}+\mu \cdot A^{(bnd)}+A^{(lin)},$$

$$A_{\mu }\approx {\tilde{A}}_{\mu }=A_{\mu }^{(aff)}+{\mathcal{I}}_{M^{(bv)}}[A^{(bv)}]+{\mathcal{I}}_{M^{(1/c)}}[A^{(1/c)}].$$

$$P_{\tilde{V}}\left\{\left[\begin{array}{c}\frac{1}{\Delta t^{(n+1)}}({\tilde{c}}_{\mu }^{(n+1)}-{\tilde{c}}_{\mu }^{(n)})\\ 0\end{array}\right]+{\tilde{A}}_{\mu }\left(\left[\begin{array}{c}{\tilde{c}}_{\mu }^{(n+1)}\\ {\tilde{\phi }}_{\mu }^{(n+1)}\end{array}\right]\right)\right\}=0,$$

$$({\tilde{c}}_{\mu }^{(n+1)},{\tilde{\phi }}_{\mu }^{(n+1)})\in \tilde{V},$$

$$\Vert {\hat{c}}_{\mu }^{(n)}-c_{\mu }^{(n)}\Vert \le \Theta \cdot \Vert {\tilde{c}}_{\mu }^{(n)}-c_{\mu }^{(n)}\Vert ,\Vert {\hat{\phi }}_{\mu }^{(n)}-{\phi }_{\mu }^{(n)}\Vert \le \Theta \cdot \Vert {\tilde{\phi }}_{\mu }^{(n)}-{\phi }_{\mu }^{(n)}\Vert$$

$$\Vert {\tilde{c}}_{\mu }^{(n)}-c_{\mu }^{(n)}\Vert \le \frac{1}{1-\Theta }\cdot \Vert {\tilde{c}}_{\mu }^{(n)}-{\hat{c}}_{\mu }^{(n)}\Vert ,\Vert {\tilde{\phi }}_{\mu }^{(n)}-{\phi }_{\mu }^{(n)}\Vert \le \frac{1}{1-\Theta }\cdot \Vert {\tilde{\phi }}_{\mu }^{(n)}-{\hat{\phi }}_{\mu }^{(n)}\Vert .$$

Full text

1 Introduction

The ubiquity and importance of rechargeable lithium-ion batteries lead to the increasing demand for physics-based simulation methods that are able to analyze and predict battery behavior. These methods can not only contribute in improving cell design and operation, but they can also greatly support battery research in its understanding of basic mechanisms, like lithium plating that determine battery life and safety, which is yet not well understood.

The electrochemical simulation of lithium-ion cells goes back to the work of Newman and his co-workers [9, 16, 43]. Their simulation methodology is based on the porous electrode theory developed by Newman [45]. This model approach neglects the details of electrode microstructures and describes them as a homogeneous medium where electrolyte and the solid material coexist at every point. The most commonly used model of Newman only considers the through-direction of the battery. It takes into account the diffusion of lithium ions into the active material by assuming a spherical, microscopic particle of average size in each discretization point in which a one-dimensional diffusion equation is solved. Hence this model is sometimes called a pseudo-2d (P2D) model [55]. There exist many applications for Newman-type models like the study of cell behavior as well as degradation [55, 3, 54, 61]. But the main drawback of these models is that the complex electrode microstructures are only approximately accounted for by a few aggregated parameters: the thickness of the electrode $L$, the porosity $\varepsilon$, the mean particle radius $r$ and the specific interface area between electrolyte and active material $a$ [9]. Furthermore, effective transport parameters need to be determined to account for the influence of the microstructure on the average species transport. While these models are able to describe the average battery behavior surprisingly well [61, 10, 12, 33], they cannot be expected to capture local microscopic effects. In particular, many degradation effects like, for instance, lithium plating depend on the local environment. Hence homogenized models cannot fully capture the interplay between microstructure and degradation phenomena with sufficient predictive power. Therefore, more fundamental, spatially resolved models should be applied that are able to take the electrode microstructure explicitly into account [30, 31]. Without simplifications like volume averaging for the P2D-models these allow the computation of quantities on the scale of the electrode microstructure and are hence better suited for plating predictions. To give an example, in [57] a microstructure-based simulation study for a LCO-graphite battery was performed concentrating on the discharge behavior for a 2D cut of one given realization of the electrodes. Although the numerical solution of these micro-models is computationally much more demanding they have been successfully applied to study cell performance [36, 64, 25], coupling to thermal effects [33, 68, 60], and to account for phase-separation dynamics within certain electrode materials [24]. A framework for these spatially resolved simulations has been implemented in the software BEST [15].

Lithium plating is one of the major degradation factors and security risks in lithium-ion batteries. Lithium plating describes the deposition of metallic lithium on the negative electrode [63]. This causes a loss of usable lithium (which reduces the cell’s capacity) and might lead to the growth of lithium dendrites which can eventually create a short-circuit between the electrodes which can favor catastrophic thermal runaways. While model extensions to account for lithium plating are typically based on the porous electrode theory [61, 46, 2, 59, 35] only very little work has been published where lithium plating models take the electrode microstructure into account [41, 1]. In a recent publication [22] a micro-scale model has been developed that is able to take the inhomogeneous electrode structure into account. The research presented in the current paper is based on this degradation model.

Spatially resolved electrochemical simulations as described above allow investigations of electrochemical behavior for realistic 3D microstructures. Thus, as input for these simulations, realistic 3D image data of battery electrodes is needed, which is already available even in-operando [14]. However, tomographic measurements of battery electrodes in 3D involve high costs and efforts. A methodology that has proven to be very promising in this context is stochastic microstructure modeling. Based on (only one or a few) tomographic measurements, a 3D parametric stochastic microstructure model can be constructed and calibrated using tools from stochastic geometry [7]. The model has been implemented in our software library GEOSTOCH [38]. It is able to describe the complex geometric microstructure in a statistical sense with only a few parameters such that each realization of the model represents the morphological characteristics of the tomographic image data (e.g., the distributions of particle size and shape, pore size distribution, etc.). Once fitted to tomographic image data, with hardly any effort an arbitrary number of realistic 3D microstructures can be generated on the computer. Moreover, systematic variation of model parameters allows the realization of virtual, but still realistic microstructures. Such an approach has been considered, for example, in the context of organic solar cells [58]. Using regression in the parameter space, microstructures that represent various manufacturing conditions could be generated on the computer and analyzed regarding their functionality. This results in an enormous reduction of complexity, as (most of) the structures do not have to be manufactured in the laboratory, but only tomographic image data of a few ones is needed. Similar examples of stochastic microstructure modeling can be found in literature [17, 67, 42]. In this work, a stochastic microstructure model for anodes in lithium-ion batteries [13] is used.

While the aforementioned microscopic battery model can be solved by relatively standard iterative numerical methods, the solution process is computationally very demanding. In order to get meaningful results a sufficiently large electrode cutout needs to be resolved in the simulation. This results in huge time-dependent discrete systems which require considerable computing resources, already for single simulation runs. Computational studies to identify critical parameters, to estimate the dependence of degradation on operating conditions or to support optimal design and control of batteries, however, require many forward simulation runs with varying material or state parameters and are thus virtually impossible. Hence, model reduction approaches for the resulting parameterized systems are indispensable for such simulation tasks. Concerning model reduction for lithium-ion battery models, we refer to the pioneering work [6] in the context of proper orthogonal decomposition (POD), and to the more recent contributions [37, 27, 65, 29, 49, 48, 47] in the context of reduced basis methods. In the work presented here, we rely on an implementation of recent model reduction methods (such as the reduced basis method, POD, and the empirical interpolation method) implemented in our model order reduction library pyMOR [39, 52].

The ability to efficiently and realistically predict the degradation behavior (here: lithium plating) of lithium-ion batteries under arbitrary load conditions relies on the following prerequisites:

A physics-based predictive microscopic battery model that includes the plating mechanism. 2. 2.

A method to create a number of virtual, yet realistic microstructures as basis to understand the correlation between structural properties and battery performance and degradation behavior. 3. 3.

A numerical method that is able to efficiently perform a considerable number of three-dimensional, microstructure-resolving simulations for a variety of operating conditions. 4. 4.

A software interface that is able to integrate these aspects into a common workflow.

Within the project MULTIBAT [26] the authors developed and technically implemented a workflow that covers all the aforementioned aspects, namely stochastic geometry generation, model extension to account for plating, numerical implementation and development of model order reduction techniques. There are numerous papers in the literature (including several ones written by the authors of the present paper) on different components of the presented workflow. However, we are not aware of any publication on an algorithm integrating all these components into one single, holistic workflow, which enables comprehensive solutions of really complex problems related to Li-ion batteries. Thus, development, implementation, and testing of a holistic algorithm / workflow which integrates all components of the above-mentioned chain is one of the main contributions of the present paper. A feasibility study for a really complex problem, such as plating, is presented in order to illustrate the capabilities of the workflow that has been developed. The investigation of the interplay between 3D microstructure and electrochemical processes during the plating processes, which up to our knowledge has not been done so far in the literature, is another main contribution of the present paper.

The authors’ developments on the individual components of the workflow have been reported earlier, details can be found in the references listed in the present paper, hence, these components are presented here relatively shortly. The emphasis in the present paper is given to the developed interfaces, to the integration of all the components into one single workflow, to the peculiarities related to the selected feasibility study, and to the parametric study of the plating process in stochastic geometry generation. Special attention is paid to computational efficiency, adapting the model reduction approach to the heavily nonlinear system of partial differential equations. The presented study reveals, that the complex information produced by the interplay between microstructure, lithium-ion transport and intercalation kinetics is hidden in a vastly reduced subspace of the full 3D information contained in the time-dependent scalar fields for lithium-ion concentration and electrochemical potential. The essential dynamics leading to plating in a complex microstructure can be represented by a sophisticated reduced model without losing spatial precision. The advantage of the model presented in this paper is the ability to perform fully 3D microstructure-resolved simulations of plating with nearly the same numerical efficiency as simulation with a P2D model i.e. a 1D volume averaged battery model, in which all structural details are lost. We report on the MULTIBAT workflow and briefly describe the details of all the individual aspects in Section 2. In Section 3 we demonstrate the application of the developed methods by showing and discussing results of a simulation study and conclude with a summary in Section 4.

2 The MULTIBAT workflow

In this section we discuss the individual components of the MULTIBAT workflow (see Fig. 1) and their realization in more detail. Based on experimental data, random electrode geometries with the same or modified structural characteristics are generated (Sec. 2.1), and a mathematical model of the relevant electrochemical effects is formulated (Sec. 2.2 and 2.3). The resulting continuum model is then discretized (Sec. 2.4) and reduced (Sec. 2.5), leading to a quickly computable microscale model of the cell dynamics on realistic electrode geometries. The software implementation and integration into a unified modeling and simulation workflow is discussed in Sec. 2.6.

It should be noted that, while we present a specific realization of the MULTIBAT workflow targeted at lithium plating, the same workflow can be applied to other questions in electrochemistry and similar problem domains. Each individual component can be further developed and optimized for other specific applications, independently of the other workflow components.

2.1 Generation of random structures

The study of local effects in the complex microstructures of battery anodes by spatially resolved models is computationally very expensive, particularly regarding random access memory. Therefore, only quite small sample sizes can be considered. As we are interested in local phenomena, there is a need for high-resolution of the images, which on the other hand means that the images typically only represent small cutouts of the material. This is why, in order to get reliable results, the computer experiments have to be carried out repeatedly using different samples. Furthermore, the imaging techniques are complex in preparation and involve long imaging times as well as high costs. This is why a suitable approach is to use randomly generated images of microstructures using parametric stochastic 3D models. This approach has already been used successfully in related applications for energy materials in fuel cells [18] and solar cells [58]. A parametric stochastic model which describes the spatial structure is developed for the material and its parameters are fitted to image data. Using the calibrated model, an arbitrary number of structures that are similar to the image data in a statistical sense can be generated with hardly any effort. ‘Similar in a statistical sense’ means that the realizations of the model do not resemble the image data exactly, but with respect to aggregated quantities and spatial properties. For example, simple characteristics like volume fraction and specific surface area can be matched, but also more complex spatial characteristics like the distribution of pore sizes or local tortuosity. Thus, realizations of a parametric stochastic microstructure model are an ideal input for spatially resolved electrochemical simulations. A further advantage is that their parameters can be changed to create virtual structures that have not been produced in the laboratory yet, and the electrochemical performance of those virtual structures can be analyzed on the computer, a procedure called virtual materials testing.

Here, we make use of a parametric stochastic 3D microstructure model for anode structures from lithium-ion battery cells [13]. Besides the validation based on structural characteristics [13] a validation using spatially resolved electrochemical simulations has been performed [21]. The variability of the modeling approach used here is demonstrated since the same model with some adaptions can be used to generate microstructures for energy cells [13] and power cells [66].

We now briefly recall some details of the stochastic 3D model that is used to generate the virtual anode microstructures used in the MULTIBAT workflow. As mentioned above the model has already been published [13] and all parameters as well as further details can be found there. Generally, the construction of the model consists of four steps that are also depicted in Fig. 2.

First, the locations, sizes and shapes of the particles are determined. Technically speaking, a Laguerre tessellation is generated (see Fig. 2) based on a random sequential adsorption process. This tessellation decomposes the region of interest into convex polytopes. Later on, a particle is placed inside each of these polytopes. Thus, the Laguerre tessellation roughly indicates the spatial arrangement of particles. For details regarding tessellations, a broad spectrum of literature is available [7, 34, 40, 62].

In the next step, a connectivity graph is generated that describes which particles are supposed to be connected, i.e., for each polytope $P$, we determine a set of neighboring polytopes $\{P_{i},\,i=1,...,N\},\,N\in\mathbb{N}$. The particles that are placed inside $\{P_{i},\,i=1,...,N\}$ have to touch the particle in $P$. Full connectivity of all particles is ensured by the usage of a minimum spanning tree [51]. Further connections are added to the minimum spanning tree depending on the size of the facet between two polytopes, as the probability of two particles being connected is larger for larger facet areas. Such a graph is depicted in Fig. 2.

Now, a particle can be realized in each polytope, fulfilling the boundary conditions, i.e., touching the particles indicated by the connectivity graph. In more detail, the particles are modeled using Gaussian random fields on the sphere. Thus, the shape of the particles can be characterized by a mean radius $\mu$ and the angular power spectrum $A:[0,\infty)\rightarrow[0,\infty)$, see [28]. The angular power spectrum is approximated by the function $A(l)=\frac{al+b}{l^{2}+cl+d}$ with coefficients $a=0.4241$, $b=0.356$, $c=-3.858$ and $d=3.903$. In more detail, we do not use the mean radius $\mu$ directly but we generate the particles in a way that their volume is proportional to the volume of corresponding Laguerre cells. The particles are sampled with the boundary conditions indicated by the connectivity graph using a special sampling algorithm that creates only realizations of the given Gaussian random field that fulfill those conditions.

The schematic depiction in Fig. 2 shows the particles with the tessellation and the connectivity graph. One can clearly see that the particles touch each other where indicated by the graph and on the other hand also fill their respective Laguerre polytopes. Fig. 2 shows the system of connected particles without the tessellation and the connectivity graph as these are auxiliary tools that are no longer needed after the creation of particles.

Finally, a morphological smoothing [8] is performed on the system of connected particles to mimic the effect of binder. In the given sample the volume fraction of the binder as well as the contrast in the tomographic images were too low to identify and model the binder as separate phase. From the known production process (slurry coating) we assume that this approach produces a similar effect as depicted in Fig. 2.

In Fig. 3 and 3, a cutout from the tomographic image data can be compared to a simulated anode structure. A very good visual accordance can be observed.

The model described so far is an excellent tool to generate virtual anode microstructures of energy cells, which are characterized by a high volume fraction of the solid phase. However, note that it can not directly be used to model the morphology of anodes in power cells, because due to the lower volume fraction of the solid phase, the boundary conditions for particles cannot be fulfilled reasonably. Therefore, an extension of the model has been proposed [66]. To account for the lower volume fraction, a Laguerre tessellation with marked polytopes is used. The polytopes are marked either as ‘filled’, i.e., a particle is placed in the polytope, or as ‘empty’, which means that no particle is placed here. Thereby, a reasonable allocation of the different polytopes as well as full connectivity of the resulting structure is ensured. Furthermore, the model is able to include anisotropy effects of the solid phase, i.e., particles can be elongated in horizontal direction rather than in vertical direction. This results in a remarkable flexibility such that the model can be used to create a broad spectrum of virtual anode microstructures with a variety of morphological properties, see Fig. 3 for some examples.

In this study, we focus on electrochemical simulations of anode microstructures in energy cells. Thus, all microstructures which are discussed in the present paper, are created using the energy cell model [13]. Electrochemical simulations on virtual structures generated by the power cell model are subject of further research.

2.2 Microscopic cell modeling

In this section we recollect the equations that describe the transport of lithium ions in a three-dimensional microstructure generated by the method as described in the previous section. The physical model has been derived based on species, charge and energy conservation to yield a set of equations that describe the spatial and temporal distribution of lithium ions, electrical potentials and temperature [31, 32]. However, in order to fit to the isothermal plating model we neglect the effects of heat production here and consider an isothermal system where temperature enters as a constant model parameter.

For the purpose of this study we restrict ourselves to half-cell simulations, i.e. we consider a porous graphite electrode modeled by the discussed stochastic method against a lithium foil as counter-electrode. The simulation domains are connected with the external operation conditions through dedicated current collector phases on each electrode. The remaining space of the computational domain is filled with ion conductive electrolyte.

Within the graphite particles we have the following equations for lithium concentration $c_{\mathrm{Gr}}$ and electrical potential $\Phi_{\mathrm{Gr}}$

[TABLE]

where $D_{\mathrm{Gr}}$ is the lithium diffusion coefficient and $\sigma_{\mathrm{Gr}}$ is the electrical conductivity of the material. The ion flux and electric current density are denoted by $\bm{N_{\mathrm{Gr}}}$ and $\bm{j_{\mathrm{Gr}}}$, respectively. Also in the domains of the lithium foil and the current collectors electronic conduction is considered and hence (2.2) is also solved in these domains (with the respective conductivities of course). Since there is no intercalation and diffusion of ions neither in the lithium counter-electrode nor in the current collectors, (2.1) is only relevant for the graphite domain.

Within the electrolyte domain, lithium concentration $c_{\mathrm{El}}$ and electrochemical potential $\varphi_{\mathrm{El}}$ are coupled through

[TABLE]

where $t_{+}$ is the transference number of lithium in the electrolyte, $\kappa_{\mathrm{El}}$ is the ionic conductivity of lithium inside the electrolyte and $F$ is the Faraday constant. The derivative of the electrolyte chemical potential is given as

[TABLE]

with $T$ denoting the temperature, $R$ the gas constant and $f_{+}$ the activity coefficient.

On the interfaces between the electrodes and electrolyte two types of reactions need to be considered: That is an intercalation reaction on the graphite side and a lithium deposition reaction on the counter-electrode side. The different phases (graphite, electrolyte and counter-electrode) are coupled via interface conditions

[TABLE]

with “So” (solid) being either graphite or metallic lithium. By convention the interface normal $\bm{n}_{\mathrm{So-El}}$ points from solid into the electrolyte. These conditions express the continuity of the current and mass fluxes between the phases. The current flow through these interfaces depends on the corresponding reactions. The intercalation reaction is described by a Butler-Volmer-like expression [32]

[TABLE]

where the overpotential is given by $\eta_{\mathrm{Gr-El}}=\Phi_{\mathrm{gr}}-U_{\mathrm{0}}^{\mathrm{Gr}}-\varphi_{\mathrm{El}}$. The electrode’s open-circuit potential $U_{\mathrm{0}}^{\mathrm{Gr}}$ is a concentration dependent material property. The rate constant $i_{\mathrm{Gr-El}}^{\mathrm{00}}$ depends on the lithium salt and the electrolyte composition. The transfer coefficients $\alpha_{c}+\alpha_{a}=1$ of the intercalation reaction were assumed to be symmetrical ($\alpha_{a,c}=0.5$). The form of (2.7) differs from the usual Butler-Volmer expression by omitting the common $\left(c_{\mathrm{Gr}}^{\mathrm{max}}-c_{\mathrm{Gr}}\right)^{\alpha}_{a}$ term, since a rigorous thermodynamically consistent derivation does not in general yield this prefactor [32]. The Butler-Volmer model and the used exchange current include a relation between the potentials and the current flux. The exponential shape of the exchange current introduces highly pronounced non-linearities into the numerical system. The reaction at the counter-electrode is described by a simple exchange current, which minimizes the effect of the counter-electrode on the simulation results

[TABLE]

where the rate constant is given by $i_{\mathrm{CE-El}}^{\mathrm{00}}$. On the lithium electrode the overpotential is simply given by $\eta_{\mathrm{CE-El}}=\Phi_{\mathrm{CE}}-\varphi_{\mathrm{El}}$. The remaining interface conditions are shown in Tab. 1.

These equations describe an ideal battery, i.e. no degradation processes are considered. This extension is outlined in the next section.

2.3 Electrochemical degradation modeling

The focus in this work is the degradation process lithium plating, where the lithium ions form an unwanted metallic phase on the surface of the intercalation material of the negative electrode. The electrochemical modeling of this process is briefly described in the following section, details can be found in the corresponding publication [22]. Two states of the lithium ions are of relevance: the lithium ions dissolved in the electrolyte $\mathrm{Li}_{\mathrm{El}}^{\mathrm{+}}$ and the metallic/plated lithium phases $\mathrm{Li}_{\mathrm{Pl}}^{\mathrm{0}}$. The transition between these two is expressed by the reaction

[TABLE]

The overpotential of the plating and stripping reaction is defined by the difference between the electrochemical potential $\tilde{\mu}$ of the two lithium phases involved[44]

[TABLE]

With the definition of the reference state $\mu_{\mathrm{Li^{0}}}^{\mathrm{Pl}}=\mu_{\mathrm{Li^{+}}}^{\mathrm{Pl}}+\mu_{\mathrm{e^{-}}}^{\mathrm{Pl}}$ and the electrochemical potentials of lithium ions inside the electrolyte and a solid phase (see [22]), the overpotential (2.10) can be rewritten to

[TABLE]

where $\varphi_{\mathrm{i}}^{\mathrm{p}}$ denotes the electrochemical potential of species $i$ in phase $p$ with respect to the reference state $\mu_{\mathrm{Li^{0}}}^{\mathrm{Pl}}$. Plating of lithium is occurring if the overpotential reaches negative values ($\eta_{\mathrm{Pl/St}}<0$). The metallic lithium phase on the surface of the anode is not in a stable configuration, even if no external current is applied to the system. As soon as lithium is plated on the surface of the active material, the lithium metal can react with its surroundings in different ways. The reaction between the plated lithium and the electrolyte results in the growth of a solid-electrolyte interphase (SEI), which leads to an irreversible loss of lithium [63]. Apart from phenomenological models no theory exists which combines lithium intercalation, lithium plating and SEI growth. Hence, this irreversible pathway is not included in the present paper. The plated lithium can also intercalate charge-neutrally into the supporting graphite. This reaction represents a reversible lithium stripping pathway. We are not aware of any literature regarding the identification and parameterization of the charge-neutral reintercalation. Hence the direct reintercalation from the plated lithium into graphite was neglected in this work.

The stripping and plating reaction of the lithium is described by a Butler-Volmer-like equation

[TABLE]

The Bulter-Volmer-like expression is derived for non-vanishing phases. But, the plated lithium phase can completely desolve during stripping. Hence, the vanishing of the plated lithium phase is considered in the exchange current by the numerical regularization function $f_{\mathrm{pre}}\left(n_{\mathrm{Li}}\right)$, which depends on the amount of plated lithium $n_{\mathrm{Li}}$

[TABLE]

Based on numerical considerations [22], we set the constant $n_{\mathrm{Li}}^{\mathrm{const}}$ to a value corresponding to a thickness of plated lithium of $0.48\,\text{nm}$. For partially covered surfaces more detailed models are necessary to capture the stripping of partially covered surfaces including the surface-driven dissolution of small lithium droplets.

At the interface between the plated lithium and the electrolyte the current through the interface is equal to the stripping current

[TABLE]

All the interface conditions which are relevant for electrochemical simulations in this paper are listed in Tab. 1.

In this paper the stripping process of plated lithium is simulated by including the plated lithium into the 3D microstructure as an additional volume phase. In Fig. 4 an example of a 3D microstructure with plated lithium is shown. The porous electrode (red/right) is generated by the stochastic generation algorithm as described in Sec. 2.1. Additionally, regions with plated lithium are positioned randomly at the separator-electrode interface.

The microstructure shown in Fig. 4 is used for the model order reduction experiment described in Sec. 3.2.

One important deviation from the microstructure model [13, 66] was to introduce a third phase in the anode structure, plated lithium. We use a fairly simple model to create a slightly plated structure as initial condition. These initial conditions are, like the microstructure model, simulated stochastically. This means that for every run the microstructure model extended by the plated lithium phase generates a new structure but with similar statistical properties.

The method used to create the plated lithium phase is a germ-grain model. This means that in a first step germs are simulated and in the second step grains are placed around the germs [7]. The parameters of this model are the intensity $\lambda$ of the Poisson process that is used to generate the germs and the grain radius $r_{s}$. Due to the lack of experimental data the values where chosen to $\lambda=0.01/$\mathrm{\SIUnitSymbolMicro m}$^{2}$ and $r_{s}=2.2\,$\mathrm{\SIUnitSymbolMicro m}$$. The general idea of the model is to place plating germs randomly on the surface of the particles in the electrode and then initialize plating around the germs on all the points on the structure’s surface within the radius given by $r_{s}$. In more detail the following is done, see Fig. 4:

• •

Select points $\{p_{i},i\in\mathbb{N}\}$ at the separator-anode interface via a Poisson point process with intensity $\lambda$.

• •

For each point $p_{i},i\in\mathbb{N}$ do the following:

• –

Find the first point $p_{i}^{*}$ where the straight line from $p_{i}$ towards the anode current collector interface meets the anode structure. Let $\Theta$ be the particle on the boundary of which $p_{i}^{*}$ is placed.

• –

Consider a sphere $B(p_{i}^{*},r_{s})$ around $p_{i}^{*}$ with radius $r_{s}$. Let $L_{i}=\{x\in\mathbb{R}^{3}:x\in B(p_{i}^{*},r_{s})\text{ and }x\in\partial\Theta\}$, where $\partial\Theta$ is the boundary of $\Theta$.

• •

Let $L=\bigcup\limits_{i\in\mathbb{N}}L_{i}$ be the plated lithium phase.

The plated lithium phase $L$ is then discretized as a one voxel thick phase on the surface of the particles. The material parameters and reaction constants of the ion transport and plating model are adapted from a previous publication [22].

2.4 Discretization and high-dimensional simulation

For the spatial discretization of the presented plating model, a cell-centered finite volume scheme on a uniform voxel grid is considered. Hence the discretization is naturally conservative. Furthermore, with the simple grid structure meshing of the complex microstructure is straight-forward. The Butler-Volmer interface conditions (2.6) are prescribed as the numerical flux across the respective domain interfaces, leading to a global space differential operator on the entire computational domain. Choosing implicit Euler time stepping for time discretization, we obtain a series of discrete nonlinear equation systems of the form

[TABLE]

Here, $V_{h}$ denotes the discrete finite volume space of locally constant grid functions, $c_{\mu}^{(n)},\varphi_{\mu}^{(n)}\in V_{h}$ denote the concentration and potential fields at time step $n$ for some $p$-tuple of model parameters $\mu$ contained in a parameter domain of interest $\mathcal{P}\subset\mathbb{R}^{p}$, and $A_{\mu}:V_{h}\oplus V_{h}\to V_{h}\oplus V_{h}$ is the finite volume space differential operator. The system is closed by $c^{(0)}_{\mu}=c_{0}$ for some fixed initial lithium distribution $c_{0}\in V_{h}$. The time step size $\Delta t^{(n)}$ is chosen adaptively for each time step to accommodate the different time scales during and after the stripping of the plated lithium. The nonlinear equation systems are solved in BEST using Newton’s method and an algebraic multigrid solver for the solution of the linear correction equations. Details on the discretization are provided in a previous publication [50].

2.5 Model order reduction and reduced simulation

The computation of a single solution trajectory $c_{\mu}^{(n)},\varphi_{\mu}^{(n)}$ requires many hours, even for relatively small geometries (cf. Sec. 3.2). In order to make parameter studies computationally feasible, reduced basis model reduction techniques[53, 23, 19] are applied which have been implemented in our model order reduction library pyMOR [39, 52]. This allows us to obtain a quickly solvable reduced order model approximating the full order model (2.4).

To construct the reduced order model, solutions of (2.4) are computed for few appropriately selected parameters $\mu_{1},\ldots,\mu_{S}$. Various advanced algorithms exist for the selection of these snapshot parameters, often based on a greedy search procedure (cf. [53, 23, 19]). In our basic test case with a one-dimensional parameter domain $\mathcal{P}$ (Sec. 3.2), a simple equidistant parameter sampling will be sufficient, however.

From this data, reduced approximation spaces $\tilde{V}_{c},\tilde{V}_{\varphi}$ for the concentration and potential fields are constructed via proper orthogonal decomposition (POD, principal component analysis) [56] of the snapshot data sets $\mathcal{S}_{c}=\{c^{(n)}_{\mu_{s}},c^{(n,i)}_{\mu_{s}}\}$, $\mathcal{S}_{\varphi}=\{\varphi^{(n)}_{\mu_{s}},\varphi^{(n,i_{k})}_{\mu_{s}}\}$. Here, $c^{(n,i)}_{\mu_{s}}$, $\varphi^{(n,i)}_{\mu_{s}}$ denote the intermediate Newton stages during the solution of (2.5), which are included for improved numerical stability. By construction, we in particular have $\tilde{V}_{c}\subseteq\operatorname{span}\mathcal{S}_{c}$ and $\tilde{V}_{\varphi}\subseteq\operatorname{span}\mathcal{S}_{\varphi}$. While $\dim V_{h}$ is in the order of $10^{6}$, we typically have $\dim\tilde{V}_{c},\dim\tilde{V}_{\varphi}<100$, which makes significant computational speedups possible.

After the reduced approximation space $\tilde{V}=\tilde{V}_{c}\oplus\tilde{V}_{\varphi}$ has been computed, the reduced order model is obtained via Galerkin projection of (2.4) onto $\tilde{V}$. I.e., we solve

[TABLE]

$\tilde{c}^{(0)}_{\mu}=P_{\tilde{V}_{c}}(c_{0})$, where $P_{\tilde{V}}$ / $P_{\tilde{V}_{c}}$ denotes the $L^{2}$-orthogonal projection onto $\tilde{V}$ / $\tilde{V}_{c}$.

However, even though (2.5) contains only $\dim\tilde{V}$ degrees of freedom, its solution requires the evaluation of the high-dimensional system operator $A_{\mu}$. This strongly limits the achievable speedup in computation time when solving (2.5) instead of (2.4).

To overcome this issue, $A_{\mu}$ is replaced by a quickly evaluable low-order approximation using the empirical interpolation technique [5, 11]: for an arbitrary (nonlinear) operator $O:X\to Y$, the EI-Greedy algorithm is used to compute a low-order interpolation space $\tilde{Y}\subseteq Y$ from evaluations of $T$ on given solution trajectories, after which the interpolated operator $\mathcal{I}_{M}[O]$ is determined by requiring it to agree with $O$ at appropriate $M=\dim\tilde{Y}$ interpolation degrees of freedom $\pi_{1},\ldots,\pi_{M}:Y\to\mathbb{R}$. I.e., for all $x\in X$ we have

[TABLE]

Due to the locality of finite volume operators, the point evaluations $\pi_{m}(O(x))$ can be computed quickly and independently from the dimension of $V_{h}$.

Since the potential part of $A_{\mu}$ vanishes identically for solutions of (2.4), a direct application of empirical interpolation to $O=A_{\mu}$ results in an unusable approximation, however. Instead, we further decompose $A_{\mu}$ and only use empirical interpolation for appropriate sub-operators.

In the following, we are interested in the behavior of the model in dependence on the applied current density. In this case, with $\mu$ being the applied current density, $A_{\mu}$ decomposes as

[TABLE]

where $A^{(bv)}$, $A^{(1\!/\!c)}$ are the parameter-independent nonlinear parts of $A_{\mu}$ corresponding to the Butler-Volmer interface terms and the summand in (2.4) containing $\partial\mu/\partial c_{El}$. Assuming constant $t_{+}$, the remainder $A^{(aff)}_{\mu}$ is affine linear and decomposes as

[TABLE]

where $A^{(const)}$ is constant and $A^{(bdn)}$, $A^{(lin)}$ are linear, non-parametric operators.

Now we apply the EI-Greedy algorithm on the training datasets $\mathcal{S}_{*}=\{A^{(*)}(c^{(n)}_{\mu_{s}},\varphi^{(n)}_{\mu_{s}}),\ A^{(*)}(c^{(n,i)}_{\mu_{s}},\varphi^{(n,i)}_{\mu_{s}})\}$, $*\in\{bv,1\!/\!c\}$, to obtain empirically interpolated operators $\mathcal{I}_{M^{(*)}}[A^{(*)}]\approx A^{(*)}$, which give us the approximation

[TABLE]

Substituting (2.20) into (2.5) we arrive at the fully reduced model

[TABLE]

with $\tilde{c}^{(0)}_{\mu}=P_{\tilde{V}_{c}}$. After pre-computation of the matrix representations of the linear (constant) operators $P_{\tilde{V}}\circ A^{(const)},P_{\tilde{V}}\circ A^{(bnd)},P_{\tilde{V}}\circ A^{(lin)}:\tilde{V}\to\tilde{V}$, as well as the projections from the interpolation spaces for $A^{(bv)}$, $A^{(1\!/\!c)}$ onto $\tilde{V}$, the solution of (2.5) can be obtained quickly for arbitrary new parameters $\mu$ with an effort that only depends on $\dim\tilde{V}$, $M^{(bv)}$ and $M^{(1\!/\!c)}$.

In the following experiments (see Sec. 3.2) we are interested in the cell potential as well as the average lithium concentration in the electrode as functions of time and the applied delithiation current density $\mu$. These quantities are linear functionals $s_{cp},s_{ac}:V_{h}\oplus V_{h}\to\mathbb{R}$, assigning to a state of the cell the respective quantity of interest. Due to their linearity, the vector representation for the evaluation of $s_{cp}$, $s_{ac}$ on $\tilde{V}$ can again be pre-computed, such that for any given solution of (2.5), $s_{cp}(\tilde{c}_{\mu}^{(n)},\tilde{\varphi}_{\mu}^{(n)})$, $s_{ac}(\tilde{c}_{\mu}^{(n)},\tilde{\varphi}_{\mu}^{(n)})$ are quickly obtained with an effort only depending on $\dim\tilde{V}$.

The model order reduction introduces an additional approximation error between full order and reduced order model that needs to be accounted for in the simulation workflow. As we are not aware of any rigorous error estimates which would provide sufficiently tight error bounds for the model under consideration, we here consider the following heuristic a posteriori error estimator [20]: In addition to $\tilde{V}$ we construct a second, larger validation space $\hat{V}=\hat{V}_{c}\oplus\hat{V}_{\varphi}\supset\tilde{V}_{c}\oplus\tilde{V}_{\varphi}=\tilde{V}$ and extended interpolation bases of dimensions $\hat{M}^{(bv)}>M^{(bv)}$ and $\hat{M}^{(1/c)}>M^{(1/c)}$, yielding a larger reduced model with solutions $(\hat{c}_{\mu}^{(n)},\hat{\varphi}_{\mu}^{(n)})\in\hat{V}$. Under the heuristical assumption that

[TABLE]

for all timesteps $n$, $\mu\in\mathcal{P}$ with a fixed $\Theta\in[0,1)$, from the triangle inequality we immediately obtain the error estimates

[TABLE]

The right-hand sides of (2.23) can be quickly computed at the expense of an additional solution of a second (slightly larger) reduced order model. Note that (2.22) is precisely the saturation assumption in the context of hierarchical a posteriori estimates for finite element schemes (see e.g. [4]).

2.6 Algorithmical integration and software interfaces

MULTIBAT aims to allow computationally fast studies of local effects in the complex microstructure of battery anodes within one software workflow. This is achieved by breaking the multi-disciplinary goal into task units and interfacing these units with BEST to varying degrees of depth.

The presented workflow resulting from these interfaces is schematically depicted in Fig. 5. It allows speeding up the numerical solution of the microscopic cell degradation modeling from Sec. 2.2 and 2.3 with discretization from Sec. 2.4 by BEST, using the randomly generated structures from Sec. 2.1 through the POD/EI based model order reduction approach from Sec. 2.5. The workflow has been used to create the results and speedups depicted in Sec. 3.

We introduce three distinct interfaces. The first is file-based and allows usage of randomly generated structures of Sec. 2.1 in BEST through a conversion tool. The conversion tool provides standard BEST geometry input which is matched with physical modeling input parameters and numerical solution parameters suitable for the models of Sec. 2.2 with the extensions from Sec. 2.3.

The second interface allows to extend the BEST numerical solution code to advanced interface flux modeling between the plated anode-lithium and the electrolyte from Sec. 2.3. The model extensions are compiled into the BEST library.

The third and most extensive interface is library-based and gives pyMOR runtime access to the solution process, vectors, discretization matrices, Jacobians, linear algebra solver and parameters of BEST through the BEST library to carry out POD/EI based MOR. The separation is strict: All MOR-related operations and the Newton methods are carried out in pyMOR and all evaluations of nonlinear operators and Jacobians are carried out by the BEST library ordered by pyMOR.

3 Workflow demonstration on lithium stripping case study

3.1 Microstructure generation and electrochemical verification

The stochastic microstructure model used in this work was parameterized on real tomographic image data [13]. As mentioned in Sec. 2.1, the validity of the structural parameterization was investigated through spatially resolved electrochemical simulations [21]. The validated stochastic microstructure model is used in this work. In the following a short summary of the electrochemical validation is given. 20 simulated realizations of the stochastic microstructure model and 20 microstructure cutouts from the tomographic image data are used as electrode structure samples for electrochemical simulations. These microstructures are delithiated with a constant current. The simulation results were compared using various electrochemical quantities, such as local current density and lithium concentration. A very good agreement between the real and virtual microstructures was found, see Fig. 6. The advantage of spatially resolved electrochemical simulations is the access to localized inhomogeneities. The spatial distribution of the electrolyte concentration for two cutouts of real and virtual microstructures is shown as an example in Fig. 6.

Both cutouts exhibit similar features: less-than-average connected pores and large particles. Apart from the visual similarity between the real and virtual cutouts, averaged quantities were used for a more quantitative comparison. The average lithium concentration in the electrolyte in through direction (from one current collector towards the other) is shown in Fig. 6.

The general shapes of the concentration functions are nearly identical. The superposition of transport within the electrolyte and deintercalation of lithium from the solid phase results in a nonlinear gradient. Without any sources of lithium a linear concentration gradient forms in the separator. More details regarding the electrochemical validation can be found in the corresponding publication [21].

3.2 Model order reduction

As a first numerical test for the entire developed modeling and simulation workflow (see Sec. 2), we simulated the full model (including plated lithium) on a randomly generated half-cell geometry of size $44$\mathrm{\SIUnitSymbolMicro m}$\times 44$\mathrm{\SIUnitSymbolMicro m}$\times 65.6$\mathrm{\SIUnitSymbolMicro m}, which is meshed with a grid of $100\times 100\times 149$ voxels (see Fig. [4](#S2.F4)). This size is required to cover a representative volume containing several particles in each direction, and, at the same time, to achieve a sufficient resolution to resolve a relevant part of the electrode’s morphology. Starting with plated lithium and a high lithium concentration in the electrode we performed a delithiation simulation and hence expect to see lithium stripping. We simulated 60 seconds with constant current densities in the interval $\mathcal{P}=[2.5,250]$\mathrm{A}\mathrm{/}\mathrm{m}^{2}, which corresponds to currents from $4.84\text{\,}\mathrm{nA}$ to $484\text{\,}\mathrm{nA}$ or to C-rates from $\mathrm{C}/10$ to $10\mathrm{C}$.

A single simulation of the full order model (2.4) requires around 16 hours (cf. Tab. 2). To generate the snapshot data for the computation of the reduced order model (2.5), the full order model (2.4) was solved for the three delithiation current densities $\min\mathcal{P}$, $\max\mathcal{P}$ and $(\min\mathcal{P}+\max\mathcal{P})/2$. The reduced spaces $\tilde{V}_{c}$, $\tilde{V}_{\varphi}$, $\hat{V}_{c}$, $\hat{V}_{\varphi}$, as well as the interpolation spaces for $\mathcal{I}_{M^{(bv)}}[A^{(bv)}]$, $\mathcal{I}_{M^{(1\!/\!c)}}[A^{(1\!/\!c)}]$ were computed using the POD and EI-Greedy algorithms. To ensure that a numerically stable reduced order model is obtained, a small relative error tolerance of $10^{-7}$ was chosen, using in each case 97% of the resulting basis vectors for the construction of the reduced order model and all basis vectors for construction of the validation model used for the error estimator (2.23). The resulting spaces are of the following dimensions: $\dim\tilde{V}_{c}=178$, $\dim\hat{V}_{c}=183$, $\dim\tilde{V}_{\varphi}=67$, $\dim\hat{V}_{\varphi}=69$, $M^{(bv)}=924$, $\hat{M}^{(bv)}=952$, $M^{(1\!/\!c)}=997$ and $\hat{M}^{(1\!/\!c)}=1027$.

To validate the resulting reduced order model (2.5), we compared the solutions of (2.5) to the full order model (2.4) for 10 random parameter values $\mu_{i}\in\mathcal{P}\,,i=1,\ldots,10$ in addition to the three snapshot parameters used for training. While achieving a relative model reduction error of at most $4.81\cdot 10^{-4}$  resp. $4.50\cdot 10^{-3}$  for the concentration and potential variables (Fig. 7), the reduced order model can be simulated in less than 8 minutes, yielding a speedup factor of 120. Since the generation of the reduced model from the high-dimensional snapshot data is faster than a single solution of the full order model, an overall saving of computation time is already achieved for one additional model simulation (Tab. 2). For $\Theta=0$ the error estimator (2.23) overestimates the real model reduction error in these 13 parameters by a factor of at most 1.08 (3.46) for the concentration (potential) and underestimates the error by a factor of at most 2.89 (1.45). The estimator was evaluated for 100 additional current densities in $\mathcal{P}$, yielding a maximum estimated relative error of $3.61\cdot 10^{-4}$ resp. $5.55\cdot 10^{-3}$ for concentration and potential.

In Fig. 8 and 8, the cell potential and average lithium concentration in the electrode have been plotted over the transferred charge for the 10 random test parameters. A short interpretation of these results is given in the subsequent section. Overall, no visual distinction between the data generated by the reduced and full order models can be made.

3.3 The lithium stripping process

In Fig. 8 the cell voltage is shown for 10 of the applied currents. The cell voltages for all applied currents exhibit a similar shape. A voltage plateau at the start of delithiation is followed by a rise, which is in turn succeeded by a region following the shape of the open-circuit potential $U_{\mathrm{0}}^{\mathrm{Graphite}}$. The initial voltage plateau results from the stripping reaction (see Eq. (2.9)). The increase in cell voltage begins as soon as the majority of the plated lithium is consumed. The apparent plateau afterwards is the cell voltage of the supporting graphite at about $75\%$ state of charge. This equilibrium potential is shifted with an overpotential, which depends on the applied current. A transferred charge of $2\text{\,}\mathrm{nA}\text{\,}\mathrm{h}$ corresponds to a variation in the state of charge of $4$%$$ since the used microstructure has a maximum capacity of $49\text{\,}\mathrm{n}\mathrm{A}\mathrm{h}$. Large applied stripping currents lead to a fast decrease of lithium concentration at the surface of graphite. The solid diffusion can not equilibrate the lithium concentration in the electrode in the same rate for large currents as for small currents. Therefore, the state of charge at the surface will vary more than the overall change in state of charge. The open-circuit potential will then increase faster for larger currents. Thus resulting in a more sloped cell voltage for the larger currents. The length of the stripping plateau in the cell voltage depends on the applied current. For smaller currents, the change from constant potential to graphite dominated region is at lower transferred charge. The intercalation of lithium during the stripping of the plated lithium leads to an increase of the lithium concentration in the solid phase, as can be seen in Fig. 8. For low applied currents a net intercalation during the lithium stripping is visible. As soon as the majority of the plated lithium is dissolved a net delithiation exists. More information about the distribution of the delithiation current on the stripping and intercalation reaction are provided in a recent publication [22].

The cell voltages and average lithium concentrations obtained from the model reduction experiments are identical to the ones obtained from the full order model. This indicates that the reduced model sufficiently represents the electrochemical relevant regions in the simulation domain.

4 Conclusion & outlook

We conclude that it is absolutely possible to do a lot of in-depth research on lithium-ion cells virtually. In this work we have shown one approach to solve many of the existing problems using simulation techniques for the investigation in lithium-ion battery cells.

First, the limitation of 1D or pseudo 2D models which consider only averaged structural quantities and thus neglect all local effects can be overcome by switching to spatially resolved 3D models. Hence we presented a physics-based model that describes the cell’s behavior on a microscopic scale and includes effects of lithium plating and stripping. Based on the software tool BEST the mathematical model was implemented and solved in a three-dimensional geometry.

Another current limitation is that the acquisition of tomographic 3D images as basis for simulations is costly and restricts the ability to simulate new structures that have not been produced experimentally. This limitation is overcome by the usage of a 3D stochastic microstructure model which has been implemented in the software library GEOSTOCH. Once the model is fitted to a material by the usage of tomographic 3D images it is possible to generate arbitrary many virtual cutouts with arbitrary sizes. By reasonable changes on the model parameters it is even possible to generate structures that have not been processed experimentally in order to investigate their properties and to test their performance.

Finally, a great problem of simulation-based parameter studies in particular including the plating/stripping behavior is the extensive simulation runtime. This problem is solved using model order reduction methods implemented in the software library pyMOR, which speed up the simulation of similar cycles significantly. Tab. 2 shows a speedup of factor 120 for the simulation of the reduced model in comparison to a full order model simulation.

Overall we have shown that the combinations of all the methods described above work well in a demo scenario and can improve the accuracy of the geometry models, increase the computational speed considerably, and extended the predictive power of the electrochemical battery models.

Acknowledgement

This work was partially funded by BMBF under grant numbers 05M13VUA, 05M13PMA, 05M13AMF and 05M13CLA in the programme “Mathematik für Innovationen in Industrie und Dienstleistungen”.