Multiphase Porous Electrode Theory
Raymond B. Smith Martin Z. Bazant




Department of Chemical Engineering, Massachusetts
Institute of Technology, Cambridge, Massachusetts 02139, USA


[Corresponding author: ]bazant@mit.edu
Department of Chemical Engineering, Massachusetts
Institute of Technology, Cambridge, Massachusetts 02139, USA
Department of Mathematics, Massachusetts Institute of
Technology, Cambridge, Massachusetts 02139, USA



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.


Multiphase Porous Electrode Theory
    Martin Z. Bazant
    December 30, 2023
==================================





§ INTRODUCTION

Lithium-based batteries have growing importance in global society <cit.> as a result of increased prevalence of portable electronic devices <cit.>, and their enabling role in the transition toward renewable energy sources <cit.>.
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 <ref>.
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 <cit.>.
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 <cit.> as a result of inhomogeneities over a length scale smaller than the electrode but many times the size of a primary active particle <cit.>.
The active materials themselves also often have highly non-trivial behavior including poor connectivity with each other and the conductive matrix <cit.> as well as complex material properties leading to deformations and accompanying stresses <cit.> and phase separation <cit.> during the intercalation process.
Coupling these behaviors in multi-particle electrode environments leads to further complexities <cit.>.

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 <cit.>.
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 <cit.>.
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. <cit.> 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 <cit.>.
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 <cit.> or shrinking annuli <cit.>, 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 <cit.> and porous electrodes <cit.>.
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 <cit.>.
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 <cit.>, as well as degradation related to mechanical stresses <cit.> and side reactions that depend on the local surface concentration profile <cit.>.

Regardless of the thermodynamic model, volume averaged simulations using porous electrode theory are carried out in a number of ways.
Newman's  code uses a finite difference method via the BAND subroutine, and it is freely available <cit.> 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 <cit.> have also been used to implement versions of porous electrode theory <cit.>, 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 <cit.>, 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 <cit.>.
We take a similar approach in this work.
More comprehensive reviews of commonly used simulation approaches can be found in refs. <cit.>.

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 <cit.>, 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 <ref>, we begin by presenting the full mathematical framework of MPET, based on the original formulation of Ferguson and Bazant <cit.>, with several modifications.
First, we incorporate the standard description of transport in concentrated electrolytes based on Stefan-Maxwell coupled fluxes and chemical diffusivities <cit.>.
Second, we capture the continued development by our group of phase field models for the active solid materials <cit.>, which have increasingly been validated by direct experimental observations of phase separation dynamics <cit.>.
Third, we provide alternatives to the empirical Butler-Volmer model of Faradaic reaction kinetics <cit.>, based on the quantum mechanical theory of electron transfer pioneered by Marcus <cit.> and extended here for concentrated solutions <cit.>, motivated by recent battery experiments <cit.>.
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 <ref>, the equations are made dimensionless, and numerical methods to solve them are presented in Section <ref>, along with the overall software structure.
In Section <ref>, 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 <ref>.



§ 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.



 §.§ Electrode Scale Equations



  §.§.§ 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 <cit.>.
We consider electrolytes of salts defined by ∑_iν_iM_i^z_i, with species M_i having valence z_i, and ν_i the number of ions M_i in solution from dissolving one molecule of the neutral salt.
For example, for CuCl_2, ν_+ = 1, M_+ = Cu, z_+ = 2, ν_- = 2, M_- = Cl, and z_- = -1.
With electrolyte species flux, 𝐅_ℓ,i defined per area of porous medium, conservation of species requires

    ∂(  c_ℓ,i)/∂ t = -𝐅_ℓ,i + R_V,i

where c_ℓ,i is the concentration of species i in the electrolyte,  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 <cit.> or conducting porous electrodes <cit.> 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 <cit.>, or hybrid pseudo-capacitors <cit.>.
In batteries, however, such effects are usually neglected <cit.>, 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 <cit.>.
Charge conservation and quasi-neutrality together require

    ∂( ρ_e )/∂ t≈ 0 = -𝐢_ℓ
        + ∑_i z_ieR_V,i,

where e is the elementary charge, the charge density ρ_e = ∑_i z_iec_ℓ,i, and 𝐢_ℓ is the current density in the electrolyte, related to a sum of ionic fluxes,

    𝐢_ℓ   = ∑_i z_ie𝐅_ℓ,i.

We will relate fluxes to both concentrations and electrostatic potentials in the electrolyte, so with constitutive flux relationships, Eqs. <ref> and <ref> 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)

    z_nc_ℓ,n = -∑_i nz_ic_ℓ,i,

which allows us to neglect Eq. <ref> 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_ℓ = c_ℓ,+/ν_+ = c_ℓ,-/ν_- and instead simulate one of

    ∂(  c_ℓ)/∂ t =
    
        1/ν_i(-𝐅_ℓ,i + R_V,i).

and

    0 = -𝐢_ℓ + ( z_+eR_V,+ + z_-eR_V,-).

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, 𝐧·𝐅_ℓ,i = 0 and 𝐧·𝐢_ℓ = 0.
However, if a foil (e.g. Li metal electrode) is used for one electrode, the boundary condition at that side is replaced by 𝐧·𝐢_ℓ = i_cell, where i_cell is the macroscopic current density of the cell, Eq. <ref>. 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, 𝐧·𝐅_ℓ,- = 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, 𝐅_ℓ,i to gradients in electrochemical potentials, μ_ℓ,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 <cit.>,

    μ_ℓ,i = k_BTln( a_ℓ,i) + μ_ℓ,i^Θ + z_ieϕ_ℓ = k_BTln( c_ℓ,i) + μ_ℓ,i^ex

where k_B is the Boltzmann constant, T is the absolute temperature, a_ℓ,i is the activity of species i, μ_ℓ,i^Θ is its reference chemical potential, and ϕ_ℓ is the (Galvani) electrostatic potential.
The final expression serves as the definition of the excess chemical potential, μ_ℓ,i^ex, and c_ℓ,i is the concentration, c_ℓ,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 <cit.>.


An alternative approach common in battery modeling <cit.> 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. <ref>, 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 ϕ_ℓ^r, and we note that it is actually measuring (a combination of) full electrochemical potentials of ions in solution.
For the Li/Li^+ example,

    eϕ_ℓ^r = μ_ℓ,Li^+ +  = k_BTln( a_ℓ,Li^+) + eϕ_ℓ + μ_ℓ,Li^+^Θ + .

Quasi-electrostatic potentials, ϕ_ℓ^q, can also be used in electrolyte models <cit.>, and they are defined by the excess chemical potential of a particular ion in solution, n,

    z_neϕ_ℓ^q = μ_ℓ,n^ex = z_neϕ_ℓ + k_BTln( γ_ℓ,n) + μ_ℓ,n^Θ

where γ_ℓ,i = a_ℓ,i/c_ℓ,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. ∑_i z_ic_ℓ,i = 0 everywhere, or (2) by solving the Poisson equation, ( εϕ_ℓ) = -ρ_e with permittivity ε, which could depend on concentration or electric field <cit.>, or capture non-local ion-ion correlations as a differential operator <cit.>.
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,

    𝐅_ℓ,i = -M_ℓ,ic_ℓ,iμ_ℓ,i

where M_ℓ,i is the species' mobility.
Using the Galvani potential form in Eq. <ref>,

    𝐅_ℓ,i   = -D_ℓ,i/T(c_ℓ,i(Tln( a_ℓ,i)) + z_ic_ℓ,iϕ_ℓ),

where we have used the Einstein relation between the diffusivity in free solution, D_ℓ,i and mobility, D_ℓ,i = M_ℓ,ik_BT, non-dimensionalized the potential by the thermal voltage at some reference temperature, ϕ_ℓ = eϕ_ℓ/k_BT_ref, and defined a non-dimensional T = T/T_ref.
With D_ℓ,chem,i = D_ℓ,i( 1 + ∂lnγ_ℓ,i/∂lnc_ℓ,i) and uniform temperature,



    𝐅_ℓ,i = -(D_ℓ,chem,i c_ℓ,i - D_ℓ,i/Tz_ic_ℓ,iϕ_ℓ).

In a porous medium we use effective transport properties, which are adjusted by the tortuosity of the electrolyte phase, τ.
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, .

    𝐅_ℓ,i   = -/τ(D_ℓ,chem,i c_ℓ,i + D_ℓ,i/Tz_ic_ℓ,iϕ_ℓ).

The tortuosity is often described as a function of the porosity, the volume fraction of the electrolyte phase, , commonly by employing the Bruggeman relation τ = ^a <cit.>.
The value of a is often set to -0.5 but can be adjusted <cit.> to account for experimentally <cit.> or theoretically <cit.> observed departures from the original derivation.
Eq. <ref> with Eqs. <ref>, <ref> and <ref> 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,

    𝐅_ℓ,i = ∑_jU_ijμ_ℓ,j.

where U_ij are the direct (i=j) and indirect (i  j) transport coefficients <cit.>.
This formulation has been used to describe relatively concentrated electrolytes and is the most commonly used model in battery simulation <cit.>.
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, 𝐯_i = 𝐅_ℓ,i/c_ℓ,i, <cit.>

    c_ℓ,iμ_ℓ,i = ∑_jK_ij( 𝐯_j - 𝐯_i )
        = k_BT∑_jc_ℓ,ic_ℓ,j/c_T𝒟_ij( 𝐯_j - 𝐯_i )

with c_T = ∑_ic_ℓ,i and K_ij = K_ji.

For a binary electrolyte in a porous medium with cations, anions, and solvent (denoted as species 0), assuming uniform temperature and that the solvent concentration varies only negligibly with salt concentration, and again neglecting convection,

    𝐅_ℓ,+   = -ν_+/τD_ℓ c_ℓ + t_+^0𝐢_ℓ/z_+e
    𝐅_ℓ,-   =  -ν_-/τD_ℓ c_ℓ + t_-^0𝐢_ℓ/z_-e

where, defining γ_ℓ,±^ν = γ_ℓ,+^ν_+γ_ℓ,-^ν_- with ν=ν_+ + ν_-,

    D_ℓ   = 𝒟c_T/c_ℓ,0( 1 + ∂lnγ_ℓ,±/∂lnc_ℓ),
        
       = _0+_0-( z_+ - z_- )/z_+_0+ - z_-_0-,
        
    
        t_+^0    = 1 - t_-^0 = z_+_0+/z_+_0+ - z_-_0-.

The current density is given by

    𝐢_ℓ = -σ_ℓ/Tτ(ϕ_ℓ^r
        + ν k_BT/e( s_+/nν_+ + t_+^0/z_+ν_+ - s_0c_ℓ/nc_ℓ,0)
        ( 1 + ∂lnγ_ℓ,±/∂lnc_ℓ)lnc_ℓ)

with

    1/σ_ℓ   =
        -( k_BT_ref/c_Tz_+z_-e^2)
        ( 1/_+-
        + c_ℓ,0t_-^0/c_ℓ,+_0-).

The values of s_i are specified by the choice for the reference electrode with reaction

    s_-M_-^z_- + s_+M_+^z_+ + s_0M_0⇌ ne^-.

For lithium-ion batteries, the typical choice for the reference electrode defining ϕ_ℓ^r is Li/Li^+, so s_+ = -1, and s_- = s_0 = 0.
Thus, Eqs. <ref> and <ref> with Eqs. <ref> and <ref> define the electrolyte model when using the Stefan-Maxwell concentrated solution theory with a binary electrolyte.



  §.§.§ 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,

    0 = -𝐢_s - ∑_iz_ieR_V,i

or

    0 = -𝐢_s - (z_+eR_V,+ + z_-eR_V,-)

for a binary electrolyte in which cations and/or anions may undergo electrochemical reactions, and where 𝐢_s is the current density in the solid phase.
The sign difference compared to Eq. <ref> 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,

    𝐢_s = -( 1-)/τσ_sϕ_s

where σ_s is the conductivity of the electronically conductive matrix, ϕ_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, 𝐧·𝐢_s = 0, and the potential at the current collector is specified by the operating voltage on the system in the macroscopic equations, ϕ_c or ϕ_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 <cit.>.
We demonstrate use of the parallel case with a distribution of contact resistances in ref. <cit.>.
In the second case of series wiring we implement a simplified version of that developed by Stephenson et al. <cit.> by imposing a finite conductance between particles in series, indexed by k, within a simulation volume, j,

    G_j,k(ϕ_j,k - ϕ_j,k+1) = I_j,k

where G_j,k is the conductance and I_j,k is the current between particle k and k+1.
From charge conservation,

    I_j,k - I_j,k+1 = ∫_S_k+1^j_j,k+1 A,

where j_j,k+1 is the intercalation rate into particle k+1.



 §.§ 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 <ref>), 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,

    G
        = ∫_V_p^ g  V + ∫_A_S^γ_S A

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 γ_S is the surface energy.
Typically, we will separate the free energy density into homogeneous, g_h and non-homogeneous, g_nh, contributions,

    g = g_h + g_nh + …

where the remaining terms could describe the stress state of the system <cit.> or other energetic contributions <cit.>.
Following van der Waals <cit.> and Cahn and Hilliard <cit.>, we use a simple gradient penalty term to describe the non-homogeneous free energy,

    g_nh = 1/21/c_s,ref^2 c_i·κ c_i

where κ is a gradient penalty tensor (assumed to be isotropic here such that κ = κ1 and 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,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,

    μ_i = δ G/δ c_i = ∂ g/∂ c_i - ∂ g/∂ c_i.




  §.§.§ Electrochemical Reactions


The electrochemical reaction can be described by a number of different models, such as the empirical Butler-Volmer equation <cit.> or quantum-mechanical models based on Marcus kinetics <cit.>, which must be consistently generalized for concentrated solutions in nonequilibrium thermodynamics <cit.>.
We describe here electrochemical reactions of the form

    S_1 = ∑_i s_i,OO_i^z_i,O + ne^- →∑_j s_j,RR_j^z_j,R = S_2,

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

    neη = μ_R - (μ_O + nμ_e) = Δμ_rxn = Δ G_rxn,

where μ_R = ∑_js_j,Rμ_j is the electrochemical potential of the reduced state, μ_O = ∑_i s_i,Oμ_i is the electrochemical potential of the oxidized state, and μ_e is the electrochemical potential of the electrons, which we relate here to the electric potential measured in the conductive matrix, μ_e = -eϕ_s.
Δμ_rxn and Δ G_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, μ_O = eϕ_ℓ^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 <cit.> to the reaction rate model, and have been recently reviewed <cit.>.
Frumkin-corrected Butler-Volmer reaction models have been applied to various electrochemical techniques including steady constant current <cit.>, voltage steps <cit.>, current steps <cit.>, and linear sweep voltammetry <cit.>, as well as nano <cit.> and porous <cit.> electrodes, with clear indication of departure from models neglecting double layers, especially at low salt concentrations with thick double layers (“Gouy-Chapman limit” <cit.>).
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” <cit.>).
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

    i = i_0( exp( -α eη_eff/k_BT)
        -exp( ( 1-α)eη_eff/k_BT))

where α 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 <cit.> or a function of species concentrations or activities.
Introduced in ref. <cit.>, the effective overpotential, η_eff, accounts for any film resistance, R_film, by

    η_eff = η + iR_film.

Bazant and co-workers proposed a form for i_0 based on reacting species' activities and a transition state activity coefficient, γ_ <cit.>, 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 <cit.>,

    i_0 = k_0ne( a_Oa_e^n )^1-αa_R^α/γ_

where k_0 is a rate constant, a_O = ∏_ia_i^s_i,O, a_e is the activity of the electrons (taken to be unity here), and a_R = ∏_ja_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 LiFePO_4, Bai et al. originally proposed γ_ = ( 1 - c/c_max)^-1 where c_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_Li^+ = c_ℓ for both electrolyte models.

It is also common to use an exchange current density form based solely on species concentrations <cit.>,

    i_0 = k_0c_ℓ^ 1-αc_i^ α( 1 - ξ_ic_i)^α,

where c_ℓ = c_ℓ/c_ℓ,ref, c_i = c_i/c_s,ref, and ξ_i = c_s,ref/c_i,max, and the ref subscripts indicate some suitable concentration scale. For the case of lithium insertion, we will choose c_s,ref = c_i,max, so ξ_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 <cit.>.
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 <cit.>, 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),

    i_red   = ek^0c_O∫_-∞^∞ρ(z)n_e(z)W_red z
        
    
        i_ox   = ek^0c_R∫_-∞^∞ρ(z)( 1-n_e(z) )W_ox z

where k^0 is some rate prefactor, c_O = ∏_ic_i^s_i,O and c_R = ∏_jc_j^s_j,R, z represents the energy level of electrons in the electrode, ρ(z) is the density of states, and n_e(z) is the Fermi function.
W_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,

    W_red   = k_wexp( -w_O/k_BT)
        exp( -( λ + eη_f )^2/4λ k_BT)
        
    
        W_ox   = k_wexp( -w_R/k_BT)
        exp( -( λ - eη_f )^2/4λ k_BT),

where k_w is a prefactor with explicit dependence on various factors, including the overlap integrals for the wave functions of the various reaction elements <cit.>.
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.
λ 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, Δ G_rxn^ex,

    eη_f = Δ G_rxn^ex = eη + k_BTln( c_O/c_R)

where we have assumed single-electron transfer, as supported by Marcus theory for elementary reaction events <cit.>.
The tildes indicate non-dimensional concentrations c_O = ∏_ic_i^s_i,O and c_R = ∏_jc_j^s_j,R, with each scaled to its reference concentration (c_ℓ,ref for ions in the electrolyte and c_s,ref for intercalated species in the active material).
Thus, η_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 ±λ, and experimental validation of this so-called “inverted” region of decreasing reaction rate at increasing driving force  (i.e. negative differential reaction resistance <cit.>) paved the way for the Nobel Prize of Marcus <cit.>.

Assuming w_R/O are independent of the electron energy level (not necessarily true <cit.>), 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 <cit.>,

    i_red/ox = k_Mc_O/Rexp( -w_O/R/k_BT)
        ∫_-∞^∞exp( -( z - λ∓ eη_f )^2/4λ k_BT) z/1 + exp( z/k_BT ),

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 <cit.>, 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 λ for the initial and final states, referred to as “asymmetric Marcus Hush” (AMH) kinetics <cit.>.
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. <cit.>.

We make one final assumption following ref. <cit.> 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 <cit.>, 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,

    w_O = w_R = μ_M,^ex

or

    exp( -w_R/O/k_BT) ∝1/γ_M,.

Thus, we arrive at our expressions for MHC kinetics,

    i = i_M( c_Ok_red - c_Rk_ox)

with

    k_red/ox = ∫_-∞^∞exp( -( z - λ∓ eη_f )^2/4λ k_BT) z/1 + exp( z/k_BT )

and

    i_M = k_M/γ_M,.

Although the derivation of the Butler-Volmer equation <cit.> 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 γ_M,.
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. <cit.>, we replace Eq. <ref> with

    k_red/ox≈√(πλ)/1 + exp( ±η_f )( λ - √(1 + √(λ) + η_f^2)/2√(λ))

where tildes indicate scaling by the thermal energy or voltage, k_BT or k_BT/e, (z) = 1 - (z) is the complementary error function, and reduction/oxidation corresponds to the top/bottom sign.



  §.§.§ 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 <cit.>.



  
Homogeneous
When transport within the solid particles is fast, it can be computationally beneficial to approximate the particles with an average concentration, c_s <cit.>.
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, j_p,

    ∂c_i/∂ t = A_p/V_pj_p,i.




  
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 LiFePO_4 <cit.>.
Then, neglecting transport, the local rate of change of the concentration at particle location 𝐫 is

    ∂ c_i/∂ t = A_p/V_pj_p,i(𝐫).




  
Cahn-Hilliard Reaction
In this model, transport of the intercalant within the particles is described by species conservation,

    ∂ c_i/∂ t = -𝐅_i.

The flux can be modeled using linear irreversible thermodynamics <cit.>, which postulates that fluxes arise from gradients in the diffusional chemical potential,

    𝐅_i = -D_i/Tc_iμ_i,

where we have again used the Einstein relation.
The diffusional chemical potential is scaled to the thermal energy at some reference temperature, k_BT_ref, and can be calculated from Eq. <ref>.
Following Bazant <cit.>, D_i can be rewritten in terms of tracer diffusivity in the dilute limit, D_0,i, the activity coefficient of species i, γ_i = a_i/c_i, and the activity coefficient of the diffusion transition state, γ_,i^d,

    D_i = D_0,iγ_i/γ_,i^d.

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, γ_i/γ_,i^d = ( 1-c_i/c_i,max), giving

    𝐅_i = -D_0,i/Tc_i( 1-c_i/c_i,max)μ_i,

although other effects could be accounted for in γ_,i^d including stresses in the transition state <cit.>.
At the particle surface, the flux is given by the electrochemical reaction, 𝐧·𝐅_i = -j_p,i where 𝐧 is a unit normal vector pointing from the active material to the electrolyte.
The natural boundary condition <cit.> imposes a constraint on the concentration at the surface as a function of the surface energy, γ_S, 𝐧·∂ g/∂c_i = 𝐧·κc_i = ∂γ_S/∂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. <ref> can be rewritten as

    𝐅_i = -D_i/Tc_i∂μ_i/∂ c_i c_i
        = -D_chem,i/T c_i

where D_chem,i = D_ic_i∂μ_i/∂ c_i.
Here, it is sufficient to prescribe the flux at the surface, as in the Cahn-Hilliard reaction model, 𝐧·𝐅_i = -j_p,i.



 §.§ 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 <cit.>.
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,

    ∂c_p,i/∂ t = 1/V_p∫_A_p^j_p,i A.

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

    R_V,i = -( 1-)P_L∑_pV_p/V_u∂c_p,i/∂ t

where P_L is the loading percent of active material in the solid phase, V_u = ∑_pV_p, and p indicates properties of particle p, and the summation is over particles at the location where R_V,i is evaluated.



 §.§ Macroscopic Equations

The overall current density per electrode area, i_cell, is defined as the integral of the net charge consumed by the reactions in the electrodes per unit cross-sectional area,

    i_cell = ∑_i∫_L_a^z_ieR_V,i L
        = -∑_i∫_L_c^z_ieR_V,i L

where L_a and L_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_max,k in electrode k by the capacity per area of the limiting electrode,

    Q_A = min_k∈{ a,c }{ eL_k( 1-_k )P_L,kc_max,k}.

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

    i_cell,C = i_cell/Q_Aτ_hr

where τ_hr is a conversion factor to ensure units of 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,

    u_k,i = 1/L_k∫_L_k^∑_pV_p/V_uξ_ic_p,i L.

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, ξ_i = c_s,ref/c_i,max, and c_p,i = c_p,i/c_s,ref, so the product is the average filling fraction of particle p.

The overall cell voltage, Δϕ_cell = ϕ_c - ϕ_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 ϕ_c = 0 in the simulations.
We account for series resistance by defining an applied voltage, Δϕ_appl, by

    Δϕ_cell = Δϕ_appl - i_cellR_ser,

where R_ser is the area specific resistance of the cell.



§ NON-DIMENSIONAL EQUATIONS


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




 §.§ Electrode Scale Equations




  §.§.§ Electrolyte Model

Defining a reference electrolyte concentration, c_ℓ,ref (e.g. c_ℓ,ref/N_A = 1 M, where N_A is the Avogadro constant),

    ∂( c_ℓ,i)/∂t   =
        1/ν_i(-𝐅_ℓ,i + R_V,i)
        
    
        0    = -𝐢_ℓ + ∑_iz_iR_V,i
    𝐢_ℓ   = ∑_iz_i𝐅_ℓ,i

with = L_ref and tildes indicating non-dimensionalization by the following scales

    F_ℓ,ref   = c_ℓ,refL_ℓ,ref/t_ref
    
        R_V,ref   = c_ℓ,ref/t_ref
    
        i_ℓ,ref   = eF_ℓ,ref.

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

    𝐅_ℓ,i   = -/τ(D_ℓ,chem,ic_ℓ,i + D_ℓ,i/Tz_ic_ℓ,iϕ_ℓ)

with diffusive scale as that chosen to define t_ref above.
The binary concentrated solution theory can be non-dimensionalized with

    𝐅_ℓ,+   = -ν_+/τD_ℓc_ℓ + t_+^0𝐢_ℓ/z_+
    𝐅_ℓ,-   =  -ν_-/τD_ℓc_ℓ + t_-^0𝐢_ℓ/z_-
    
        𝐢_ℓ   = -σ_ℓ/Tτ(ϕ_ℓ^r
        + νT( s_+/nν_+ + t_+^0/z_+ν_+ - s_0c_ℓ/nc_ℓ,0)
        ( 1 + ∂lnγ_±/∂lnc_ℓ)lnc_ℓ)

with conductivity scale

    σ_ref = e^2c_ℓ,refD_ℓ,ref/k_BT_ref.




  §.§.§ Solid phase electronic model

The solid phase current density and conductivity are scaled to i_ref and σ_ref respectively such that

    0 = -𝐢_s - ( z_+R_V,+ + z_-R_V,-)

and

    𝐢_s = -( 1-)/τσϕ_s.









 §.§ 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,

    F_s,ref   = L_s,refc_s,ref/t_ref

where L_s,ref and c_s,ref are relevant scales to the particle. We will use L_s,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_max, we choose the solid reference concentration to be the maximum filling for lithium insertion cells, c_s,ref = c_max.



  §.§.§ Electrochemical Reactions

The various forms of the overpotential are all scaled to the thermal voltage, k_BT_ref/e, and the rate prefactors and current density are both scaled to

    i_s,ref   = eF_s,ref.

The film resistance is scaled to R_film,ref = k_BT_ref/ei_s,ref, so

    η_eff = η + iR_film

and the Butler-Volmer expression is

    i = i_0( exp( -αη_eff/T)
        - exp( ( 1-α)η_eff/T))

while the MHC expression is

    i = i_M( c_Ok_red - c_Rk_ox).




  §.§.§ 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,ref, such that, for homogeneous particles,

    ∂c_i/∂t = δ_p,Lj_p,i

where c_i = c_i/c_s,ref, and

    δ_p,L = A_pL_p/V_p.

For Allen-Cahn reaction particles,

    ∂c_i/∂t = δ_p,Lj_p,i( 𝐫 )

where 𝐫 is non-dimensionalized by L_p.
For Cahn-Hilliard reaction particles,

    ∂c_i/∂t = -𝐅_i.

The non-dimensional flux is given by

    𝐅_s,i = -D_i/Tc_iμ_i

where the diffusivity is scaled to

    D_s,ref = L_s,ref^2/t_ref.

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

    𝐅_i = -D_0,i/Tc_i( 1-ξ_ic_i )μ_i.

The natural boundary condition is scaled using

    κ_ref   = k_BT_refL_s,ref^2c_s,ref
    γ_S,ref   = κ_ref/L_s,ref

such that

    𝐧·c_i = 1/κ∂γ_S/∂c_s

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

    𝐧·𝐅_i = -j_p,i.

For solid solution particles,

    𝐅_i = -D_s,chem,i/Tc_s,i

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



 §.§ Coupling Equations

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

    ∂c_p,i/∂t
        = δ_p,L∫_A_p^j_p,iA

and

    R_V,i = -β( 1-)P_L∑_pV_p/V_u∂c_p,i/∂t

with β = c_s,ref/c_ℓ,ref.



 §.§ Macroscopic Equations

Keeping the same scales as in the electrode model,

    i_cell = ∑_i∫_L_a^z_iR_V,iL
        = -∑_i∫_L_c^z_iR_V,iL

and

    u_k,i = 1/L_k∫_L_k^∑_pV_pξ_ic_p,iL.

where k indicates either the anode or cathode, L_k = L_k/L_ℓ,ref, and V_p = V_p/V_u.
The series resistance can be scaled by the electrode model scales, R_ser,ref = k_BT_ref/ei_ref, and the potentials can all be scaled by the thermal voltage at the reference temperature, k_BT_ref/e, so

    Δϕ_cell = Δϕ_appl - i_cellR_ser.




§ 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 <cit.>.
We use the IDAS time stepper of the SUNDIALS integration suite <cit.> to solve the resulting DAE's with a backward differentiation formula approach.
The model and simulation are defined within the DAE Tools software package <cit.>, which provides a modeling language environment within Python similar to that of more specialized modeling languages like gProms <cit.> or Modelica <cit.>.
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 <cit.> 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.



 §.§ 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 <ref>.
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 <ref>, 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 <cit.> or hierarchical structures <cit.>.
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 <cit.>.
Other variations such as 2D particles of arbitrary geometry discretized using the finite volume <cit.> or finite element method would also be straightforward, especially with the DAE Tools interface to the deal.ii finite element library <cit.>.
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.



 §.§ 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. <cit.>).
With minor changes to the code, completely arbitrary functions are possible.
Currents are imposed relative to a 1 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 <ref>.
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 <ref>.
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 <cit.>, 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 <ref>, 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 <cit.>.



 §.§ Numerical Methods

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


  §.§.§ Electrolyte

The electrolyte equations are discretized using a typical 1D finite volume scheme with cell centers <cit.>.
They are non-dimensionalized using the same scales as in the electrode scale equations, Section <ref>.
For example, for species conservation in cell j with width Δx and neglecting subscripts for electrolyte/species identifiers for clarity,

    ∂c_j/∂t = F_j-1/2 - F_j+1/2/Δx

where the j subscript represents an average quantity within the discretized volume, and the F_j± 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 <cit.>.



  §.§.§ Active Material

The active material can be defined as spherical, cylindrical, or a rectangular grid approximation of platelet-like particles common in LiFePO_4 <cit.>.
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 <cit.>.
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. <cit.>.
They are non-dimensionalized using the same scales as in the single particle scale equations, Section <ref>.
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

    𝐌𝐕∂𝐜/∂t = 𝐛,

where 𝐜 is a vector of concentrations at positions going from the center of the particle to the surface of the particle with N sub-volumes, 𝐌 is a tridiagonal matrix given by

    𝐌   = [ 3/4 1/8   0   0   ⋯   0   0   0   0; 1/4 3/4 1/8   0   ⋯   0   0   0   0;   0 1/8 3/4 1/8   ⋯   0   0   0   0;   ⋮   ⋮   ⋮   ⋮   ⋱   ⋮   ⋮   ⋮   ⋮;   0   0   0   0   ⋯ 1/8 3/4 1/8   0;   0   0   0   0   ⋯   0 1/8 3/4 1/4;   0   0   0   0   ⋯   0   0 1/8 3/4;     ],

and 𝐕 is a diagonal matrix defined by the volumes of the shells scaled to L_p^3.
Indexing the volumes with j=1… N, and with a uniformly spaced radial vector, r_j, with N points going from the particle center to its surface, then for a sphere 𝐕 has components defined by

    V_jj^sphere = 
            4π( Δr^3/24),     j = 1 
    
            4π( r_j^2Δr + Δr^3/12),     j = 2… N-1 
    
            4π( r_j^3/3 - ( r_j - Δr/2 )^3/3),     j = N 
    ,

and for a cylinder,

    V_jj^cylinder = πh_p Δr^2/4,     j = 1 
    
            2πh_pr_jΔr,     j = 2… N-1 
    πh_p( r_jΔr - Δr^2/4),     j = N 
    .

The right hand side, 𝐛, is defined in relation to the radial components of the fluxes evaluated at the interfaces between the shells, F_j± 1/2, and the areas of the interfaces between the shells, A_j± 1/2.
The fluxes and areas are scaled to F_s,ref and L_p^2 respectively.
For each geometry,

    b_j = A_j-1/2F_j-1/2 - A_j+1/2F_j+1/2.

The shell areas differ for the two geometries, given scaled to L_p^2,

    A_1 - 1/2^cylinder   = A_1 - 1/2^sphere = 0,
        
    A_j + 1/2^sphere   = 
            4π( r_j + Δr/2)^2,     j = 1… N-1 
    
            4πr_j^2,     j = N 
    ,
        
    A_j + 1/2^cylinder   = 
            2πh_p( r_j + Δr/2),     j = 1… N-1 
    
            2πh_p r_j,     j = N 
    .


The flux is calculated using the two-point centered finite difference approximation using the concentrations in the adjacent cells.
Unlike in Zeng et al. <cit.>, 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 <cit.> 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 <cit.>.
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 <cit.>.
For interior points, j = 2… N-1,

    ∇^2c_j^ sphere   = 2/r_j.∂c/∂r|_r_j
        + .∂^2c/∂r^2|_r_j≈2/r_jc_j+1 - c_j-1/2Δr
        + c_j-1 - 2c_j + c_j+1/Δr^2,
        
    ∇^2c_j^ cylinder   = 1/r_j.∂c/∂r|_r_j
        + .∂^2c/∂r^2|_r_j≈1/r_jc_j+1 - c_j-1/2Δr
        + c_j-1 - 2c_j + c_j+1/Δr^2.

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

    ∇^2c_1^ sphere   = 3.∂^2c/∂r|_r_1≈ 32c_2 - 2c_1/Δr^2,
        
    ∇^2c_1^ cylinder   = 2.∂^2c/∂r|_r_1≈ 22c_2 - 2c_1/Δr^2.

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

    𝐧·c_N
           = 1/κ∂γ_S/∂c_s≈c_N+1 - c_N-1/2Δr

giving c_N+1≈2Δr/κ∂γ_S/∂c_s + c_N-1, such that

    ∇^2c_N^ sphere   ≈2/r_N1/κ∂γ_S/∂c_s
        + 2c_N-1 - 2c_N + 2Δr/κ∂γ_S/∂c_s/Δr^2,
        
    ∇^2c_N^ cylinder   ≈1/r_N1/κ∂γ_S/∂c_s
        + 2c_N-1 - 2c_N + 2Δr/κ∂γ_S/∂c_s/Δr^2.




 §.§ DAE Consistent Initial Conditions

Because the discretized equations are coupled DAE's, we must begin the system from consistent initial conditions <cit.>.
Although various approaches for initializing this type of model have been developed <cit.>, 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 V material against a lithium metal or graphite anode is at ∼ 100 thermal volts compared to k_BT/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.



§ 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.



 §.§ 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 <cit.>.
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,

    μ_RS = k_BTln( c/1-c) + Ω( 1-2c) - κ/c_max∇^2 c + μ^Θ

where c = c/c_max is a filling fraction, c_max is the maximum concentration of lithium in the active material, and with μ^Θ arbitrarily set to -2 eV defined relative to Li/Li^+.
The regular solution parameter, Ω, is set to 3k_BT_ref, where T_ref = 298 K is the absolute temperature at which the simulation is carried out.
We examine an electrode with particles with radius, R = 1 μm, and we choose the interfacial gradient penalty, κ = 1.16×10^-7 J/m such that the phase interface width is several times the finite volume grid size, 20 nm.
The maximum filling fraction in the active material, c_max = 25 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. <ref>, a piece-wise continuous function defined as equal to μ_RS outside the miscibility gap and μ^Θ within it, the red dashed curve in Figure <ref>.
Note, this could differ from a particular equilibrium filling/emptying path because a particle described by Eq. <ref> can enter metastable regions, leading to equilibrium hysteresis <cit.> 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 (α = 0.5) Butler-Volmer model with a constant exchange current density, i_0 = 1 A/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 <cit.>.
The transport of solid solution material can be described using Fickian diffusion.
For the solid solution, we use constant D_chem=1×10^-16 m^2/s, and for the phase separating case, we use Eq. <ref>, which approximately matches the concentration-independence of the chemical diffusivity in the solid solution regimes <cit.>.
We adjust D_0 until we get a similar dependence of available battery capacity as a function of current (Fig. <ref> (b)), where capacity is defined as the electrode filling fraction when the voltage reaches a 1.5 V cutoff.
We find D_0 = 8×10^-16 m^2/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. <ref>).



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 μm and 20 μ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 <ref> using parameters as determined by Valo̸en and Reimers <cit.> but replacing the conductivity with that measured by Bernardi and Go <cit.>.

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 <ref>), 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 <cit.>.
Of note, in the companion paper, we used Eq. <ref> with constant prefactor, i.e. D_0,iγ_ic_i/γ_,i^d =, 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 <ref>, we note that the models give very different predictions of the behavior at the particle level, shown in Figure <ref> 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 <ref> that when both particles are simulated using the concentration-dependent exchange current density in Eq. <ref>, 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 <cit.>.




 §.§ 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 <cit.>, and we briefly compare the dilute approach previously used by Bazant and co-workers <cit.> with the concentrated model here.
We use the same phase separating material as that studied Section <ref> with constant exchange current density but simulate it in a battery with a longer electrode of 200 μm.
The concentrated electrolyte model is identical to that used in Section <ref>, 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 <cit.>, neglecting concentration dependence of the electrolyte transport coefficients, and using D_ℓ,+ = 2.42×10^-10 m^2/s and D_ℓ,- = 3.95×10^-10 m^2/s, which corresponds to the concentrated solution model with a typical value of D_ℓ = 3×10^-10 m^2/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 <ref>.
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 <cit.>.
As in Figure <ref> (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.



 §.§ Full Cell LiFePO_4/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 LiFePO_4 (LFP) which was developed by Bazant and co-workers <cit.> and used to explain in situ measurements of mosaic phase transformations in porous electrodes <cit.>.
In the anode, we simulate graphite using a 2-layer Cahn-Hilliard reaction model <cit.> as used to capture the experimental intercalation of a single crystal of graphite <cit.>.
Although we have found a variation of this model to better describe the behavior of porous graphite electrodes <cit.>, 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 μm, 20 μm, and 150 μ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 <cit.>,

    μ_LFP = k_BTln( c/1-c) + Ω( 1-2c) + B/c_max,LFP( c - c) - κ_LFP/c_max,LFP∇^2c + μ^Θ_LFP

where Ω = 4.51k_BT, the stress coefficient B = 0.19 GPa, c_max,LFP = 23 M, c is the average filling fraction in the particle, κ_LFP = 5×10^-10 J/m, and μ^Θ_LFP = -3.4 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 <cit.>.
The flux within each layer is given by Eq. <ref>, and the diffusional chemical potential of each is given by

    μ_G,i = k_BTln( c_i/1-c_i)
        + Ω_a( 1-2c_i ) - 2κ_G/c_max,G∇^2c_i
        + Ω_bc_j + Ω_c( 1-2c_i )c_j( 1-c_j )
        + μ^Θ_G

where j i and i ∈{ 1,2 } and c_i represents the filling fraction of layer i.
The parameters for graphite come from refs. <cit.> and are Ω_a = 3.4k_BT, Ω_b = 1.4k_BT, Ω_c = 20k_BT, c_max,G = 28.2 M, κ_G = 4×10^-7 J/m, and μ^Θ_G = -0.12 eV with respect to Li/Li^+.
In both materials, the reactions are governed by the symmetric (α=0.5) Butler-Volmer, Eq. <ref>, using Eq. <ref> for the exchange current density.
For LFP, we set k_0 = 0.16 A/m^2 and use γ_ = 1/( 1-c) as proposed by Bai et al. <cit.>.
For graphite, we set k_0=10 A/m^2 and use γ_,i = 1/( c_i( 1-c_i ) ) as we initially used in Guo et al. <cit.> and reexamined recently to suggest alternate models for practical battery simulation <cit.>.



We present overall cell polarization curves for a selection of currents in Figure <ref>.
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 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 <cit.> and explained theoretically <cit.>.
This behavior can be seen in Figure <ref> (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 <cit.>.
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 <cit.> but is unable to capture experiments on a full porous electrode, where we also used a simplified version of the graphite model <cit.>.



 §.§ Electrochemical reactions


Battery models are typically constructed assuming Butler-Volmer reaction kinetics <cit.>, although recent evidence suggests Marcus kinetics may better describe some electrochemical reactions <cit.>.
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 <cit.>.
A second is that the expression for calculating Marcus style kinetics at an electrode involves an improper integral <cit.>, 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 <cit.>, 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. <ref> 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 <ref>, 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 A/m^2 and R_film = 0.001 Ω m^2, a value similar to that used to match electrochemical data using porous electrode modeling <cit.>.
The MHC kinetics use the same exchange current density and a reorganization energy, λ = 18k_BT, between that approximated for LiFePO_4 <cit.> and values used in other interfacial reactions <cit.>.
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 <cit.> 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 <ref> 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, η_f, is offset from the activation overpotential.
In addition, the prefactor, i_M = k_M/γ_M,, may also have concentration dependence related to concentrated solution effects captured in γ_M,.
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 γ_M,, 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 <ref> (a) and note that it has a dependence on the reduced species (intercalated lithium) that is approximately a square root.
This applies to both c_R and 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 √(c_Rc_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. <ref> with Eq. <ref>.
Thus, for the MHC model, we choose

    i_M = k_Mn( a_Oa_e^n )^1-αa_R^α/γ_√(c_Rc_O),

with n=1 and α=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, η.
We use γ_ = 1/( 1-c_R ) following previous work on LFP <cit.> and supported by recent experiments <cit.>.
This gives an exchange current density with broken symmetry around half filling, as depicted in Figure <ref> (b).

We choose a reaction reorganization energy for the MHC model λ = 18k_BT as in Figure <ref>.


In order to minimize the effect of solid phase transport losses, we modify the phase separating particles from Section <ref> to have D_0 = 1×10^-14 m^2/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 μm with porosity 0.8 and a thin cathode of length 25 μ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 μm with k_0 = 1 A/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 <ref> to R_film = 0.02 Ω 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 <ref>.
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. <cit.>.

We expose the cell to the current profile and corresponding state of charge profile in Figure <ref> (b), and the output voltage profiles are shown in Figure <ref> (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 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 -μ^Θ/e = 2 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 min, the surface concentrations rise sharply, and the exchange current density decreases (Figure <ref> (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 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 <ref> (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η≪λ, the Butler-Volmer and MHC models predict identical results.
The value of λ we used here is between a value calculated for LFP <cit.> and values commonly found in calculations or fits to data for other systems <cit.>.
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.



§ 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 <cit.> applied to porous electrodes <cit.>.
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 <cit.>, adapted for electrochemical systems <cit.>.
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 <cit.> 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 <cit.>, 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 <cit.> 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. <cit.> 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 <cit.> or the double layers of charged porous separators <cit.>, which also activates additional mechanisms for ion transport by surface conduction and electro-osmotic flow <cit.>, 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.



§ 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.








