Multiphase Porous Electrode Theory

TL;DR

This paper introduces a comprehensive multiphase porous electrode theory framework for batteries, incorporating phase separation, advanced reaction kinetics, and thermodynamics, implemented in an open-source software package called MPET.

Contribution

It extends traditional porous electrode models to include multiphase materials and complex reaction kinetics using electrochemical thermodynamics, enabling more accurate battery simulations.

Findings

Demonstrates the software's ability to model phase-separating materials.

Shows integration of Marcus-Hush-Chidsey kinetics for electron transfer.

Provides example simulations illustrating the model's capabilities.

Abstract

Porous electrode theory, pioneered by John Newman and collaborators, provides a useful macroscopic description of battery cycling behavior, rooted in microscopic physical models rather than empirical circuit approximations. The theory relies on a separation of length scales to describe transport in the electrode coupled to intercalation within small active material particles. Typically, the active materials are described as solid solution particles with transport and surface reactions driven by concentration fields, and the thermodynamics are incorporated through fitting of the open circuit potential. This approach has fundamental limitations, however, and does not apply to phase-separating materials, for which the voltage is an emergent property of inhomogeneous concentration profiles, even in equilibrium. Here, we present a general theoretical framework for "multiphase porous electrode theory" implemented in an open-source software package called "MPET", based on electrochemical nonequilibrium thermodynamics. Cahn-Hilliard-type phase field models are used to describe the solid active materials with suitably generalized models of interfacial reaction kinetics. Classical concentrated solution theory is implemented for the electrolyte phase, and Newman's porous electrode theory is recovered in the limit of solid-solution active materials with Butler-Volmer kinetics. More general, quantum-mechanical models of Faradaic reactions are also included, such as Marcus-Hush-Chidsey kinetics for electron transfer at metal electrodes, extended for concentrated solutions. The full equations and numerical algorithms are described, and a variety of example calculations are presented to illustrate the novel features of the software compared to existing battery models.

Figures

Equations

$$\displaystyle\frac{\partial\mathopen{}\mathclose{{}\left(\epsilon c_{\ell,i}}\right)}{\partial t}=-\boldsymbol{\nabla}\cdot\mathbf{F}_{\ell,i}+R_{V,i}$$

$$\displaystyle\frac{\partial\mathopen{}\mathclose{{}\left(\epsilon\rho_{e}}\right)}{\partial t}\approx 0=-\boldsymbol{\nabla}\cdot\mathbf{i}_{\ell}+\sum_{i}z_{i}eR_{V,i},$$

$${\mathbf{i}}_{\ell }$$

$$z_n{c}_{\ell ,n}=-\sum _{i\ne n}{z}_i{c}_{\ell ,i},$$

$$\displaystyle\frac{\partial\mathopen{}\mathclose{{}\left(\epsilon c_{\ell}}\right)}{\partial t}=\frac{1}{\nu_{i}}\Big{(}-\boldsymbol{\nabla}\cdot\mathbf{F}_{\ell,i}+R_{V,i}\Big{)}.$$

$$\displaystyle 0=-\boldsymbol{\nabla}\cdot\mathbf{i}_{\ell}+\mathopen{}\mathclose{{}\left(z_{+}eR_{V,+}+z_{-}eR_{V,-}}\right).$$

$$\displaystyle\mu_{\ell,i}=k_{\mathrm{B}}T\ln\mathopen{}\mathclose{{}\left(a_{\ell,i}}\right)+\mu_{\ell,i}^{\Theta}+z_{i}e\phi_{\ell}=k_{\mathrm{B}}T\ln\mathopen{}\mathclose{{}\left(\widetilde{c}_{\ell,i}}\right)+\mu_{\ell,i}^{\mathrm{ex}}$$

$$\displaystyle e\phi_{\ell}^{r}=\mu_{\ell,\mathrm{Li}^{+}}+\mathrm{constant}=k_{\mathrm{B}}T\ln\mathopen{}\mathclose{{}\left(a_{\ell,\mathrm{Li}^{+}}}\right)+e\phi_{\ell}+\mu_{\ell,\mathrm{Li}^{+}}^{\Theta}+\mathrm{constant}.$$

$$\displaystyle z_{n}e\phi_{\ell}^{q}=\mu_{\ell,n}^{\mathrm{ex}}=z_{n}e\phi_{\ell}+k_{\mathrm{B}}T\ln\mathopen{}\mathclose{{}\left(\gamma_{\ell,n}}\right)+\mu_{\ell,n}^{\Theta}$$

$${\mathbf{F}}_{\ell ,i}=-M_{\ell ,i}{c}_{\ell ,i}\nabla {\mu }_{\ell ,i}$$

$${\mathbf{F}}_{\ell ,i}$$

$$\displaystyle\mathbf{F}_{\ell,i}=-\mathopen{}\mathclose{{}\left(D_{\ell,\mathrm{chem},i}\boldsymbol{\nabla}c_{\ell,i}-\frac{D_{\ell,i}}{\widetilde{T}}z_{i}c_{\ell,i}\boldsymbol{\nabla}\widetilde{\phi}_{\ell}}\right).$$

$${\mathbf{F}}_{\ell ,i}$$

$${\mathbf{F}}_{\ell ,i}=\sum _j{U}_{ij}\nabla {\mu }_{\ell ,j}.$$

$$\displaystyle c_{\ell,i}\boldsymbol{\nabla}\mu_{\ell,i}=\sum_{j}K_{ij}\mathopen{}\mathclose{{}\left(\mathbf{v}_{j}-\mathbf{v}_{i}}\right)=k_{\mathrm{B}}T\sum_{j}\frac{c_{\ell,i}c_{\ell,j}}{c_{T}\mathscr{D}_{ij}}\mathopen{}\mathclose{{}\left(\mathbf{v}_{j}-\mathbf{v}_{i}}\right)$$

$${\mathbf{F}}_{\ell ,+}$$

$${\mathbf{F}}_{\ell ,-}$$

$$D_{\ell }$$

$$\mathcal{D}$$

$$t_{+}^0$$

$$\displaystyle\mathbf{i}_{\ell}=-\frac{\sigma_{\ell}\epsilon}{\widetilde{T}\tau}\mathopen{}\mathclose{{}\left(\boldsymbol{\nabla}\phi_{\ell}^{r}+\frac{\nu k_{\mathrm{B}}T}{e}\mathopen{}\mathclose{{}\left(\frac{s_{+}}{n\nu_{+}}+\frac{t_{+}^{0}}{z_{+}\nu_{+}}-\frac{s_{0}c_{\ell}}{nc_{\ell,0}}}\right)\mathopen{}\mathclose{{}\left(1+\frac{\partial\ln\gamma_{\ell,\pm}}{\partial\ln\widetilde{c}_{\ell}}}\right)\boldsymbol{\nabla}\ln\widetilde{c}_{\ell}}\right)$$

$$\frac{1}{{\sigma }_{\ell }}$$

$$s_{-}{M}_{-}^{z_{-}}+s_{+}{M}_{+}^{z_{+}}+s_0{M}_0\rightleftharpoons n{\mathrm{e}}^{-}.$$

$$0=-\nabla \cdot {\mathbf{i}}_s-\sum _i{z}_{i}e{R}_{V,i}$$

$$\displaystyle 0=-\boldsymbol{\nabla}\cdot\mathbf{i}_{s}-\mathopen{}\mathclose{{}\left(z_{+}eR_{V,+}+z_{-}eR_{V,-}}\right)$$

$$\displaystyle\mathbf{i}_{s}=-\frac{\mathopen{}\mathclose{{}\left(1-\epsilon}\right)}{\tau}\sigma_{s}\boldsymbol{\nabla}\phi_{s}$$

$$\displaystyle G_{j,k}\mathopen{}\mathclose{{}\left(\phi_{j,k}-\phi_{j,k+1}}\right)=I_{j,k}$$

$$I_{j,k}-I_{j,k+1}={\int }_{S_{k+1}}{j}_{j,k+1}\mathrm{d}A,$$

$$G={\int }_{V_p}g\mathrm{d}V+{\int }_{A_S}{\gamma }_S\mathrm{d}A$$

$$g=g_{\mathrm{h}}+g_{\mathrm{nh}}+\dots$$

$$g_{\mathrm{nh}}=\frac{1}{2}\frac{1}{c_{s,\mathrm{ref}}^2}\nabla c_i\cdot \mathbf{\kappa }\nabla c_i$$

$${\mu }_i=\frac{\delta G}{\delta c_i}=\frac{\partial g}{\partial c_i}-\nabla \cdot \frac{\partial g}{\partial \nabla c_i}.$$

Full text

I Introduction

Lithium-based batteries have growing importance in global society Scrosati and Garche (2010) as a result of increased prevalence of portable electronic devices Tarascon and Armand (2001), and their enabling role in the transition toward renewable energy sources Whittingham (2014). For example, lithium batteries can help mitigate intermittency of renewable energy sources such as solar power, and lithium battery powered electric vehicles are facilitating movement away from liquid fossil fuels for transportation. Each of these growing areas demands high performance batteries, with requirements specific to the particular needs of the application driving specialized battery design for sub-markets. Thus, it is critical that battery models be based on the underlying physics, enabling them to greatly facilitate cell design to take best advantage of the existing battery technologies.

Lithium-ion batteries are generally constructed using two porous electrodes and a porous separator between them. The porous electrodes consist of various interpenetrating phases including electrolyte, active material, binder, and conductive additive. A schematic is shown in Figure 1. In a charged state, most of the lithium in the cell is contained in the active material within the negative electrode. During discharge, the lithium undergoes transport to the surface of the active material, electrochemical reaction to move from the active material to the electrolyte, transport through the electrolyte to the positive electrode, and reaction and transport to move into the active material of the positive electrode Winter and Brodd (2004); Newman and Thomas-Alyea (2004). Physical models must capture each of these behaviors accurately.

Complicating the situation further, the microstructure of the interpenetrating porous media within the electrodes can have a strong effect on the cell behavior Smith et al. (2009) as a result of inhomogeneities over a length scale smaller than the electrode but many times the size of a primary active particle García et al. (2005); Kehrwald et al. (2011); Novák et al. (2000); Ebner et al. (2014); Harris and Lu (2013). The active materials themselves also often have highly non-trivial behavior including poor connectivity with each other and the conductive matrix Stephenson et al. (2007); Gaberscek et al. (2007); Orvananos et al. (2014a); Ding et al. (2014) as well as complex material properties leading to deformations and accompanying stresses Aziz et al. (1991); Christensen and Newman (2006); Tang et al. (2009); Cheng and Verbrugge (2010); Bower et al. (2011); Di Leo et al. (2015); Woodford et al. (2010); Sethuraman et al. (2012); Qi and Harris (2010); Cogswell and Bazant (2012); Ulvestad et al. (2014); Cannarella and Arnold (2014) and phase separation Whittingham (2004); Reimers and Dahn (1992); Padhi et al. (1997); Chen et al. (2006); Ariyoshi et al. (2004); Ohzuku et al. (1995); Dahn (1991); Singer et al. (2014) during the intercalation process. Coupling these behaviors in multi-particle electrode environments leads to further complexities Christensen and Newman (2006); Dreyer et al. (2010); Dargaville and Farrell (2010); Li et al. (2014); Zhang et al. (2014); Ferguson and Bazant (2012); Dargaville and Farrell (2013); Bai and Tian (2013); Levi et al. (2013); Ferguson and Bazant (2014); Orvananos et al. (2015, 2014b).

Because capturing all the relevant physical processes from the atomistic length scale to the cell pack level within a single simulation is computationally intractable, various approaches have been developed to simulate aspects of battery behavior Grazioli et al. (2016); Miranda et al. (2015). In particular, porous electrode theory, pioneered by John Newman and co-workers over the past fifty years, has proven highly successful in describing the practical scale of individual cells Newman and Tiedemann (1975); Newman and Tobias (1962); Doyle et al. (1993); Fuller et al. (1994); Doyle et al. (1996); Darling and Newman (1997); Arora et al. (1999); Srinivasan and Newman (2004); Christensen and Newman (2006); De Vidts and White (1997); Verbrugge and Koch (2003); Zhang and White (2007); Dees et al. (2005, 2009); Baker and Verbrugge (2012); Bernardi and Go (2011); Gallagher et al. (2012); Lai and Ciucci (2011). This approach is based on volume averaging over a region of the porous electrode large enough to treat it as overlapping, homogeneous, continuous phases to describe the behavior of the electrons in the conductive matrix and the ions in the electrolyte. The behavior of the active material is treated by defining representative particles and placing them within the simulation domain as volumetric source/sink terms for ions and electrons according to the actual volume fraction of active material in the electrode. In this way, details on the length scale of transport within small active material particles can be consistently coupled to volume averaged transport over much larger length scales. Ref. Thomas et al. (2001) provides an excellent overview of the fundamentals of the theory.

As a result of volume averaging, heterogeneities over intermediate length scales are lost, causing inaccuracies in predictions. Efforts to capture these heterogeneities in simulations have had success in characterizing the consequences of the volume averaging procedure and providing more accurate alternatives, including refinements to microstructural parameters such as tortuosity in the volume averaged approach García et al. (2005); García and Chiang (2007); Bertei et al. (2013); Ebner et al. (2014); Latz and Zausch (2015). Nevertheless, simulations including complete microstructure information are much more computationally expensive than volume averaged approaches, so the simpler approach retains value in situations requiring faster model calculation or development.

Porous electrode theory has been developed and tested for decades for a variety of battery materials, but strictly speaking, it can only describe solid-solution active materials, whose thermodynamics are uniquely defined by fitting the open circuit voltage versus state of charge, or average composition. Active materials with more complex thermodynamics, resulting in multiple stable phases of different equilibrium concentrations, cannot be described, except by certain empirical modifications. Phase separating materials, such as lithium iron phosphate and graphite, can be accommodated by introducing artificial phase boundaries, such as shrinking cores Srinivasan and Newman (2004) or shrinking annuli Heß and Novák (2013), respectively, but this approach masks the true thermodynamic behavior.

Instead, the open circuit voltage of a battery is an emergent property of multiphase materials, which reflects phase separation in single particles Bai et al. (2011); Cogswell and Bazant (2012, 2013); de Klerk et al. (2016) and porous electrodes Lai and Ciucci (2011); Ferguson and Bazant (2012, 2014). It can only be predicted by modeling the free energy functional, rather than the voltage directly, and consistently defining electrochemical activities, overpotentials, and reaction rates using variational nonequilibrium thermodynamics Bazant (2013, 2017). This is the approach of “multiphase porous electrode theory” presented below. In principle, such models are required to predict multiphase battery performance over a wide range of temperatures and currents Ferguson and Bazant (2012, 2014); Smith et al. (2017a); Bazant (2017), as well as degradation related to mechanical stresses Cogswell and Bazant (2012, 2013) and side reactions that depend on the local surface concentration profile Thomas-Alyea et al. (2017).

Regardless of the thermodynamic model, volume averaged simulations using porous electrode theory are carried out in a number of ways. Newman’s dualfoil code uses a finite difference method via the BAND subroutine, and it is freely available Newman (1998) and commonly used. It has been developed and tested for decades, but it requires analytical derivative information about model equations to form the Jacobian, which makes modifications to the code less straightforward. The free energy approach we take here naturally describes both phase separating and solid solution materials using the same mathematical framework in a more user-friendly implementation. Popular commercial software packages such as COMSOL com have also been used to implement versions of porous electrode theory Stephenson et al. (2007); Bernardi and Go (2011); Dees et al. (2009); Cai and White (2011), usually using the finite element method. This has the advantage of being quick to set up, but it can involve costly software licenses. More importantly, the closed source means that detailed inspection of the software is impossible for the purpose of verifying, modifying, and improving the numerical methods for the particular problems investigated. A number of authors have also written custom versions of porous electrode software Ferguson and Bazant (2012); Dargaville and Farrell (2010); Torchio et al. (2016), for example using a manually implemented finite volume method and some general differential algebraic equation (DAE) solver for time advancement. This approach provides significant flexibility, but it is not common to share the code to facilitate use and inspection by a broader community, although Torchio et al. recently published their MATLAB implementation using the finite volume method with a variable time stepper Torchio et al. (2016). We take a similar approach in this work. More comprehensive reviews of commonly used simulation approaches can be found in refs. Botte et al. (2000); Ramadesigan et al. (2012).

Here, we present the equations and algorithms for a finite-volume based simulation software package, which implements multiphase porous electrode theory (“MPET”). The code is freely available Smith (2017), and it is developed in a modular way to facilitate modification and re-use. It is based only on open-source software and is written in Python, a modern, high-level language commonly used in the scientific computing community. Computationally expensive aspects of the code are all done using standard and freely available, open source numerical libraries. This takes advantage of Python’s ease of use for the model definitions while retaining the fast and vetted computation of libraries written in lower level languages like C and Fortran.

The paper is organized as follows. In Section II, we begin by presenting the full mathematical framework of MPET, based on the original formulation of Ferguson and Bazant Ferguson and Bazant (2012), with several modifications. First, we incorporate the standard description of transport in concentrated electrolytes based on Stefan-Maxwell coupled fluxes and chemical diffusivities Newman and Thomas-Alyea (2004); Valo̸en and Reimers (2005). Second, we capture the continued development by our group of phase field models for the active solid materials Bazant (2013); Singh et al. (2008); Burch and Bazant (2009); Bai et al. (2011); Cogswell and Bazant (2012, 2013); Ferguson and Bazant (2014); Smith et al. (2017a); de Klerk et al. (2016), which have increasingly been validated by direct experimental observations of phase separation dynamics Ferguson and Bazant (2014); Li et al. (2014); Lim et al. (2016); Bazant (2017); Guo et al. (2016); Thomas-Alyea et al. (2017); de Klerk et al. (2016). Third, we provide alternatives to the empirical Butler-Volmer model of Faradaic reaction kinetics Bard and Faulkner (2001); Newman and Thomas-Alyea (2004); Dreyer et al. (2016); Heubner et al. (2015), based on the quantum mechanical theory of electron transfer pioneered by Marcus Marcus (1956, 1993, 1965); Chidsey (1991) and extended here for concentrated solutions Smith et al. (2017b), motivated by recent battery experiments Bai and Bazant (2014). Fourth, the porous electrode model is modified to allow for different network connections between active particles, as well as half-cells with Li-foil counter electrodes and full two-electrode cells. In Section III, the equations are made dimensionless, and numerical methods to solve them are presented in Section IV, along with the overall software structure. In Section V, a variety of example simulations are presented to highlight the novel features of MPET compared to previous models, and the paper concludes with an outlook for future developments in Section VI.

II Model

As discussed above, the basic structure of the model involves volume averaging over a region larger than the particles of active material. Because the details of dynamics within the active material particles can strongly affect model predictions, they are simulated at their small length scale and treated as a source term for equations defined over the larger, electrode length scales. Thus, the model can be broken down into a number of scales. For the overall cell, we specify either the cell current density or voltage input profiles along with any series resistance. The unspecified current or voltage is an output of the simulation. At the electrode scale, we solve electrolyte transport equations, potential losses in the electron-conducting matrix, and potential drop between particles. At the particle scale, we simulate how they react with the electrolyte and the internal concentration dynamics. We will assume uniform (though not necessarily constant) temperature in the model derivation. The software currently only supports uniform and constant room temperature simulations, but we retain temperature factors here for generality. We will use the convention in referring to electrodes that the electrode which is negative/positive at open circuit and charged conditions is referred to as the anode/cathode.

II.1 Electrode Scale Equations

II.1.1 Electrolyte Model

The general form of the electrolyte model equations arises from statements of conservation of species and conservation of charge within the electrolyte phase of a quasi-neutral porous medium Newman and Tobias (1962); De Vidts and White (1997). We consider electrolytes of salts defined by $\sum_{i}\nu_{i}M_{i}^{z_{i}}$, with species $M_{i}$ having valence $z_{i}$, and $\nu_{i}$ the number of ions $M_{i}$ in solution from dissolving one molecule of the neutral salt. For example, for CuCl2, $\nu_{+}=1$, $M_{+}=\mathrm{Cu}$, $z_{+}=2$, $\nu_{-}=2$, $M_{-}=\mathrm{Cl}$, and $z_{-}=-1$. With electrolyte species flux, $\mathbf{F}_{\ell,i}$ defined per area of porous medium, conservation of species requires

[TABLE]

where $c_{\ell,i}$ is the concentration of species $i$ in the electrolyte, $\epsilon$ is the electrolyte volume fraction (porosity), and $R_{V,i}$ describes a volume-averaged reaction rate which is the result of interfacial electrochemical reactions with the active materials in which ions are added to or removed from the electrolyte.

In general, electrolyte transport models for porous media, the macroscopic charge density may be nonzero and varies in response to current imbalances. Diffuse electrolyte charge screens internal charged surfaces of porous membrane materials Dydek et al. (2011); Dydek and Bazant (2013); Schmuck and Bazant (2015); Yaroshchuk (2012) or conducting porous electrodes Biesheuvel et al. (2011); Mirzadeh et al. (2014) and provides additional pathways for ion transport by electromigration (“surface conduction”) and electro-osmotic flows. In addition to Faradaic reactions, porous electrodes can also undergo capacitive charging by purely electrostatic forces, as in electric double layer capacitors and capacitive deionization systems Biesheuvel et al. (2011), or hybrid pseudo-capacitors Biesheuvel et al. (2012). In batteries, however, such effects are usually neglected Newman and Tiedemann (1975); Thomas et al. (2002); Newman and Thomas-Alyea (2004), since the focus is on electrochemical, rather than electrostatic, energy storage, using highly concentrated electrolytes.

In these electrolyte models, we will assume they satisfy quasi-neutrality, i.e. there is no net charge in the electrolyte over the simulated length scales Dickinson et al. (2011). Charge conservation and quasi-neutrality together require

[TABLE]

where $e$ is the elementary charge, the charge density $\rho_{e}=\sum_{i}z_{i}ec_{\ell,i}$, and $\mathbf{i}_{\ell}$ is the current density in the electrolyte, related to a sum of ionic fluxes,

[TABLE]

We will relate fluxes to both concentrations and electrostatic potentials in the electrolyte, so with constitutive flux relationships, Eqs. 1 and 2 fully specify the system for both the set of concentrations and the electrostatic potential field. However, it is convenient to use quasi-neutrality to eliminate one of the species conservation equations using (with arbitrary $n$)

[TABLE]

which allows us to neglect Eq. 1 for one species and post-calculate the missing concentration profile. For example, for the case of a binary electrolyte of cations, $+$, and anions, $-$, we define the neutral salt concentration, $c_{\ell}=\frac{c_{\ell,+}}{\nu_{+}}=\frac{c_{\ell,-}}{\nu_{-}}$ and instead simulate one of

[TABLE]

and

[TABLE]

In the case of Li-ion batteries for which we will assume a binary electrolyte in which only the Li+ ions react, $R_{V,-}=0$, so it is particularly convenient to simulate the anion species conservation equation.

The boundary conditions for a porous electrode simulation relate the fluxes of the simulated species and current at the anode and cathode current collectors and depend on the simulation. At the current collectors of porous electrodes, $\widehat{\mathbf{n}}\cdot\mathbf{F}_{\ell,i}=0$ and $\widehat{\mathbf{n}}\cdot\mathbf{i}_{\ell}=0$. However, if a foil (e.g. Li metal electrode) is used for one electrode, the boundary condition at that side is replaced by $\widehat{\mathbf{n}}\cdot\mathbf{i}_{\ell}=i_{\mathrm{cell}}$, where $i_{\mathrm{cell}}$ is the macroscopic current density of the cell, Eq. 60. The species flux boundary conditions depend on which species is being eliminated in the set of species conservation equations, but for the case of anion conservation in a Li-ion battery, $\widehat{\mathbf{n}}\cdot\mathbf{F}_{\ell,-}=0$ at all current collectors (neglecting side reactions).

The transport in the electrolyte can be simulated using either simple dilute Nernst-Planck equations or a concentrated solution model based on Stefan-Maxwell coupled fluxes. In both cases, we will neglect convection and assume binary electrolytes. Both formulations relate a species flux, $\mathbf{F}_{\ell,i}$ to gradients in electrochemical potentials, $\mu_{\ell,i}$. The electrochemical potential of a given ion generally has both chemical and electrostatic contributions, and these can be separated a number of ways. For example, it can be separated using an “inner” or Galvani potential Vetter (1967); Newman and Thomas-Alyea (2004),

[TABLE]

where $k_{\mathrm{B}}$ is the Boltzmann constant, $T$ is the absolute temperature, $a_{\ell,i}$ is the activity of species $i$, $\mu_{\ell,i}^{\Theta}$ is its reference chemical potential, and $\phi_{\ell}$ is the (Galvani) electrostatic potential. The final expression serves as the definition of the excess chemical potential, $\mu_{\ell,i}^{\mathrm{ex}}$, and $\widetilde{c}_{\ell,i}$ is the concentration, $c_{\ell,i}$, scaled to some suitable reference. The excess chemical potential contains all of the entropic and enthalpic contributions to the free energy of a concentrated solution beyond that of a dilute solution of a neutral species, such as short-ranged forces, long-ranged electrostatics, and excluded volume effects for finite sized ions Bazant et al. (2009); Liu and Monroe (2014).

An alternative approach common in battery modeling Doyle et al. (1993); Thomas et al. (2002); Verbrugge and Koch (2003) for separating these contributions is to define the electrostatic potential in the electrolyte as that measured with a suitable reference electrode at a position of interest in the electrolyte with respect to another reference at a fixed position in the solution. This has the advantage, unlike Eq. 7, of being defined entirely by measurable quantities, unlike the activities of individual ions in solution. For example, in a lithium ion battery, this is typically done with a Li/Li+ reference electrode. We will refer to this potential as $\phi_{\ell}^{r}$, and we note that it is actually measuring (a combination of) full electrochemical potentials of ions in solution. For the Li/Li+ example,

[TABLE]

Quasi-electrostatic potentials, $\phi_{\ell}^{q}$, can also be used in electrolyte models Newman and Thomas-Alyea (2004), and they are defined by the excess chemical potential of a particular ion in solution, $n$,

[TABLE]

where $\gamma_{\ell,i}=\frac{a_{\ell,i}}{\widetilde{c}_{\ell,i}}$ is the activity coefficient of species $i$.

We also note here that the electric potential field is determined in electrolytes either (1) by an assumption of quasi-neutrality, i.e. $\sum_{i}z_{i}c_{\ell,i}=0$ everywhere, or (2) by solving the Poisson equation, \boldsymbol{\nabla}\cdot\mathopen{}\mathclose{{}\left(\varepsilon\boldsymbol{\nabla}\phi_{\ell}}\right)=-\rho_{e} with permittivity $\varepsilon$, which could depend on concentration or electric field Bazant et al. (2009), or capture non-local ion-ion correlations as a differential operator Bazant et al. (2011). Throughout this work, we assume quasi-neutrality but make a few comments about the alternative here. Use of the Poisson equation enables physical boundary conditions on the electric potential such as specified surface charge densities at interfaces, and the resulting electric potential is the potential of mean force acting on a test charge in the solution. It captures double-layers of diffuse charge with thickness characterized by the Debye length at interfaces, outside of which the net charge approaches zero. The assumption of quasi-neutrality leads to a different potential field which cannot capture the effects of electric double layers at interfaces. Of note, the form of the chemical potential which can most easily accommodate models with both the quasi-neutral and the Poisson equations is the Galvani potential. Of course, this can lead to further confusion because the potential field obtained by assuming quasi-neutrality need not satisfy the Poisson equation with zero charge density and generally will not; rather, quasi-neutrality in this case is the result of the “outer” solution to a singular perturbation in which the Debye length approaches zero. Thus, we could further distinguish between Galvani potentials obtained via each of those methods, but refrain from complicating the notation here, as the quasi-neutral Galvani potential closely resembles the Poisson-satisfying Galvani potential in the large majority of systems with dimensions much larger than the Debye length.

(Semi-)Dilute Electrolyte

In the simpler model, we neglect couplings between species fluxes, such that the flux of a given ionic species is a function only of gradients in its own electrochemical potential,

[TABLE]

where $M_{\ell,i}$ is the species’ mobility. Using the Galvani potential form in Eq. 7,

[TABLE]

where we have used the Einstein relation between the diffusivity in free solution, $D_{\ell,i}$ and mobility, $D_{\ell,i}=M_{\ell,i}k_{\mathrm{B}}T$, non-dimensionalized the potential by the thermal voltage at some reference temperature, $\widetilde{\phi}_{\ell}=\frac{e\phi_{\ell}}{k_{\mathrm{B}}T_{\mathrm{ref}}}$, and defined a non-dimensional $\widetilde{T}=T/T_{\mathrm{ref}}$. With D_{\ell,\mathrm{chem},i}=D_{\ell,i}\mathopen{}\mathclose{{}\left(1+\frac{\partial\ln\gamma_{\ell,i}}{\partial\ln\widetilde{c}_{\ell,i}}}\right) and uniform temperature,

[TABLE]

In a porous medium we use effective transport properties, which are adjusted by the tortuosity of the electrolyte phase, $\tau$. In addition we define the flux per area of porous medium rather than per area of electrolyte requiring a prefactor of the porosity of the electrolyte phase, $\epsilon$.

[TABLE]

The tortuosity is often described as a function of the porosity, the volume fraction of the electrolyte phase, $\epsilon$, commonly by employing the Bruggeman relation $\tau=\epsilon^{a}$ Bruggeman (1935). The value of $a$ is often set to $-0.5$ but can be adjusted Shen and Chen (2007) to account for experimentally Thorat et al. (2009); Shearing et al. (2010); Wilson et al. (2011); Ender et al. (2011); Ebner et al. (2014) or theoretically Vijayaraghavan et al. (2012); Torquato (2013); García-García and García (2016) observed departures from the original derivation. Eq. 13 with Eqs. 3, 5 and 6 define the (semi-)dilute electrolyte model. Although the electrolyte transport model developed in the following section is more reasonable for battery models, we present and retain the (semi-)dilute model here for a number of reasons. Retaining it facilitates comparisons between the models, and the dilute model is easier to extend for electrolytes with more components, even if doing so loses the information related to the extra transport parameters associated with Stefan-Maxwell transport theories. In addition, as mentioned above, it is straightforward to connect the Galvani potential used here to extensions using the Poisson equation to investigate behaviors at interfaces.

Stefan-Maxwell Concentrated Electrolyte

The formulation above assumes that gradients in the electrochemical potential of species $i$ lead to fluxes only of species $i$. However, the framework can be generalized by assuming that fluxes of a given species are related to gradients in electrochemical potentials of each species in the system,

[TABLE]

where $U_{ij}$ are the direct ($i=j)$ and indirect ($i\neq j$) transport coefficients Balluffi et al. (1954). This formulation has been used to describe relatively concentrated electrolytes and is the most commonly used model in battery simulation Newman and Thomas-Alyea (2004); Thomas et al. (2002). Noting that not all electrochemical potentials are independent (from the Gibbs-Duhem relationship), the above can be reorganized to the more commonly written form in terms of species velocities, $\mathbf{v}_{i}=\mathbf{F}_{\ell,i}/c_{\ell,i}$, Newman and Thomas-Alyea (2004)

[TABLE]

with $c_{T}=\sum_{i}c_{\ell,i}$ and $K_{ij}=K_{ji}$. For a binary electrolyte in a porous medium with cations, anions, and solvent (denoted as species [math]), assuming uniform temperature and that the solvent concentration varies only negligibly with salt concentration, and again neglecting convection,

[TABLE]

where, defining $\gamma_{\ell,\pm}^{\nu}=\gamma_{\ell,+}^{\nu_{+}}\gamma_{\ell,-}^{\nu_{-}}$ with $\nu=\nu_{+}+\nu_{-}$,

[TABLE]

The current density is given by

[TABLE]

with

[TABLE]

The values of $s_{i}$ are specified by the choice for the reference electrode with reaction

[TABLE]

For lithium-ion batteries, the typical choice for the reference electrode defining $\phi_{\ell}^{r}$ is Li/Li+, so $s_{+}=-1$, and $s_{-}=s_{0}=0$. Thus, Eqs. 17 and 21 with Eqs. 5 and 6 define the electrolyte model when using the Stefan-Maxwell concentrated solution theory with a binary electrolyte.

II.1.2 Solid phase electronic model

The solid phase of a porous electrode is composed of several length scales and percolating phases. The active material stores the reduced species (e.g. lithium), conductive additive improves electronic wiring, and binder is added to keep all the components connected. Lithium transport occurs within individual particles, and electrons must reach the surface of those particles via traveling over the length of the electrode. These equations approximately capture the variation in electric potential in the matrix of conductive material, which may be assumed to be identical to that of (well connected) active material. Otherwise, additional relations can be added to describe losses between the conductive matrix and the active materials. Over the length scale of the electrode, we describe conservation of charge as in the electrolyte,

[TABLE]

or

[TABLE]

for a binary electrolyte in which cations and/or anions may undergo electrochemical reactions, and where $\mathbf{i}_{s}$ is the current density in the solid phase. The sign difference compared to Eq. 6 comes from the observation that charge entering the liquid phase must be leaving the solid phase. The current density in the solid phase is given by assuming a bulk Ohm’s law,

[TABLE]

where $\sigma_{s}$ is the conductivity of the electronically conductive matrix, $\phi_{s}$ is its electrostatic potential, and we have assumed that the volume fraction of the conductive phase is given by the space not occupied by the electrolyte. The boundary conditions for this come from observing no electronic current can flow into the separator, $\widehat{\mathbf{n}}\cdot\mathbf{i}_{s}=0$, and the potential at the current collector is specified by the operating voltage on the system in the macroscopic equations, $\phi_{c}$ or $\phi_{a}$ ($c$ stands for cathode while $a$ stands for anode). The Bruggeman relation can again be used to estimate the tortuosity in terms of the volume fraction of the conductive matrix.

To avoid the computational and practical difficulties of simulating full microstructures, we describe the behavior of particle interactions with a small number of representative particles, both along the length of the electrode and also in parallel with each other in terms of electrolyte access to capture the effects of particle size distributions. The particles at the same electrode position (interacting in common with the local electrolyte) could be in parallel or in series electronically where parallel wiring would describe them each having direct access via a single resistance to a conductive network and series wiring might describe a comb structure in which some particles (perhaps the edge of a secondary particle) are connected to the conductive backbone, but electrons must pass through poorly conducting particles to get to particles without good contact to the conductive backbone, similar in concept to the hierarchical model of Dargaville and Farrell Dargaville and Farrell (2010). We demonstrate use of the parallel case with a distribution of contact resistances in ref. Thomas-Alyea et al. (2017). In the second case of series wiring we implement a simplified version of that developed by Stephenson et al. Stephenson et al. (2007) by imposing a finite conductance between particles in series, indexed by $k$, within a simulation volume, $j$,

[TABLE]

where $G_{j,k}$ is the conductance and $I_{j,k}$ is the current between particle $k$ and $k+1$. From charge conservation,

[TABLE]

where $j_{j,k+1}$ is the intercalation rate into particle $k+1$.

II.2 Single Particle Equations

A single particle interacts with the electrolyte via an electrochemical reaction, leading to intercalation of neutral species into the solid phase. However the electrochemistry is modeled (see Section II.2.1), the reaction serves as a source/removal of species into/from the particle. We will describe several different solid models here. Generally, we begin by postulating a free energy functional describing the important physics of the particle,

[TABLE]

where $G$ is the total system free energy, $V_{p}$ is the particle volume, $g$ is the free energy density, $A_{S}$ is the particle surface area, and $\gamma_{S}$ is the surface energy. Typically, we will separate the free energy density into homogeneous, $g_{\mathrm{h}}$ and non-homogeneous, $g_{\mathrm{nh}}$, contributions,

[TABLE]

where the remaining terms could describe the stress state of the system Cogswell and Bazant (2012); Bazant (2013) or other energetic contributions García et al. (2004). Following van der Waals Rowlinson (1979) and Cahn and Hilliard Cahn and Hilliard (1958), we use a simple gradient penalty term to describe the non-homogeneous free energy,

[TABLE]

where $\kappa$ is a gradient penalty tensor (assumed to be isotropic here such that $\boldsymbol{\kappa}=\kappa\mathbf{1}$ and $\mathbf{1}$ is the second-order identity tensor) related to interfacial energy between phases, $c_{i}$ is the concentration of species $i$ within a single particle, and $c_{s,\mathrm{ref}}$ is a suitable concentration scale for the insertion species in the active material. The diffusional chemical potential can then be obtained form a variational derivative of the free energy,

[TABLE]

II.2.1 Electrochemical Reactions

The electrochemical reaction can be described by a number of different models, such as the empirical Butler-Volmer equation Bard and Faulkner (2001) or quantum-mechanical models based on Marcus kinetics Marcus (1993); Chidsey (1991); Zeng et al. (2014), which must be consistently generalized for concentrated solutions in nonequilibrium thermodynamics Bazant (2013); Smith et al. (2017b). We describe here electrochemical reactions of the form

[TABLE]

and we will describe the reactions as a function of the activation overpotential,

[TABLE]

where $\mu_{R}=\sum_{j}s_{j,R}\mu_{j}$ is the electrochemical potential of the reduced state, $\mu_{O}=\sum_{i}s_{i,O}\mu_{i}$ is the electrochemical potential of the oxidized state, and $\mu_{\mathrm{e}}$ is the electrochemical potential of the electrons, which we relate here to the electric potential measured in the conductive matrix, $\mu_{e}=-e\phi_{s}$. $\Delta\mu_{\mathrm{rxn}}$ and $\Delta G_{\mathrm{rxn}}$ indicate the total free energy change of the reaction. In the case of Li+ insertion and the Stefan-Maxwell concentrated electrolyte model using Li/Li+ as a reference electrode, $\mu_{O}=e\phi_{\ell}^{r}$. We adopt the convention here that a positive electrochemical reaction current corresponds to the net rate of reduction. The net reduction current, $i$, can be related to the intercalation flux, $j_{i}$, for a given species by the reaction stoichiometry. For example, for lithium intercalation, $j_{i}=i/e$. Extension to multiple reactions simply involves describing each $j_{i}$ as a sum over the relevant reactions.

It is worth noting that the reaction models currently implemented and discussed below follow the trend in battery modeling to neglect the impact of diffuse charge within double layers on the reaction kinetics. Accounting for this involves using a model of the double layer to adjust the surface concentrations from those outside the double layer to those at the distance of closest approach to the electrode surface (the Stern layer), which actually drive the reaction, as well as accounting for the local electric field driving electron transfer. These changes constitute the Frumkin correction Frumkin (1933); Bazant et al. (2005) to the reaction rate model, and have been recently reviewed Biesheuvel et al. (2009); Yan et al. (2017). Frumkin-corrected Butler-Volmer reaction models have been applied to various electrochemical techniques including steady constant current Itskovich et al. (1977); Kornyshev and Vorotyntsev (1981); Bazant et al. (2005); Biesheuvel et al. (2009), voltage steps Streeter and Compton (2008), current steps van Soestbergen et al. (2010), and linear sweep voltammetry Yan et al. (2017), as well as nano He et al. (2006) and porous Biesheuvel et al. (2011, 2012) electrodes, with clear indication of departure from models neglecting double layers, especially at low salt concentrations with thick double layers (“Gouy-Chapman limit” Bazant et al. (2005); Biesheuvel et al. (2009)). To be used consistently with the models developed here, Frumkin reaction kinetics would need to be extended to concentrated electrolyte solutions, including models of individual ionic activities within the double layers, although Frumkin effects are reduced for very thin double layers at high salt concentration (“Helmholtz limit” Bazant et al. (2005); Biesheuvel et al. (2009)). On the other hand, Frumkin effects that dominate in dilute solutions could be important for practical battery operation at high rates, where severe electrolyte depletion can occur and limit the achievable power density.

Butler-Volmer kinetics

Butler-Volmer reaction kinetics are described by exponential dependence on the activation overpotential, and the net reduction current can be written as

[TABLE]

where $\alpha$ is a symmetry coefficient and $i_{0}$ is the exchange current density, the rate of reaction in the forward and reverse directions when the reaction is in equilibrium. Depending on the system considered, the exchange current density could be modeled as constant Cogswell (2015) or a function of species concentrations or activities. Introduced in ref. Doyle et al. (1996), the effective overpotential, $\eta_{\mathrm{eff}}$, accounts for any film resistance, $R_{\mathrm{film}}$, by

[TABLE]

Bazant and co-workers proposed a form for $i_{0}$ based on reacting species’ activities and a transition state activity coefficient, $\gamma_{\ddagger}$ Bai et al. (2011); Cogswell and Bazant (2012); Ferguson and Bazant (2012), derived by assuming thermally activated transitions in an excess chemical potential energy surface with an electric field across the reaction coordinate contributing to the transition state Bazant (2013),

[TABLE]

where $k_{0}$ is a rate constant, $a_{O}=\prod_{i}a_{i}^{s_{i,O}}$, $a_{e}$ is the activity of the electrons (taken to be unity here), and $a_{R}=\prod_{j}a_{j}^{s_{j,R}}$. The transition state activity coefficient is a postulate about the structure and characteristics of the transition state. For example, for lithium intercalation into LiFePO4, Bai et al. originally proposed \gamma_{\ddagger}={\mathopen{}\mathclose{{}\left(1-c/c_{\mathrm{\max}}}\right)}^{-1} where $c_{\mathrm{\max}}$ is the maximum concentration of lithium within the solid, to indicate the transition state excludes one site. Of note, although the above expression is defined in terms of the activities of individual ions within the electrolyte, these quantities are difficult to directly measure. Because we have not implemented models to estimate ion activities, the software currently assumes $a_{\mathrm{Li}^{+}}=\widetilde{c}_{\ell}$ for both electrolyte models.

It is also common to use an exchange current density form based solely on species concentrations Newman and Thomas-Alyea (2004); Thomas et al. (2002),

[TABLE]

where $\widetilde{c}_{\ell}=c_{\ell}/c_{\ell,\mathrm{ref}}$, $\widetilde{c}_{i}=c_{i}/c_{s,\mathrm{ref}}$, and $\xi_{i}=c_{s,\mathrm{ref}}/c_{i,\mathrm{\max}}$, and the ref subscripts indicate some suitable concentration scale. For the case of lithium insertion, we will choose $c_{s,\mathrm{ref}}=c_{i,\mathrm{\max}}$, so $\xi_{i}=1$.

Marcus-Hush-Chidsey kinetics

Marcus-Hush electron transfer kinetics describe an electron transfer event in terms of a reaction coordinate corresponding to the collective rearrangement of species involved in the transition state Marcus (1956, 1965); Hush (1958); Marcus (1993); Hush (1999); Kuznetsov and Ulstrup (1999). The electron transfer event can occur when it is energetically equivalent to occupy the donor or acceptor species, and the fluctuations which lead to this are related to dielectric rearrangement of molecules and charges around the donor-acceptor pair. The Franck-Condon condition is satisfied because the electron transfer event is much faster than the rearrangements of the nearby molecules contributing to the energy of the reacting species. Following the overview of Fedorov and Kornyshev Fedorov and Kornyshev (2014), the reductive and oxidative currents at an electrode can be calculated as an integral over the energy levels of the electrons in the electron conducting phase (the electrode),

[TABLE]

where $k^{0}$ is some rate prefactor, $c_{O}=\prod_{i}c_{i}^{s_{i,O}}$ and $c_{R}=\prod_{j}c_{j}^{s_{j,R}}$, $z$ represents the energy level of electrons in the electrode, $\rho(z)$ is the density of states, and $n_{e}(z)$ is the Fermi function. $W_{\mathrm{red/ox}}$ are the transition probabilities of the elementary reduction/oxidation electron transfer processes. Marcus theory considers the collective motion and reorganization of molecules near the electron transfer event, causing higher or lower energy contributions to the initial or final states of the reaction event via the electrostatic interactions between nearby polarization and the electron donor/acceptor. By treating this collective motion as approximately parabolic around their lowest energy configurations, the following can be derived for the transition probabilities,

[TABLE]

where $k_{\mathrm{w}}$ is a prefactor with explicit dependence on various factors, including the overlap integrals for the wave functions of the various reaction elements Kuznetsov and Ulstrup (1999). We assume reaction symmetries by the equality of the reductive and oxidative prefactors and treat it as a lumped constant here. $w_{R/O}$ are energies describing the probabilities of the species arriving and orienting at the reaction site. $\lambda$ is the reorganization energy, related to the force constants or the curvature of the reaction parabolas, which we take to be the same for forward and reverse reactions, i.e. symmetric Marcus theory. It is defined as the energy required to perturb the system to the stable configuration of the product without allowing electron transfer. The driving force is the reaction change in excess free energy, $\Delta G_{\mathrm{rxn}}^{\mathrm{ex}}$,

[TABLE]

where we have assumed single-electron transfer, as supported by Marcus theory for elementary reaction events Bazant (2013). The tildes indicate non-dimensional concentrations $\widetilde{c}_{O}=\prod_{i}\widetilde{c}_{i}^{s_{i,O}}$ and $\widetilde{c}_{R}=\prod_{j}\widetilde{c}_{j}^{s_{j,R}}$, with each scaled to its reference concentration ($c_{\ell,\mathrm{ref}}$ for ions in the electrolyte and $c_{\mathrm{s,\mathrm{ref}}}$ for intercalated species in the active material). Thus, $\eta_{f}$ is also the departure of the electrode potential from the formal potential. Interestingly, for non-electrode reactions in solution, this predicts a maximum in reaction rate at driving forces given by $\pm\lambda$, and experimental validation of this so-called “inverted” region of decreasing reaction rate at increasing driving force (i.e. negative differential reaction resistance Bazant (2017)) paved the way for the Nobel Prize of Marcus Marcus (1993).

Assuming $w_{R/O}$ are independent of the electron energy level (not necessarily true Kuznetsov and Ulstrup (1999); Fedorov and Kornyshev (2014)), neglecting variation in the density of states, modifying the driving force to account for the energy levels of electrons along the Fermi distribution, and lumping constants into $k_{M}$, we can arrive at the Marcus-Hush-Chidsey (MHC) electron transfer reaction model Chidsey (1991),

[TABLE]

where $z$ is related to the energy level of the electronic states in the metal, and integration is over the Fermi distribution. Reduction/oxidation correspond to the top/bottom signs.

Curiously, Marcus-style kinetics were not used to describe electron transfer reactions in batteries until very recently Bai and Bazant (2014), possibly in part because of the complexity of the expressions involving improper integrals, making their computational evaluation cumbersome. But recently, various approaches have been developed to facilitate the evaluation of the MHC expression, including an asymmetric variant with different values of $\lambda$ for the initial and final states, referred to as “asymmetric Marcus Hush” (AMH) kinetics Oldham and Myland (2011); Bieniasz (2012); Migliore and Nitzan (2012); Zeng et al. (2014, 2015). We will focus exclusively on the symmetric variant here. As a result, MHC kinetics are now as easy to simulate as Butler-Volmer kinetics, so we include them in this simulation software to facilitate comparison in porous electrode modeling and data fitting. For reviews of (asymmetric) MHC kinetics and applications to experimental data, see refs. Henstridge et al. (2012); Laborda et al. (2013).

We make one final assumption following ref. Smith et al. (2017b) before arriving at the form we will use. The $w_{R/O}$ functions are related to the probabilities of the system arriving at the state described by the parabolic minima in the theory. As in the derivation of the Butler-Volmer expression by Bazant Bazant (2013), we postulate that this could capture effects of concentrated solutions beyond simple probabilistic occupation related to the concentration prefactors above. In other words, they are related to the excess chemical potential of the state corresponding to the parabolic minima. Assuming a symmetric approach to the reacting state from the forward and reverse directions,

[TABLE]

or

[TABLE]

Thus, we arrive at our expressions for MHC kinetics,

[TABLE]

with

[TABLE]

and

[TABLE]

Although the derivation of the Butler-Volmer equation Bazant (2013) and the approach followed above both lead to some prefactor related to the excess chemical potential of the transition state, we suggest that the two terms are not capturing the same physical phenomena and thus may differ in their functional forms. In this approach to the microscopic Marcus theory, there is some separation between the energetic contributions accounted for in the microscopic reorganization and the “approach” contributions lumped into $\gamma_{M,\ddagger}$. The Butler-Volmer derivation does not make this distinction, suggesting that the associated contributions in the two theories need not necessarily be equal. Finally, following Zeng et al. Zeng et al. (2014), we replace Eq. 48 with

[TABLE]

where tildes indicate scaling by the thermal energy or voltage, $k_{\mathrm{B}}T$ or $k_{\mathrm{B}}T/e$, $\operatorname{erfc}(z)=1-\operatorname{erf}(z)$ is the complementary error function, and reduction/oxidation corresponds to the top/bottom sign.

II.2.2 Species Conservation

Conservation of the intercalant within the solid particles should be specialized to describe the particular physics of the material being studied. Several options are described below. We will simulate the individual particles by describing neutral species transport within them and their mass exchange with the electrolyte via the electrochemical reactions. This assumes that electron mobility within the active materials is much larger than that of the inserted species Thomas et al. (2002).

Homogeneous

When transport within the solid particles is fast, it can be computationally beneficial to approximate the particles with an average concentration, $\overline{c}_{s}$ Ferguson and Bazant (2014). Then, given a reacting surface area, $A_{p}$, and volume, $V_{p}$, the dynamics of intercalant $i$ can be described simply by the average intercalation rate from the electrochemical reaction, $\overline{j}_{p}$,

[TABLE]

Allen-Cahn Reaction

For Allen-Cahn reaction particles, the reaction occurs as a volumetric source term. This arises when particles are assumed to be either homogeneous or depth averaged in the direction normal to the reacting surface, as Bazant and co-workers have employed in models of LiFePO4 Bai et al. (2011); Cogswell and Bazant (2012); Li et al. (2014). Then, neglecting transport, the local rate of change of the concentration at particle location $\mathbf{r}$ is

[TABLE]

Cahn-Hilliard Reaction

In this model, transport of the intercalant within the particles is described by species conservation,

[TABLE]

The flux can be modeled using linear irreversible thermodynamics Groot and Mazur (1962), which postulates that fluxes arise from gradients in the diffusional chemical potential,

[TABLE]

where we have again used the Einstein relation. The diffusional chemical potential is scaled to the thermal energy at some reference temperature, $k_{\mathrm{B}}T_{\mathrm{ref}}$, and can be calculated from Eq. 32. Following Bazant Bazant (2013), $D_{i}$ can be rewritten in terms of tracer diffusivity in the dilute limit, $D_{0,i}$, the activity coefficient of species $i$, $\gamma_{i}=a_{i}/\widetilde{c}_{i}$, and the activity coefficient of the diffusion transition state, $\gamma_{\ddagger,i}^{d}$,

[TABLE]

Assuming simple diffusion on a lattice in which the diffusion transition state has similar enthalpic contributions as the diffusing species in lattice sites but in which the transition state excludes two adjacent sites, \gamma_{i}/\gamma_{\ddagger,i}^{d}=\mathopen{}\mathclose{{}\left(1-c_{i}/c_{i,\mathrm{\max}}}\right), giving

[TABLE]

although other effects could be accounted for in $\gamma_{\ddagger,i}^{d}$ including stresses in the transition state Aziz et al. (1991). At the particle surface, the flux is given by the electrochemical reaction, $\widehat{\mathbf{n}}\cdot\mathbf{F}_{i}=-j_{p,i}$ where $\widehat{\mathbf{n}}$ is a unit normal vector pointing from the active material to the electrolyte. The natural boundary condition Cahn (1977); Bazant (2013) imposes a constraint on the concentration at the surface as a function of the surface energy, $\gamma_{S}$, $\widehat{\mathbf{n}}\cdot\frac{\partial g}{\partial\boldsymbol{\nabla}\widetilde{c}_{i}}=\widehat{\mathbf{n}}\cdot\kappa\boldsymbol{\nabla}\widetilde{c}_{i}=\frac{\partial\gamma_{S}}{\partial\widetilde{c}_{i}}$.

Solid Solution

If the free energy can be described as a function of only the concentration (disregarding the effect of gradients and other contributions), then Eq. 54 can be rewritten as

[TABLE]

where $D_{\mathrm{chem},i}=D_{i}c_{i}\frac{\partial\widetilde{\mu}_{i}}{\partial c_{i}}$. Here, it is sufficient to prescribe the flux at the surface, as in the Cahn-Hilliard reaction model, $\widehat{\mathbf{n}}\cdot\mathbf{F}_{i}=-j_{p,i}$.

II.3 Coupling Equations

The general approach we take to simulate the two coupled phases is to simulate, within the same physical and simulated space as the electrolyte, a representative sample of particles, which are duplicated to the appropriate filling fraction of solids within the electrode. Instead of simulating discrete particles, an alternative approach could involve directly simulating the distribution of particles at a given state using a population balance model Kumar and Ramkrishna (1996); Marchisio and Fox (2005); Ramkrishna (2000). This may improve accuracy at fixed computational cost, especially for simple particle models like the homogeneous approximation. However, for the more complicated particle models which explore only a tiny fraction of their state space, the current approach may be more efficient.

For example, to simulate an electrode with a very narrow particle size distribution, we simulate one particle at each electrode position. For wide distributions, we use multiple particles sampled randomly from a given input distribution, and at each location, multiple simulated particles interact with the same electrolyte. The representative simulated particles at each electrode position are scaled to occupy the specified volume fraction of active material per electrode volume. We allow the simulated particle size distribution to vary as a function of position, as this allows us to sample a broader distribution within the full electrode, but loses accuracy in situations with significant transport losses in the bulk electrolyte or electron-conducting phases. When setting up simulations, the primary consideration is to have a good representation of the full particle size distribution within the macroscopic length scales of interest, which may be the full electrode at low currents or small regions at high currents with strong electrolyte depletion or electronic losses.

The two phases are coupled via the electrochemical reactions. The volumetric reaction rate of species $i$, $R_{V,i}$, is related to the electrochemical reactions and the flux of that species out of the solid particles. This can either be done by integrating the reaction rate over particle surfaces or by applying the divergence theorem to relate that to their average rate of filling. For particle $p$, at a given position in the electrode,

[TABLE]

Thus, the volumetric source term can be cast as a sum over all the particles at a given position,

[TABLE]

where $P_{L}$ is the loading percent of active material in the solid phase, $V_{u}=\sum_{p}V_{p}$, and $p$ indicates properties of particle $p$, and the summation is over particles at the location where $R_{V,i}$ is evaluated.

II.4 Macroscopic Equations

The overall current density per electrode area, $i_{\mathrm{cell}}$, is defined as the integral of the net charge consumed by the reactions in the electrodes per unit cross-sectional area,

[TABLE]

where $L_{\mathrm{a}}$ and $L_{\mathrm{c}}$ are the lengths of the anode and cathode. We also find it useful to define a current in terms of a C-rate, which is determined by the electrode with the smaller capacity. The capacity of the cell per area, $Q_{A}$, is given for the case of an electrode with a single intercalating species with maximum concentration $c_{\mathrm{\max},k}$ in electrode $k$ by the capacity per area of the limiting electrode,

[TABLE]

With that, we can define a C-rate current,

[TABLE]

where $\tau_{\mathrm{hr}}$ is a conversion factor to ensure units of $\mathrm{hr}^{-1}$.

The overall utilization (state of charge) of electrode $k$ for species $i$, $u_{k,i}$, is defined in terms of the average filling fraction of all the particles in the electrode,

[TABLE]

where $k$ indicates either the anode or cathode and the summation is evaluated as a function of position, given the selection of particles in the simulated electrode at that position. From above, $\xi_{i}=c_{s,\mathrm{ref}}/c_{i,\mathrm{\max}}$, and $\widetilde{\overline{c}}_{p,i}=\overline{c}_{p,i}/c_{s,\mathrm{ref}}$, so the product is the average filling fraction of particle $p$.

The overall cell voltage, $\Delta\phi_{\mathrm{cell}}=\phi_{c}-\phi_{a}$ is defined as the difference in the electric potential at the cathode and anode current collectors. There is an arbitrary datum for the potential, and we set $\phi_{c}=0$ in the simulations. We account for series resistance by defining an applied voltage, $\Delta\phi_{\mathrm{appl}}$, by

[TABLE]

where $R_{\mathrm{ser}}$ is the area specific resistance of the cell.

III Non-dimensional Equations

All times in the simulation are scaled to a reference time scale $t_{\mathrm{ref}}$. We choose some representative diffusive time scale within the cell, $t_{\mathrm{ref}}=L_{\ell,\mathrm{ref}}^{2}/D_{\ell,\mathrm{ref}}$ where $L_{\ell,\mathrm{ref}}$ is a characteristic length within the cell (we will use the cathode length) and $D_{\ell,\mathrm{ref}}$ is a suitably chosen scale for the electrolyte diffusion coefficient.

III.1 Electrode Scale Equations

III.1.1 Electrolyte Model

Defining a reference electrolyte concentration, $c_{\ell,\mathrm{ref}}$ (e.g. $c_{\ell,\mathrm{ref}}/N_{\mathrm{A}}=1\ \mathrm{M}$, where $N_{A}$ is the Avogadro constant),

[TABLE]

with $\widetilde{\boldsymbol{\nabla}}=L_{\mathrm{ref}}\boldsymbol{\nabla}$ and tildes indicating non-dimensionalization by the following scales

[TABLE]

The (semi-)dilute fluxes are non-dimensionalized by

[TABLE]

with diffusive scale as that chosen to define $t_{\mathrm{ref}}$ above. The binary concentrated solution theory can be non-dimensionalized with

[TABLE]

with conductivity scale

[TABLE]

III.1.2 Solid phase electronic model

The solid phase current density and conductivity are scaled to $i_{\mathrm{ref}}$ and $\sigma_{\mathrm{ref}}$ respectively such that

[TABLE]

and

[TABLE]

III.2 Single Particle Equations

In the solid particles, we will use the same time scale as in the electrode scale equations but different length and concentration scales,

[TABLE]

where $L_{s,\mathrm{ref}}$ and $c_{s,\mathrm{ref}}$ are relevant scales to the particle. We will use $L_{s,\mathrm{ref}}=L_{p}$, where $L_{p}$ is a characteristic length scale of the active material particle and which may vary by particle. For the case of lithium insertion electrodes with a single intercalating species with maximum concentration given by $c_{\mathrm{\max}}$, we choose the solid reference concentration to be the maximum filling for lithium insertion cells, $c_{s,\mathrm{ref}}=c_{\mathrm{\max}}$.

III.2.1 Electrochemical Reactions

The various forms of the overpotential are all scaled to the thermal voltage, $\frac{k_{\mathrm{B}}T_{\mathrm{ref}}}{e}$, and the rate prefactors and current density are both scaled to

[TABLE]

The film resistance is scaled to $R_{\mathrm{film,ref}}=\frac{k_{\mathrm{B}}T_{\mathrm{ref}}}{ei_{s,\mathrm{ref}}}$, so

[TABLE]

and the Butler-Volmer expression is

[TABLE]

while the MHC expression is

[TABLE]

III.2.2 Species Conservation

We scale the rate of intercalation of each species in particle $p$, $j_{p,i}$, to the active material reference flux, $N_{s,\mathrm{ref}}$, such that, for homogeneous particles,

[TABLE]

where $\widetilde{\overline{c}}_{i}=\overline{c}_{i}/c_{s,\mathrm{ref}}$, and

[TABLE]

For Allen-Cahn reaction particles,

[TABLE]

where $\widetilde{\mathbf{r}}$ is non-dimensionalized by $L_{p}$. For Cahn-Hilliard reaction particles,

[TABLE]

The non-dimensional flux is given by

[TABLE]

where the diffusivity is scaled to

[TABLE]

For the case of diffusion on a lattice with a simple excluded site model for the transition state, Eq. 56 becomes

[TABLE]

The natural boundary condition is scaled using

[TABLE]

such that

[TABLE]

at the particle surface, and the flux boundary condition is simply

[TABLE]

For solid solution particles,

[TABLE]

with the same flux boundary condition as the Cahn-Hilliard reaction particles.

III.3 Coupling Equations

Here, we retain the same scales as above, such that

[TABLE]

and

[TABLE]

with $\beta=c_{s,\mathrm{ref}}/c_{\ell,\mathrm{ref}}$.

III.4 Macroscopic Equations

Keeping the same scales as in the electrode model,

[TABLE]

and

[TABLE]

where $k$ indicates either the anode or cathode, $\widetilde{L}_{k}=L_{k}/L_{\ell,\mathrm{ref}}$, and $\widetilde{V}_{p}=V_{p}/V_{u}$. The series resistance can be scaled by the electrode model scales, $R_{\mathrm{ser,ref}}=\frac{k_{\mathrm{B}}T_{\mathrm{ref}}}{ei_{\mathrm{ref}}}$, and the potentials can all be scaled by the thermal voltage at the reference temperature, $k_{\mathrm{B}}T_{\mathrm{ref}}/e$, so

[TABLE]

IV Model implementation

To solve the system of PDE’s in the model, we take the general approach of discretizing each in space using some variant of the finite volume method to obtain a system of differential algebraic equations (DAE’s), and then stepping in time using a variable-order adaptive time stepper. We discretize in space using finite volume methods both for their robustness to steep gradients and also their mass conservation to within numerical accuracy Eymard et al. (2007); Moukalled et al. (2016). We use the IDAS time stepper of the SUNDIALS integration suite Hindmarsh et al. (2005) to solve the resulting DAE’s with a backward differentiation formula approach. The model and simulation are defined within the DAE Tools software package Nikolić (2016), which provides a modeling language environment within Python similar to that of more specialized modeling languages like gProms Barton and Pantelides (1994) or Modelica Fritzson and Engelson (1998). This allows use of the full depth of the general purpose Python environment while adding helpful concepts from modeling languages such as logical separation of model definition from simulation setup and equation oriented model definition.

For example, to define a model, the user writes the spatially discretized equations in a familiar form within a model class and defines simulation particulars such as parameter values and initial conditions elsewhere. The DAE Tools software handles the formation of all underlying system matrices and interactions with other numerical libraries involved in the actual time advancement. As an extra convenience, it wraps IDAS with the ADOL-C automatic differentiation library Griewank et al. (1996) to form the analytical-accuracy Jacobian matrix, greatly facilitating the solution of non-linear systems of equations involved in the implicit time stepping. This eases additions and modifications to the model, which can be made without user input of analytical derivatives or reliance on numerical approximation of the Jacobian. This approach enables MPET developers to work entirely within a Python environment while using the fast and vetted lower level numerical libraries for computationally expensive or analytically tedious aspects of the simulation.

IV.1 Software Organization and Structure

To make MPET flexible enough to simulate arbitrary active material models in the context of either porous electrodes or single particles to study their dynamics directly, we chose to develop the software to logically and structurally isolate the material models from the overall cell model using an object oriented framework, see Figure 2. Objects containing local information about active material models exchange only the required information with the cell model, reducing assumptions built into each section of the software and leading to a more modular, extensible structure. Parameters are specified within input files using a standard config-file syntax for ease of scripting multiple simulations. These input files are split according to information used by the overall system (e.g. specified current profile and system dimensions) and those defining properties of the simulated active materials (e.g. material thermodynamics, transport, and reaction properties). This enables a user to reference standard material property files while editing only a system input file related to cell design. These inputs are then non-dimensionalized for the simulation, as described in Section III, and passed to the system model.

Depending on the simulation to be carried out (e.g. perfect electrolyte bath, half-cell, full-cell, with or without particle size variability), the system model creates as many instances of the appropriate active material model as necessary for the representative active material particles and establishes communication between them to only exchange key information. The active material particles get access only to local concentrations and potentials in the bulk phases near the particular active material particle, and the bulk surrounding phases only need access to the total integrated reaction rate of each particle. This isolation makes it relatively straightforward to extend MPET’s capabilities by adding new material models or modifying existing ones to add relevant physics such as stresses and strains García et al. (2004); Christensen and Newman (2006); Van der Ven et al. (2009); Cogswell and Bazant (2012); Bazant (2013); Sethuraman et al. (2012); Ulvestad et al. (2014) or hierarchical structures Dargaville and Farrell (2010). For example, to add and simulate a new model, a user must define the non-dimensional equations as a model class, specify the initial conditions via the simulation class, and update the handling of parameters to read and non-dimensionalize any new inputs. Material models can easily use distinct discretization and/or geometries such as the homogeneous or CHR particles. We have also already included the two-repeating-layers model developed and applied to single graphite particles Guo et al. (2016); Smith et al. (2017a). Other variations such as 2D particles of arbitrary geometry discretized using the finite volume Dargaville and Farrell (2015) or finite element method would also be straightforward, especially with the DAE Tools interface to the deal.ii finite element library Nikolić (2016); Bangerth et al. (2007). The structure also makes it a straightforward addition to implement multiple particle types to be simulated within each electrode region. Finally, it makes it easier to modify specific parts of any given model with confidence that side-effects will be minimized through the isolation of the logical parts.

In order to facilitate good reproducibility of scientific computations, MPET defaults to storing each simulation output within time-stamped directories, along with both input files and a snapshot of the source code which ran that simulation. Users can disable this feature to use their own system or take advantage of it by keeping a log matching the time-stamped directories with notes about each simulation. Simulation outputs are stored by default in a binary format which is readable using common scientific computing software including Python (with SciPy), R, and MATLAB. MPET includes a script to perform some basic plotting using this output and also a script to convert the output to comma-separated value (CSV) text files, which can be nearly universally interpreted.

IV.2 General Options

As discussed above, the software’s structure makes it flexible enough for a range of possible simulation options, and the key options are described here. Either the current or voltage can be imposed as piecewise-constant functions with arbitrary numbers of steps (simulated with fast but continuous steps, similar to the method proposed in ref. Boovaragavan et al. (2010)). With minor changes to the code, completely arbitrary functions are possible. Currents are imposed relative to a $1\ \mathrm{C}$ (dis)charge based on the simulated battery’s limiting electrode, and optional cut-off voltages terminate the simulation if the system voltage reaches their values. Because it is sometimes convenient to continue an old simulation, this can be done by specifying the location of a stored output. Discretization is specified in terms of the number of finite volumes in each region of the electrode. If the anode is set to have zero volumes, the simulation uses a flat lithium metal counter-electrode with specified Butler-Volmer reaction rate kinetics. If only a single cathode volume is simulated, it places any simulated particles within a perfect electrolyte bath without any transport limitations for either single- or multi-particle studies neglecting electrolyte and bulk electrode losses. Within each discretized porous electrode volume, the chosen number of particles are simulated, drawing from a specified log-normal size distribution. Conductivity losses in the bulk electron-conducting phase or between individual particles within each volume of simulated electrodes are optionally neglected or simulated as described in Section II.1.2. Electrode lengths, electrolyte volume fractions, loading percents (volume fraction of active material within the electron-conducting phase), and Bruggeman exponents describe the geometry characteristics of effective transport within the porous electrolyte. The electrolyte can be specified either as a dilute model or using a full Stefan-Maxwell concentrated solution theory model as described in Section II.1.1. Arbitrary functions can be specified for the electrolyte transport properties by defining them as Python functions.

Active material particles are specified by the basic model for their transport processes (e.g. Fickian diffusion, homogeneous, Cahn-Hilliard reaction), the thermodynamic function describing them (also specified as arbitrary Python functions), and the remaining properties associated with the particular model. The specified transport model determines the type of material model (defining the discretization method and equations solved) used in the simulations. When applicable, the transport flux prefactor can be specified as an arbitrary function. Electrochemical reactions are defined as Butler-Volmer, Marcus-Hush-Chidsey, or (experimental) Marcus kinetics as formulated by Bazant Bazant (2013), and a reaction film resistance can be added to any of these. The Butler-Volmer exchange current density is taken as one of a few built-in options, as outlined in Section II.2.1, but can be specified as arbitrary functions within the particle model source code. To facilitate stability of some models, we provide the options to replace the most extreme regions of the log terms of ideal-solution parts of thermodynamic models with linear extrapolations. Finally, to ensure that materials phase separate when they naturally would based on minor fluctuations, we include the option to add Langevin noise to the rate of change of concentration within each solid finite volume Cogswell and Bazant (2012).

IV.3 Numerical Methods

For each of the methods detailed below, we will describe the discretization in non-dimensional form using the scales presented in Section III. We will also drop subscripts indicating species and location (electrolyte or active material particle) for clarity with discretization subscripts.

IV.3.1 Electrolyte

The electrolyte equations are discretized using a typical 1D finite volume scheme with cell centers Moukalled et al. (2016). They are non-dimensionalized using the same scales as in the electrode scale equations, Section III.1. For example, for species conservation in cell $j$ with width $\Delta\widetilde{x}$ and neglecting subscripts for electrolyte/species identifiers for clarity,

[TABLE]

where the $j$ subscript represents an average quantity within the discretized volume, and the $F_{j\pm 1/2}$ represent components of the flux normal to the faces on the right/left of cell $j$ in the positive $x$ direction, where $x$ is along the length of the cell going between the current collectors. The flux is approximated using a finite difference two-point formula based on the adjacent cell centers. Where required (e.g. for electromigration terms), face values of field variables such as concentration are approximated using harmonic means, which greatly enhances stability in regions of strong electrolyte depletion. The harmonic mean also better represents variation in transport prefactors Russell and Wheeler (1983).

IV.3.2 Active Material

The active material can be defined as spherical, cylindrical, or a rectangular grid approximation of platelet-like particles common in LiFePO4 Smith et al. (2012). Those on a rectangular grid are discretized like the electrolyte, and the ratio of reacting area to volume is calculated as a function of the length and thickness, $h_{p}$, of the particle assuming reactions only occur on the top and bottom surfaces Li et al. (2014). For spheres and cylinders, the ratio of particle volume to reacting area is specified fully by the particle radius. For cylinders, we use the particle thickness, $h_{p}$, for clarity in this section. The spherical and cylindrical particles are discretized using a variant of finite volumes directly taken from Zeng et al. Zeng et al. (2013); Zeng and Bazant (2014). They are non-dimensionalized using the same scales as in the single particle scale equations, Section III.2. For the cylindrical particles, we follow the same method as that in Zeng and Bazant, but modify it for the cylindrical geometry by changing the calculation of the volumes and areas of each sub-shell. That is, for both geometries, the systems of equations can be represented as

[TABLE]

where $\widetilde{\mathbf{c}}$ is a vector of concentrations at positions going from the center of the particle to the surface of the particle with $N$ sub-volumes, $\widetilde{\mathbf{M}}$ is a tridiagonal matrix given by

[TABLE]

and $\widetilde{\mathbf{V}}$ is a diagonal matrix defined by the volumes of the shells scaled to $L_{p}^{3}$. Indexing the volumes with $j=1\dots N$, and with a uniformly spaced radial vector, $\widetilde{r}_{j}$, with $N$ points going from the particle center to its surface, then for a sphere $\widetilde{\mathbf{V}}$ has components defined by

[TABLE]

and for a cylinder,

[TABLE]

The right hand side, $\widetilde{\mathbf{b}}$, is defined in relation to the radial components of the fluxes evaluated at the interfaces between the shells, $\widetilde{F}_{j\pm 1/2}$, and the areas of the interfaces between the shells, $\widetilde{A}_{j\pm 1/2}$. The fluxes and areas are scaled to $F_{s,\mathrm{ref}}$ and $L_{p}^{2}$ respectively. For each geometry,

[TABLE]

The shell areas differ for the two geometries, given scaled to $L_{p}^{2}$,

[TABLE]

The flux is calculated using the two-point centered finite difference approximation using the concentrations in the adjacent cells. Unlike in Zeng et al. Zeng et al. (2013); Zeng and Bazant (2014), we use harmonic means for face value approximations of field variables as in the electrolyte. This discretization scheme for the cylinders and spheres has the advantage of increased accuracy of simulated surface concentrations Zeng et al. (2013) while retaining mass conservation to within numerical precision. It is slightly less robust than a more typical cell-centered finite volume scheme under certain situations, as the coupling between the adjacent volumes via the mass matrix can cause minor oscillations during the formation of very steep concentration gradients, but this typically only presents problems when beginning with very high currents from nearly full or empty particles. These oscillations could be eliminated by using a flux limiter, commonly used in higher order discretization schemes and especially for hyperbolic problems LeVeque (2002). This approach also has the advantage, much like typical finite volume schemes, of naturally facilitating the application of flux boundary conditions, as applied to the active material at the particle center (from symmetry) and surface (from reaction).

When applied to CHR particles, to evaluate the diffusional chemical potential at grid points we must calculate the Laplacian of the concentration, and we follow Zeng and Bazant Zeng and Bazant (2014). For interior points, $j=2\dots N-1$,

[TABLE]

For the center of the particle at $j=1$, isotropy gives

[TABLE]

At the surface, we use a ghost point at $j=N+1$ and impose the concentration gradient within the Laplacian of the concentration,

[TABLE]

giving $\widetilde{c}_{N+1}\approx\frac{2\Delta\widetilde{r}}{\widetilde{\kappa}}\frac{\partial\widetilde{\gamma}_{S}}{\partial\widetilde{c}_{s}}+\widetilde{c}_{N-1}$, such that

[TABLE]

IV.4 DAE Consistent Initial Conditions

Because the discretized equations are coupled DAE’s, we must begin the system from consistent initial conditions Brown et al. (1998); Ramadesigan et al. (2012). Although various approaches for initializing this type of model have been developed Wu and White (2001); Boovaragavan and Subramanian (2007); Methekar et al. (2011), we found it to be relatively robust to begin from a known equilibrium state and quickly ramp the applied current or voltage to the desired values. Although this contributes to the stiffness of the system, the higher order, adaptive backward difference formula (BDF) time stepper we use mitigates the cost of this, obviating the need to use more sophisticated methods.

To run simulations with specified current, we begin at a zero-current state, where we can easily calculate the equilibrium potentials in each of the phases (because we also begin with uniform concentration profiles). Similarly, for simulations with specified applied voltages, we begin from the calculated equilibrium applied voltage before ramping to the desired voltage. In addition, to make the initialization more robust, we offset the diffusional chemical potentials in the solid phases throughout the simulation by their initial values, such that the initial equilibrium potentials everywhere throughout the simulation are zero. We found this to substantially facilitate initialization compared to leaving diffusional chemical potential at their physical values, which can leave differences in potentials across the system on the order of a hundred thermal volts (e.g. a $3\ \mathrm{V}$ material against a lithium metal or graphite anode is at $\sim 100$ thermal volts compared to $k_{\mathrm{B}}T/e$ at room temperature). Because these potentials appear within exponentials of rate expressions, offsetting them to begin the simulation at zero substantially facilitates the numerical determination of consistent initial conditions and has no impact on the simulation output upon post-applying the reverse-offset. Finally, when performing a simulation continuation, we take the final state of the previous output as initial guesses for the new initial state.

V Examples

As a detailed examination of all the applications and model variations of the software is beyond the scope of the present work, we instead focus here on a few examples highlighting some of the key distinguishing capabilities of MPET.

V.1 Solid solutions and phase separation in porous electrodes

We begin with a comparison of solid solution and phase separating models within porous electrode simulations, which highlights the ability of the software to compare two different approaches to modeling the same material, variations of both of which are commonly used Baker and Verbrugge (2012); Levi et al. (2005). A complete comparison of the different models is beyond the scope of this work, so we present only a simplified example here as a demonstration of the model behavior. To keep the system as simple as possible, we choose to investigate a fictitious material in which lithium intercalates and has a diffusional chemical potential described by a regular solution,

[TABLE]

where $\widetilde{c}=c/c_{\mathrm{\max}}$ is a filling fraction, $c_{\mathrm{\max}}$ is the maximum concentration of lithium in the active material, and with $\mu^{\Theta}$ arbitrarily set to $-2\ \mathrm{eV}$ defined relative to Li/Li+. The regular solution parameter, $\Omega$, is set to $3k_{\mathrm{B}}T_{\mathrm{ref}}$, where $T_{\mathrm{ref}}=298\ \mathrm{K}$ is the absolute temperature at which the simulation is carried out. We examine an electrode with particles with radius, $R=1\ \mu\mathrm{m}$, and we choose the interfacial gradient penalty, $\kappa=1.16\times 10^{-7}\ \mathrm{J}/\mathrm{m}$ such that the phase interface width is several times the finite volume grid size, $20\ \mathrm{nm}$. The maximum filling fraction in the active material, $c_{\mathrm{\max}}=25\ \mathrm{M}$ is chosen arbitrarily.

To describe the thermodynamics of the solid solution material, we impose a diffusional chemical potential associated with the stable equilibrium diffusional chemical potential defined by Eq. 116, a piece-wise continuous function defined as equal to $\mu_{\mathrm{RS}}$ outside the miscibility gap and $\mu^{\Theta}$ within it, the red dashed curve in Figure 3. Note, this could differ from a particular equilibrium filling/emptying path because a particle described by Eq. 116 can enter metastable regions, leading to equilibrium hysteresis Levi and Aurbach (1997); Dreyer et al. (2010); Ferguson and Bazant (2014) which a solid solution model cannot capture without direction-dependent thermodynamic models.

In order to minimize conflating factors from the reaction kinetics, we model both phase separating and solid solution particles using the symmetric ($\alpha=0.5$) Butler-Volmer model with a constant exchange current density, $i_{0}=1\ \mathrm{A}/\mathrm{m}^{2}$, leading to minimal reaction losses in either system. We will consider the case that intercalation occurs at the surface and phase separation occurs within the bulk, indicating it can be described by the Cahn-Hilliard reaction model Bazant (2013); Zeng and Bazant (2014). The transport of solid solution material can be described using Fickian diffusion. For the solid solution, we use constant $D_{\mathrm{chem}}=1\times 10^{-16}\ \mathrm{m}^{2}/\mathrm{s}$, and for the phase separating case, we use Eq. 56, which approximately matches the concentration-independence of the chemical diffusivity in the solid solution regimes Smith et al. (2017a). We adjust $D_{0}$ until we get a similar dependence of available battery capacity as a function of current (Fig. 5 (b)), where capacity is defined as the electrode filling fraction when the voltage reaches a $1.5\ \mathrm{V}$ cutoff. We find $D_{0}=8\times 10^{-16}\ \mathrm{m}^{2}/\mathrm{s}$ reflecting the larger value of transport prefactor required to obtain similar fluxes in phase separating systems because of the smaller gradients which can arise in their end member phases compared to solid solution particles (Fig. 4).

We use a cell geometry with thin electrodes to minimize electrolyte limitations and assume no electron limitations in the porous electrode or reaction resistance on the lithium foil counter electrode. The porous separator and cathode are $15\ \mu\mathrm{m}$ and $20\ \mu\mathrm{m}$ with porosities of $0.8$ and $0.2$ respectively. The loading fraction in the cathode is taken to be 0.7 with the Bruggeman exponent $a=-0.5$. We model the electrolyte using concentrated solution theory as described in Section II.1.1 using parameters as determined by Valo̸en and Reimers Valo̸en and Reimers (2005) but replacing the conductivity with that measured by Bernardi and Go Bernardi and Go (2011).

Despite very different models for the transport in the solid phase, it is interesting to note that some macroscopic characteristics of the modeled cell are quite similar. For example a representative cell voltage polarization curve and the dependence of cell capacity as a function of current show similar behavior (Figure 5), indicating that macroscopic cell data alone would not differentiate between these models. We find in a companion paper that this observation holds for similar models even in the case of non-negligible reaction resistances which are dominated by film resistance Thomas-Alyea et al. (2017). Of note, in the companion paper, we used Eq. 54 with constant prefactor, i.e. $D_{0,i}\gamma_{i}\widetilde{c}_{i}/\gamma_{\ddagger,i}^{d}=\mathrm{constant}$, which we found here unable to produce a dependence of capacity on current similar to the solid solution with constant chemical diffusivity. Despite the similarity in Figure 5, we note that the models give very different predictions of the behavior at the particle level, shown in Figure 4 in which we show both a snapshot of a representative concentration profile for each of the simulations and also a representative profile of the evolution of surface concentration in one of the particles. Because the surface concentration also affects other model behavior, including the reaction rate, we show in Figure 6 that when both particles are simulated using the concentration-dependent exchange current density in Eq. 38, they diverge much more significantly in their predictions. A consistent coupling to stresses and strains would further differentiate the models, as the concentration profile throughout the particle directly couples to the stress profile García et al. (2004); Bower et al. (2011); Cogswell and Bazant (2012); Di Leo et al. (2014).

V.2 Concentrated and dilute solutions with phase separating materials

Here, to focus on the electrolyte, we consider the same (phase separating) electrode materials in dilute and concentrated electrolyte solutions. Dargaville and Farrell first considered CHR and ACR phase field particles in porous electrodes using Stefan-Maxwell electrolytes Dargaville and Farrell (2013), and we briefly compare the dilute approach previously used by Bazant and co-workers Ferguson and Bazant (2012, 2014); Li et al. (2014) with the concentrated model here. We use the same phase separating material as that studied Section V.1 with constant exchange current density but simulate it in a battery with a longer electrode of $200\ \mu\mathrm{m}$. The concentrated electrolyte model is identical to that used in Section V.1, and the dilute model is analogous to that model but replacing the fit function of the electrolyte conductivity with the linear dependence predicted by dilute solutions Newman and Thomas-Alyea (2004), neglecting concentration dependence of the electrolyte transport coefficients, and using $D_{\ell,+}=2.42\times 10^{-10}\ \mathrm{m}^{2}/\mathrm{s}$ and $D_{\ell,-}=3.95\times 10^{-10}\ \mathrm{m}^{2}/\mathrm{s}$, which corresponds to the concentrated solution model with a typical value of $D_{\ell}=3\times 10^{-10}\ \mathrm{m}^{2}/\mathrm{s}$ and $t_{+}^{0}=0.38$.

In this situation, the electrolyte transport is limiting, so we expect substantial differences between the two models, which we see in Figure 7. Because the electrolyte transport is much more limiting than the solid phase transport, the results are qualitatively similar to those of Valo̸en and Reimers despite their use of a solid solution active material model instead of the phase separating model here Valo̸en and Reimers (2005). As in Figure 7 (a), they found the concentrated model leads to a more gradual decrease in voltage over the cell discharge but a larger overall capacity. This highlights the importance of using the more consistent concentrated solution theory approach in situations with substantial electrolyte polarization.

V.3 Full Cell LiFePO4/C

To demonstrate the versatility of the software, we present here the result of a battery simulation with two porous electrodes using two very different solid material models for each one. In the cathode, we use an Allen-Cahn reaction model of LiFePO4 (LFP) which was developed by Bazant and co-workers Burch et al. (2008); Bai et al. (2011); Cogswell and Bazant (2012); Ferguson and Bazant (2014) and used to explain in situ measurements of mosaic phase transformations in porous electrodes Li et al. (2014). In the anode, we simulate graphite using a 2-layer Cahn-Hilliard reaction model Ferguson and Bazant (2014); Smith et al. (2017a) as used to capture the experimental intercalation of a single crystal of graphite Guo et al. (2016). Although we have found a variation of this model to better describe the behavior of porous graphite electrodes Thomas-Alyea et al. (2017), we use it here to highlight the ease of implementing and using very different active material models within the same simulation using MPET.

We use the same concentrated electrolyte model as in the previous two sections. The lengths of the anode, separator, and cathode are $100\ \mu\mathrm{m}$, $20\ \mu\mathrm{m}$, and $150\ \mu\mathrm{m}$. The porosities are $0.15$, $0.4$, and $0.2$, and the loading percents of the anode and cathode are $0.9$ and $0.7$. Because the solid phase models are described in detail in the previous work, we only briefly introduce them here. The cathode particles are ACR particles described by a regular solution with an effective stress term, designed to approximately capture the tilting of the single-particle voltage plateau resulting from stresses Cogswell and Bazant (2012, 2013); Ferguson and Bazant (2014),

[TABLE]

where $\Omega=4.51k_{\mathrm{B}}T$, the stress coefficient $B=0.19\ \mathrm{GPa}$, $c_{\mathrm{\max,LFP}}=23\ \mathrm{M}$, $\overline{c}$ is the average filling fraction in the particle, $\kappa_{\mathrm{LFP}}=5\times 10^{-10}\ \mathrm{J}/\mathrm{m}$, and $\mu^{\Theta}_{\mathrm{LFP}}=-3.4\ \mathrm{eV}$ with respect to Li/Li+. The anode particles are described by two-layer CHR transport with a two-parameter regular solution model including inter-layer repulsion terms to capture the staged structure found in intercalated graphite Dresselhaus and Dresselhaus (1981). The flux within each layer is given by Eq. 56, and the diffusional chemical potential of each is given by

[TABLE]

where $j\neq i$ and i\in\mathopen{}\mathclose{{}\left\{1,2}\right\} and $\widetilde{c}_{i}$ represents the filling fraction of layer $i$. The parameters for graphite come from refs. Ferguson and Bazant (2014); Guo et al. (2016) and are $\Omega_{a}=3.4k_{\mathrm{B}}T$, $\Omega_{b}=1.4k_{\mathrm{B}}T$, $\Omega_{c}=20k_{\mathrm{B}}T$, $c_{\mathrm{\max,G}}=28.2\ \mathrm{M}$, $\kappa_{\mathrm{G}}=4\times 10^{-7}\ \mathrm{J}/\mathrm{m}$, and $\mu^{\Theta}_{\mathrm{G}}=-0.12\ \mathrm{eV}$ with respect to Li/Li+. In both materials, the reactions are governed by the symmetric ($\alpha=0.5$) Butler-Volmer, Eq. 35, using Eq. 37 for the exchange current density. For LFP, we set $k_{0}=0.16\ \mathrm{A}/\mathrm{m}^{2}$ and use $\gamma_{\ddagger}$ = 1/\mathopen{}\mathclose{{}\left(1-\widetilde{c}}\right) as proposed by Bai et al. Bai et al. (2011). For graphite, we set $k_{0}=10\ \mathrm{A}/\mathrm{m}^{2}$ and use \gamma_{\ddagger,i}=1/\mathopen{}\mathclose{{}\left(\widetilde{c}_{i}\mathopen{}\mathclose{{}\left(1-\widetilde{c}_{i}}\right)}\right) as we initially used in Guo et al. Guo et al. (2016) and reexamined recently to suggest alternate models for practical battery simulation Smith et al. (2017a).

We present overall cell polarization curves for a selection of currents in Figure 8. The thick electrodes and low porosity make the cell experience strong electrolyte polarization at the moderate simulated currents, causing most of the losses in the cell. The cell capacity is limited by the cathode, evidenced by the graphite-caused shift in the $\mathrm{C}/10$ voltage profile occurring at a cathode filling fraction larger than $0.5$. In addition, we see erratic profiles at all simulated rates. At the lowest rates, this is primarily related to the mosaic transformation of the LFP particles along the length of the electrolyte, with one representative particle typically receiving almost all the current in the electrode, as seen experimentally Chueh et al. (2013); Li et al. (2014) and explained theoretically Dreyer et al. (2010); Ferguson and Bazant (2014). This behavior can be seen in Figure 9 (b) in which the cathode particles show sharp filling processes related only to their position in the electrode, because they do not have size variation as used in Ferguson and Bazant Ferguson and Bazant (2014). In contrast, the anode particles behave more like the CHR particles studied in the sections above, with more gradual filling processes governed more weakly by their position in the electrode. At higher rates, the cathode still demonstrates a mosaic transformation, but the voltage curves smooth out because other losses dominate. The voltage spikes result from the form of the reaction resistance causing short sections of low resistance. This rate expression helped capture single particle experiments Guo et al. (2016) but is unable to capture experiments on a full porous electrode, where we also used a simplified version of the graphite model Thomas-Alyea et al. (2017).

V.4 Electrochemical reactions

Battery models are typically constructed assuming Butler-Volmer reaction kinetics Botte et al. (2000); Thomas et al. (2002), although recent evidence suggests Marcus kinetics may better describe some electrochemical reactions Chidsey (1991); Bai and Bazant (2014). One possible reason for the use of the Butler-Volmer form is that the models do not deviate until moderate to large reaction driving forces, depending on the reaction Henstridge et al. (2012). A second is that the expression for calculating Marcus style kinetics at an electrode involves an improper integral Chidsey (1991), which is cumbersome to evaluate, especially within porous electrode simulations in which reaction rates are evaluated many times over the course of a simulation. Nevertheless, this Marcus-Hush-Chidsey (MHC) reaction model has been approximated with a mathematically simple and accurate expression Zeng et al. (2014), enabling its practical use in porous electrode battery models, which we demonstrate here. The primary effect of MHC kinetics is a downward curvature instead of the linear Tafel slope associated with Butler-Volmer reaction kinetics.

In order to accurately capture electrochemical data, the Butler-Volmer expression is commonly modified with a film resistance as in Eq. 36 to capture a series resistance associated with any film at the surface. Curiously, this modified expression could be interpreted as a way to introduce curvature to the Tafel plot in a way that can look similar to MHC kinetics over a range of driving forces. For example, in Figure 10, we compare Butler-Volmer with and without film resistance (BV and BV-film respectively) with the MHC expression, all using the same exchange current density and at fixed concentrations. In the figure, we use constant reaction rate prefactors, $i_{0}=i_{M}=1\ \mathrm{A/m}^{2}$ and $R_{\mathrm{film}}=0.001\ \Omega\ \mathrm{m}^{2}$, a value similar to that used to match electrochemical data using porous electrode modeling Thomas and Newman (2003); Srinivasan and Newman (2004); Albertus et al. (2009). The MHC kinetics use the same exchange current density and a reorganization energy, $\lambda=18k_{\mathrm{B}}T$, between that approximated for LiFePO4 Bai and Bazant (2014) and values used in other interfacial reactions Chidsey (1991); Henstridge et al. (2012). The limiting behaviors for the film resistance and MHC models differ substantially, as the film model approaches a linear resistor at high driving forces, whereas the MHC current saturates at a constant value. However, over a relatively broad range of driving forces, the two predict similar currents when neglecting concentration effects, which suggests that some variant of the MHC model may be able to capture some electrochemical data using a microscopically derived model. We should note that the concentration dependence of the MHC model differs from that of the Butler-Volmer model with film resistance, so comparisons are not completely straightforward. Nevertheless, use of the MHC expression has the advantage that it also makes clear connections to surface properties Kuznetsov and Ulstrup (1999) which could be engineered, suggesting possible improvement paths for materials with slow kinetics.

To illustrate the application of MHC kinetics in a porous electrode, we simulate a current-pulse discharge process of the idealized phase separating material studied in Section V.1 using both the Butler-Volmer and MHC models as well as the Butler-Volmer model with a film resistance. As mentioned above, the concentration dependence of the MHC model is non-trivial, as the departure from the formal potential, $\eta_{f}$, is offset from the activation overpotential. In addition, the prefactor, $i_{M}=k_{M}/\gamma_{M,\ddagger}$, may also have concentration dependence related to concentrated solution effects captured in $\gamma_{M,\ddagger}$. In order to study the two with the least impact from the different concentration dependence, we choose to modify $i_{M}$ here to match the exchange current densities of the Butler-Volmer and MHC models. This could be interpreted as choosing a particular model for $\gamma_{M,\ddagger}$, but we refrain from assigning physical meaning to the choice used here and make the selection to enable the simplest comparison between the two models.

First, we plot the exchange current density of the MHC model using constant $i_{M}$ in Figure 11 (a) and note that it has a dependence on the reduced species (intercalated lithium) that is approximately a square root. This applies to both $\widetilde{c}_{R}$ and $\widetilde{c}_{O}$. Thus, in order to approximately specify an arbitrary exchange current density for the MHC model, we can simply divide a chosen function by $\sqrt{\widetilde{c}_{R}\widetilde{c}_{O}}$ and scale the magnitude accordingly. To connect to our previous work, we choose to compare the (symmetric) Butler-Volmer form based on species activities, Eq. 35 with Eq. 37. Thus, for the MHC model, we choose

[TABLE]

with $n=1$ and $\alpha=0.5$, which very nearly matches the concentration and activity dependence of the exchange current densities of the two models, leaving the primary difference in the dependence on the activation overpotential, $\eta$. We use \gamma_{\ddagger}=1/\mathopen{}\mathclose{{}\left(1-\widetilde{c}_{R}}\right) following previous work on LFP Bai et al. (2011); Cogswell and Bazant (2012) and supported by recent experiments Lim et al. (2016). This gives an exchange current density with broken symmetry around half filling, as depicted in Figure 11 (b). We choose a reaction reorganization energy for the MHC model $\lambda=18k_{\mathrm{B}}T$ as in Figure 10.

In order to minimize the effect of solid phase transport losses, we modify the phase separating particles from Section V.1 to have $D_{0}=1\times 10^{-14}\ \mathrm{m}^{2}/\mathrm{s}$. We simulate a half-cell using the concentrated Stefan-Maxwell electrolyte model and with a lithium foil counter-electrode with fast kinetics. Although fast lithium foil kinetics may be a poor assumption because of the small relative surface area, under constant current conditions, it would only add a current-dependent, constant additional voltage drop to both simulated systems. The simulated cell has a separator of length $20\ \mu\mathrm{m}$ with porosity $0.8$ and a thin cathode of length $25\ \mu\mathrm{m}$ with porosity $0.2$ to minimize electrolyte transport losses and focus on the effect of the reaction kinetics. The cathode loading percent, as above, is set to $0.7$, and we simulate uniform particles of radius $1\ \mu\mathrm{m}$ with $k_{0}=1\ \mathrm{A}/\mathrm{m}^{2}$ and $k_{M}$ chosen such that the magnitudes of the exchange current densities match. We adjust the value of the film resistance from that used in Figure 10 to $R_{\mathrm{film}}=0.02\ \Omega\ \mathrm{m}^{2}$ to achieve a closer match between the MHC and Butler-Volmer model with the film while using the concentration-dependent exchange current densities neglected in Figure 10. This value is higher than is commonly used and closer to that used in the introduction of the film resistance to battery modeling by Doyle et al. Doyle et al. (1996).

We expose the cell to the current profile and corresponding state of charge profile in Figure 12 (b), and the output voltage profiles are shown in Figure 12 (a). Although the cell experiences moderate electrolyte polarization at the highest currents, those losses are similar for both models and do not explain the differences in the overall output voltage. The initial current pulse at $2\ \mathrm{C}$ brings all the electrode particles to a state of charge within the miscibility gap, so the equilibrium voltage is given by the value of $-\mu^{\Theta}/e=2\ \mathrm{V}$ which we see the cell relax to after the pulse. We also can see the overlap of the BV and MHC profiles in the initial region of this pulse, where the surface concentrations are changing most substantially, confirming the match of both the magnitude and concentration dependence of the two reaction models in this lower-overpotential range when the two reaction models should overlap. However, once the particles begin to phase separate at approximately $7.5\ \mathrm{min}$, the surface concentrations rise sharply, and the exchange current density decreases (Figure 11 (b)). This causes increased reaction resistance in both models and an associated increase in required driving force, but the MHC model requires a larger increase in driving force, causing the initial departure of the models near $7.5\ \mathrm{min}$. When the higher current pulses begin, we see further departure of the two reaction models with the MHC model showing substantially higher reaction losses in Figure 12 (a) than the Butler-Volmer model. The Butler-Volmer model with the film resistance leads to behavior between the two, predicting more resistance than the other models at low rates and intermediate resistance at higher rates, indicating that the model results depart enough to clearly distinguish in these conditions. Still, the effect of the film resistance and MHC kinetics are qualitatively similar, suggesting that some of the system behavior which has been attributed to reaction film resistance may instead be caused by fundamental departures from Butler-Volmer reaction kinetics.

In the limit of $e\eta\ll\lambda$, the Butler-Volmer and MHC models predict identical results. The value of $\lambda$ we used here is between a value calculated for LFP Bai and Bazant (2014) and values commonly found in calculations or fits to data for other systems Chidsey (1991); Henstridge et al. (2012). Using this intermediate value and an MHC prefactor adjusted to give a more similar concentration dependence to Butler-Volmer expressions, we can see that the two models give noticeably different predictions on the simulated macroscopic system output. This further highlights the importance of continuing the study of the two models in practical battery models, and the MPET software can help facilitate this.

VI Summary and Conclusions

Volume averaged porous electrode simulations provide insightful and industrially relevant ways to characterize battery behavior. Despite shortcomings associated with information loss from the volume averaging process, the simplicity of the approach and associated speed of model development, implementation, and simulation run time motivate its continued use. In this work, we have presented an open-source software package called MPET (Multiphase Porous Electrode Theory), which builds on foundations laid by John Newman and many others by describing the active materials with variational nonequilibrium thermodynamics Bazant (2013, 2017) applied to porous electrodes Ferguson and Bazant (2012, 2014). Despite the prevalence of this modeling approach, few open source options are available for simulating the model, particularly ones that are easy to modify with new thermodynamic models based on the powerful phase field formalism Chen (2002), adapted for electrochemical systems Bazant (2013). With MPET, we aimed to address this gap by providing a software platform implementing nonequilibrium thermodynamics of porous electrodes with an open source code, written with a modular design to encourage use, modifications, and improvements.

Through a variety of examples, we have demonstrated some of key features of the MPET software. First, we compared solid solution and phase field approaches to modeling active materials and demonstrated that the models can give similar macroscopic outputs under some situations, but deviate at the particle scale. This leads to different predictions when the surface concentrations strongly affect reaction behavior. Second, we reexamined the comparison of Stefan-Maxwell concentrated solution theory and dilute solution models of electrolyte transport in the context of electrodes made of simple phase separating particles. Under strong electrolyte limitation, the difference in the model predictions is similar to that found in electrodes with solid solution particles. Third, to highlight the ability of the software to describe unique and distinct solid models, we implemented a full cell simulation using models of graphite and iron phosphate presented in previous works. Finally, we considered the effect of changing the reaction model from the typically used Butler-Volmer to Marcus-Hush-Chidsey (MHC) reaction kinetics, based on microscopic theories of electron transfer. We demonstrated that, for some reasonably typical parameter values, the MHC reaction kinetics can look similar to Butler-Volmer reactions with a film resistance and can lead to substantial differences from the Butler-Volmer model in predicted battery behavior at high rates.

Natural extensions of this work involve implementing some of the features that other software options have and which researchers have found helpful in explaining data or better describing real systems. For example, thermal effects can substantially affect cell behavior Pals and Newman (1995); Gu and Wang (2000); Thomas and Newman (2003); Lai and Ciucci (2011); Prada et al. (2012); Northrop et al. (2015); Latz and Zausch (2015); Torchio et al. (2016) and exploration of their coupling with Marcus kinetics would be interesting. More complete thermal descriptions rely on temperature profiles over more than an individual cell layer Forgez et al. (2010); Pesaran (2002), so isothermal but non-constant temperature dependence would be a starting point. Addition of material models properly coupling the stresses to the concentration profile would also be an opportunity to study the effects of electro-mechanical models with phase separation Cogswell and Bazant (2012); Di Leo et al. (2014) in porous electrodes. Other capabilities, such as simulating electrochemical impedance outputs or others of the many additions that have been made to the original model implemented by Doyle et al. Doyle et al. (1993) could also be added. For capacitive energy storage and desalination, the electrolyte model could also be extended to allow for diffuse charge in the electrode/electrolyte interface Biesheuvel et al. (2011, 2012) or the double layers of charged porous separators Dydek et al. (2011); Dydek and Bazant (2013); Schmuck and Bazant (2015); Yaroshchuk (2012), which also activates additional mechanisms for ion transport by surface conduction and electro-osmotic flow Dydek et al. (2011); Deng et al. (2013); Mirzadeh et al. (2014), which are neglected in traditional battery models.

In summary, MPET provides some new capabilities for battery simulation focused on recent developments in the modeling of active materials based on nonequilibrium thermodynamics. It can also serve as a starting point for other researchers beginning to study in the area to make their own modifications and investigations. By highlighting its capabilities, we have shown the value of a flexible simulation package to expand on the existing porous electrode theory and begun to examine the impact of those developments.

VII Acknowledgments

The research was supported by the Samsung-MIT Program for Materials Design in Energy Applications, and in part by the D3BATT program of the Toyota Research Institute. We thank E. Khoo and K. Conforti for proofreading the manuscript.