Charge transport modelling of lithium ion batteries

This paper presents the current state of mathematical modelling of the electrochemical behaviour of lithium-ion batteries as they are charged and discharged. It reviews the models developed by Newman and co-workers, both in the cases of dilute and moderately-concentrated electrolytes and indicates the modelling assumptions required for their development. Particular attention is paid to the interface conditions imposed between the electrolyte and the active electrode material; necessary conditions are derived for one of these, the Butler-Volmer relation, in order to ensure physically realistic solutions. Insight into the origin of the differences between various models found in the literature is revealed by considering formulations obtained by using different measures of the electric potential. Materials commonly used for electrodes in lithium ion batteries are considered and the various mathematical models used to describe lithium transport in them discussed. The problem of up-scaling from models of behaviour at the single electrode particle scale to the cell scale is addressed using homogenisation techniques resulting in the pseudo 2D model commonly used to describe charge transport and discharge behaviour in lithium-ion cells. Numerical solution to this model is discussed and illustrative results for a common device are computed.


Introduction
Lithium-ion batteries are currently one of the most hopeful prospects for large scale efficient storage of electricity for mobile devices from phones to cars. Crucial to their continued improved performance is to understand how novel materials might be effectively exploited in their design. Excellent reviews of the current status of such materials 2 G. W. Richardson et al. are given by Bruce et al. (2008), Choi et al. (2012), Blomgren (2017). Understanding how these materials affect macroscopic battery behaviour is greatly aided by good mathematical models of the transport processes within the battery.
The purpose of this paper is to serve as a guide to charge transport modelling in lithium-ion batteries. Much of the work in this area is due to John Newman and his co-workers who, in a series of seminal publications Newman (1973), Newman & Thomas-Alyea (2012), Doyle et al. (1993Doyle et al. ( , 1996, Fuller et al. (1994b), Ma et al. (1995), Newman et al. (2003), Srinivasan & Newman (2004b), introduced and applied models for these devices that account both for charge transport in the electrolyte as well as solid lithium ion diffusion in the active electrode materials, and use Butler-Volmer reaction kinetics for lithium intercalation and de-intercalation on the electrode/electrolyte surface to couple these two processes together. While these models have been remarkably successful in describing the behaviour of real batteries they are not easily extracted from the literature and, in addition, improvements in understanding of electrode materials have led to significant recent advances in the modelling of lithium transport in electrode particles that can incorporated into this framework. This work aims to provide a relatively concise guide to the subject while at the same time highlighting some of the common modelling pitfalls.
A typical lithium-ion battery (LIB) cell has three regions: i) a porous negative electrode, ii) a porous positive electrode and iii) an electron-blocking separator (see figure  1). Typically the electrodes are comprised of particles, typically a few microns in size, of different (solid) active materials (AMs), that are capable of absorbing lithium into their structure and therefore act as lithium reservoirs. These particles are interspersed with an inert porous polymer binder material, combined with highly conducting carbon black, that acts to hold the electrode particles in place and form conducting links between electrode particles. The AM particle and polymer binder regions are permeated by a lithium electrolyte that serves to transport charge and lithium ions between the AMs of the two electrodes, with direct electrical contact between the negative and positive Charge transport modelling of lithium ion batteries 3 electrodes being prevented by the presence of the porous separator. On the interfaces between the electrode particles and the electrolyte (de-)intercalation reactions, in which lithium ions transfer between the electrolyte and the AM, take place. The AM of the positive electrode shows a greater affinity for lithium than that of the negative electrode so that in a charged cell there is a propensity for a current of positively charged lithium ions to flow from the negative to the positive electrode and thereby establishing a useful potential difference between the electrodes. The rates at which these (de-)intercalation reactions take place on the electrode/electrolyte interfaces are key to the electrical behaviour of the battery, and are typically described by a Butler-Volmer relation which gives the interfacial current density, from AM to electrolyte, in terms of the potential jump between the AM and the electrolyte and the lithium concentrations in the AM and the electrolyte.
In addition to discussing transport models we also propose a new formulation for the Butler-Volmer relation, for current flow between the AM and electrolyte, with the aim of ensuring a model that is able to simulate extreme cases, where classical formulations lead to physically unrealistic lithium distributions. More precisely, we address the issues of the limiting conditions in which the electrodes are close to being fully intercalated (or fully depleted) or in which the electrolyte concentration is very low.
The LIB is great examplar of multiscale and multiphysics systems. Accurately predicting battery pack (formed from many individual cells) behaviour depends on having appropriate represenatations of the physics and chemistry at a range of smaller scales, including that of individual cells within a pack, individual electrodes within each cell, individual particles within electrodes and at the level of the atomic structure of the materials making up the particles. The behaviour at many of these lengthscales is dictated by a myriad of interacting phenomena including electrochemical ones, but also thermal and mechanical ones. The seminal models developed by Newman and his co-workers span the scales of individual electrodes particles upto individual cells and are largely focussed on the electrochemical behaviour. Many authors including Ranom (2014), Schmuck (2017) and Franco (2013) discuss the challenges associated with formulating mathematical models that couple phenomena pccuring at vastly disparate lengthscales. There is also current impetus to experimentally characerise and develop new modelling approaches to better understand chemical structure within electrode materials Harris et al. (2017) and Kim et al. (2004), mechanical and thermal effects Oh et al. (2014), chemical degradation Sethurajan et al. (2019), Wang et al. (2018), Birkl et al. (2017), Monroe & Newman (2003) and Nishikawa et al. (2011), and to also develop battery management systems to optimally control batteries for use in electric vehicles Lu et al. (2013) and Kim et al. (2014).
The outline of this work is as follows. In §2 we discuss charge transport models for the electrolyte. We begin, in §2.1, with the simplest description, dilute electrolyte theory, and show that it cannot adequately describe electrolyte data at the typical concentrations encountered in real batteries. This motivates us to consider Newman's moderatelyconcentrated electrolyte theory Newman (1973), Newman & Thomas-Alyea (2012) in §2.2, which forms the basis for much of the battery electrolyte modelling currently being undertaken and to show how this fits to real data for the common electrolyte LiPF 6 Valøen & Reimers (2005). In §3 lithium transport within AM electrode particles is briefly reviewed while in §3.2 a new formulation for the Butler-Volmer relation is proposed. In §4 the various strands of the battery chemistry, described in the previous sections, are brought together to formulate a macroscopic device scale model for an entire cell. This is accomplished via homogenisation method set out in §4.2. In §5 we present a selection of solutions of the device scale model for a common modern device configuration. Finally in §6 we review the key insights of this work.

Modelling the electrolyte
Here we begin, in §2.1, by considering the theory of very dilute electrolytes, often termed Poisson-Nernst-Planck (PNP) theory. Because the theory is commonly encountered in modelling semiconductors, is relatively straightforward and physically appealing, it is useful to highlight some of the peculiarities associated with charge transport modelling in batteries. These include, charge neutrality, the use of electrochemical potentials, and the measurement of the electric potential with respect to a lithium reference electrode. However the price of simplicity is that dilute theory does not describe battery electrolyte behaviour particularly well. Many of these limitations are overcome in Newman's theory of moderately-concentrated electrolytes, which we review in §2.2. This theory is considerably more involved than the dilute theory and in practice requires that various functions be fitted to data directly measured from the electrolyte under consideration. However, within these limitations, it does provide a good description of most electrolytes formed by dissolution of a salt in a solvent. We also hope that by introducing the peculiarities of notation associated with battery electrolyte modelling in the context of the simpler dilute theory it will make it easier for the reader to follow the more complex moderatelyconcentrated theory.

Dilute electrolytes
We consider an electrolyte composed of a solvent, a negative ion with molar concentration c n and charge z n e, and a positive positive ion with molar concentration c p and charge z p e (where e is the elementary charge and z n , z p are integers accounting for the charge state). This general binary electrolyte can be easily studied but because the purpose of this section is to give a simple introduction in the rest of this article we focus on a 1:1 electrolyte with a generic negative ion and a positive lithium ion Li + , so that z n = −1, z p = 1.
Because of the long timescales over which batteries are typically charged and discharged it is entirely reasonable to neglect magnetic effects and assume that the electric field E is irrotational (i.e. ∇ × E = 0) so that it can be written in terms of an electric potential φ, via the relation E = −∇φ. (2.1) Considering the charge within the system then gives Poisson's equation where ε is the permittivity of the electrolyte.
Charge transport modelling of lithium ion batteries

5
Since the ions in battery electrolytes do not react with each other or with the solvent, conservation of the two ion species implies that ∂c n ∂t + ∇ · q n = 0, and ∂c p ∂t where q p and q n are the fluxes of positive and negative ions, respectively. We can also write this using q p = c p v p and q n = c n v n where v p and v n are the average velocities of the respective species. In the dilute limit the electrolyte solvent is assumed to be stationary and the ions to move in response to thermal diffusion and electric fields; interactions between ions are neglected. The component of the average velocity of a lithium ion due to the electric field, v ep , is given by balancing the force eE, exerted on it by the electric field, with the viscous drag force v ep /M p , exerted on it by the solvent (here M p is the mobility of the lithium ion). Thus the advective lithium ion velocity due to the electric field is v ep = −eM p ∇φ and in a similar manner the average negative ion velocity can be shown to be v en = eM n ∇φ. In addition to the advective fluxes (c p v ep and c n v en ) both ion species diffuse, in response to random thermal excitations. This gives rise to Fickian fluxes (for positive and negative ions) of size −D p ∇c p and −D n ∇c n , respectively, where D p and D n are the respective diffusion coefficients. To highlight the difference between these diffusion coefficients and those used in the Stefan-Maxwell theory that we will review in §2.2, let us point out that D p (respectively, D n ) is the diffusion coefficient for positive (respectivley, negative) ions in a mixture of solvent and negative (respectively, positive) ions. The total ion fluxes, q n and q p , are obtained by summing their advective and diffusive components, so that q n = c n v n = −D n ∇c n − e kT c n ∇φ , (2.4) Here we have substituted for ion mobilities in terms of the diffusion coefficients by using the Einstein relations M p = D p /kT and M n = D n /kT (where k is Boltzmann's constant).
Since this theory is applied in a chemical setting it is more usual to write these equations in terms of Faraday's constant F and the universal gas constant R which, on noting that e/k = F/R, leads to the alternative expressions The equations governing the three variables, c n , c p and φ, are thus (2.2), (2.3) and (2.6).

Double layers and charge neutrality.
It has long been recognised that (see e.g. Newman & Thomas-Alyea (2012)), at the concentrations typically encountered in practical electrolytes, there is almost exact charge neutrality. This implies that there is a balance between the concentrations of positive and negative charges, throughout the vast majority of the electrolyte. The exception to this rule is in the so-called double layers which lie along the boundaries of the electrolyte region and are typically extremely thin, with widths less than a few nanometres. This observation can be justified mathematically by non-dimensionalising equations (2.2), (2.3) and (2.6) and conducting a boundary layer analysis in terms of the small dimensionless parameter which measures the ratio of the Debye length (i.e. the typical width of a double layer) to the typical dimension of the electrolyte. 1 The result of such an analysis (see, for example, Richardson (2009)) is that, with the exception of the double layers, c n must be almost exactly equal to c p . The physical meaning of this fact is that the attraction between charges is very strong compared to any space charge that the electric field may create. Therefore, a very good approximation to (2.2) is which is usually called the charge neutrality condition. As we shall discuss later, because (2.7) has neglected the derivatives that were in (2.2), the model needs fewer boundary conditions at the edges of the electrolyte (i.e. only two boundary conditions are required on the electrolyte 'surface' rather than the three needed for the full system).

The approximate equations.
In line with the discussion above we introduce a single concentration c by taking (2.8) and substitute this into (2.3) and (2.6) to obtain the approximate charge-neutral equations ∂c ∂t + ∇ · q n = 0, and q n = −D n ∇c − F RT c∇φ , (2.9) ∂c ∂t + ∇ · q p = 0, and q p = −D p ∇c + F RT c∇φ . (2.10) A common approach to studying this problem is to assume that the ionic diffusivities D n and D p are constant and rewrite the system by adding (2.9a) multiplied by D p to (2.10a) multiplied by D n . On substituting for q n and q p this yields a diffusion equation for c, of the form where D eff is termed the effective ionic diffusivity. Useful physical insight can be found using an alternative formulation of (2.9) and (2.10) by introducing the electric current density, j, defined in terms of the ion fluxes by (2.12) Using this concept, a version of Ohm's Law may be obtained by subtracting (2.9b) from (2.10b), while a charge conservation equation may be found by subtracting (2.9a) from (2.10a). These may be written in the form Here t + is referred to as the transference number andκ(c) is referred to as the electrical conductivity of the electrolyte (c.f. the standard form of Ohm's law is j = −κ∇φ). We use the notationκ to distinguish the electrical conductivity in the dilute limit from the same quantity in moderately concentrated theory, see (2.46). The alternative formulation of (2.9) and (2.10) mentioned above is then (2.11), (2.13) and (2.14).
Assuming that the ionic diffusivities D n and D p are constant implies (i) that transference number t + is constant, (ii) the effective ionic diffusivity D eff is constant and (iii) electrolyte conductivity κ grows linearly with electrolyte concentration c. All three of these quantities are readily measured experimentally for real electrolytes. For most electrolytes transference number is usually found to remain close to constant (with the exception of some polymer electrolytes e.g. Doeff et al. (2000), Fauteux (1988)), electrolyte diffusivity usually decreases relatively weakly with concentration, except at very dilute concentrations, but the growth of electrical conductivity with concentration is far from linear. Examples of the experimentally measured concentration dependence of D eff and κ (from Valøen & Reimers (2005)) are plotted in Figure 2 for the battery electrolyte LiPF 6 . Notably at the typical concentrations used in batteries (roughly 1 molar for LiPF 6 ) electrical conductivity is nearly constant and often lies close to its maximum value, and is thus not well-approximated by the linear expression in (2.15). The explanation given for this poor fit is usually that even at relatively dilute concentrations there is a significant drag between ions of opposite charge. The reason for this is that two ions of opposite charge that lie close to each other experience a significant electrostatic attraction that negates, to a large extent, the effects of the global electric field which is trying to drive the ions in opposite directions. This observation has motivated Newman Newman & Thomas-Alyea (2012) to use Stefan-Maxwell theory for a multi-component solution to describe the charge transport behaviour of electrolytes. A summary of the modelling assumptions made in this theory is given in §2.2.

2.1.3
The dilute theory in terms of electrochemical potentials.
The dilute model (2.9)-(2.10) can also be written in terms of the electrochemical potentials, µ n and µ p , of negative and positive ions respectively. This is the preferred notation for the ion conservation equations in the electrochemical literature. For an electrolyte that is formed by ideal salt solution the electrochemical potentials are given by where c T is the total molar concentration of the electrolyte and, for a dilute solution, is approximately equal to the solvent concentration. The first term on the right-hand sides of both expressions in (2.16) is the standard state potential (per mole of the species), while the second is the entropy of mixing (per mole of the species) and the final term is the electrostatic potential (per mole of the species). Using this notation the conservation equations (2.9) and (2.10) can be written in the form ∂c ∂t + ∇ · (cv n ) = 0, and v n = − D n RT ∇µ n , (2.17) Here the average velocities of the two species, v n and v p , are obtained by multiplying the gradient of the electrochemical potentials by the species mobilities, D n /RT and D p /RT . This formalism extends to nonideal salt solutions and to multicomponent systems. In §2.2 this approach of using electrochemical potentials is extended to moderately-concentrated electrolytes.
2.1.4 The potential measured with respect to a lithium electrode.
The dilute theory as formulated above is at odds with the electrolyte theory used by Newman & Thomas-Alyea (2012); here, the factor in front of the concentration gradient in (2.13), the constitutive law for the current, is 2(RT /F )(1 − t + ) rather than (RT /F )(1 − 2t + ) as above. As pointed out in Ramos (2016), this has generated some confusion in the literature. The explanation for this discrepancy (as initially demonstrated by Ranom (2014) and subsequently in Bizeray et al. (2016)) is that the theory in Newman & Thomas-Alyea (2012) is formulated in terms of ϕ, the electric potential measured with respect to a reference lithium electrode, rather than φ, the true electric potential. In electrochemical applications the potential in an electrolyte is typically measured by inserting a reference electrode of a pure compound. The potential measured depends on the composition of the reference electrode through its chemical potential. Since in lithium battery applications the reference electrode used is nearly always made of lithium, and since much of the data used to calibrate battery models is collected using a lithium reference electrode, it makes sense to use the potential measured with respect to a lithium electrode. Note that φ, the true electric potential, is not a readily measured quantity. In order to switch between ϕ and φ we recall that there is a reversible reaction that occurs on the surface of the electrode between intercalated lithium in the electrode and lithium ions in the electrolyte, given by Here the subscript (l) denotes a reactant within the electrolyte and (s) one within the electrode. Typically the current flow into a reference electrode can be assumed to be sufficiently small that this reaction is in quasi-equilibrium. It follows that the electrochemical potentials of compounds on both sides of this equation are equal (i.e. µ p + µ e − = µ Li ).
Since neither the concentration of electrons nor the lithium within the electrode change, this equality implies Charge transport modelling of lithium ion batteries 9 (with +µ 0 e − being the Fermi level of the lithium electrode), which rearranges to (2.20) Using (2.20) to substitute for φ in (2.13) yields the electrolyte Ohm's law found in the Newman theory (2.21) The effect of the difference in the two different potentials is now readily seen by comparing (2.21) with (2.13). Therefore (2.21), together with (2.11) and (2.14) is a system of equations equivalent to the charge-neutral equations (2.9) and (2.10).

Moderately-concentrated electrolytes
Motivated by the confusion in the literature highlighted in Ramos (2016), this section reviews, in detail, a commonly used model for moderately-concentrated electrolytes, presented in Newman & Thomas-Alyea (2012), which is applicable to most electrolytes consisting of a salt dissolved in a solvent but not to ionic liquids. In most practical battery systems ion transport takes place through electrolytic solutions which do not behave as ideal dilute materials as demonstrated primarily by the concentration dependence of their conductivity, see Valøen & Reimers (2005), but also by activity coefficient measurements, e.g. those in Samson et al. (1999). In order to capture this non-ideal behaviour it is necessary to consider not only ion/solvent interactions (as is done in the PNP theory of ideal electrolytes as covered in §2.1) but also interactions between the ionic species.
Inter-ionic interactions are significant, in even relatively dilute solutions, because the local attraction between oppositely charged ions result in a propensity for ions of opposite charge to lie close to each other, which reduces their mobility in an electric field. This, in turn, reduces the ionic conductivity, see Samson et al. (1999). Models of batteries that use electrolytes in this moderately-concentrated regime have been pioneered, and applied successfully to a variety of systems. For example see the series of seminal works by John Newman and his co-workers: Newman (1973), Newman & Thomas-Alyea (2012), Doyle et al. (1993), Newman et al. (2003), Doyle et al. (1996), Srinivasan & Newman (2004b). The electrolyte theory reviewed below is based on the Stefan-Maxwell equations (see, for example, Bird et al. (2002)), which describe transport in a mixture (including diffusion) in terms of the drag coefficients between its various components.

Stefan-Maxwell equations
We start by briefly considering the general case in which the electrolyte is comprised of N (ionic and solvent) species. Newman & Thomas-Alyea (2012) uses the Stefan-Maxwell equations as the foundation of concentrated electrolyte theory. These relate the drag force acting on a component in a mixture to its relative velocity with the other components. It is by balancing these drag forces with gradients of electrochemical potential of a species and gradients in the fluid pressure that we obtain the average velocity of each species and in turn its flux. This combined with statements of conservation of species form the equations of moderately concentrated electrolyte theory. We note that other mechanisms in addition to the interspecies drag, electrochemical potential and pressure forces can also cause mass transfer and we briefly mention these without detailed dicussion.
In order to derive the force acting on each component of the mixture as a consequence of the pressure gradient it is neccessary first to obtain an equation of state that relates the concentrations of the N species forming the mixture. According to Liu & Monroe (2014) for most electrolytes it is usually a good approximation to assume that each species has constant molar volume. This is equivalent to the assumption that the volume occupied by one mole a given species remains fixed whatever the composition of the mixture. On denoting the molar volume of the i'th species by H i we obtain the following equation of state relating the molar concentrations (2.22) As we shall see this relation allows us to write down the force on a species arising from the pressure gradient in the mixture. We note that electrolytes may have mechanical properties, such as acting as a viscous fluid or an elastic solid, which can create additional forces. We do not consider these but they can contriubute significantly in certain situations.
The mutual friction force between species i and j is assumed to be proportional to the friction forces arising from velocity differences between the species. Furthermore this force is proportional to the mole fraction, χ k , of each species (see Bothe (2011)) as defined by Here c k is the molar concentrations of species k and c T is the total molar concentration of all species in the solution. The Stefan-Maxwell equations give a relation betweend i , the drag force exerted on species i, per unit volume of mixture, and the velocities of the various species. In light of the above comments the drag force on the i'th species (per mole unit volume) is taken to depend linearly on the velocity differences between species, and modelled (see Bothe (2011)) by the expression where v k is the velocity of species k and RT k ij is the drag coefficient on one mole of species i moving through pure species j. Note that k ij is symmetric (i.e. k ij = k ji ) because the drag exerted on species i by species j is equal and opposite to that exerted on species j by speciesi. Here the Maxwell-Stefan inter-species diffusivity is related to k ij by the Einstein relation so that D ij = 1/k ij . The drag forced i , is balanced by motive forces (per unit volume) down gradients in the electrochemical potential µ i and down gradients in the pressure p Here the electrochemical potentials µ i may be rewritten in terms of the chemical potentialsμ i and the electric potential φ in the standard fashion where z i is the valence of species i. The force balance (2.25) equations are supplemented by the standard conservation equations, which can be written as A force balance on the entire mixture can be obtained by adding together the N relations (2.25) and noting that By substituting for N i=1 H i c i from the equation of state (2.22) and for the electrochemical potential from (2.26) we obtain the following expression for the total force balance (2.28) It is straightforward to show from the Gibbs-Duhem relation between the chemical potentials, namely This result implies that gradients in the chemical potential, in isolation, do not, as might be expected, lead to a net force on the mixture. In turn this means that the pressure equation (2.28) can be simplified to in which ρ represents the charge density of the mixture. It is worth making some brief comments about (2.30) which gives an intuitively appealing balance between electrostatic forces and pressure forces acting on the mixture. By taking the curl of (2.30), it is clear that it can only be satisfied if ∇ρ × ∇φ = 0 (i.e. the gradient of the charge density lies parallel to the electric field E = −∇φ). There are two special cases where this relation is automatically satisfied which are particularly relevant here. The first of these is where ρ = 0 (or is negligible), and in this instance no pressure gradient is required to balance the electric force on the mixture and so results in a spatially uniform pressure; this is the case that applies to charge neutral elecrolytes and so is pertinent to battery modelling. The second special case is where the problem is strictly one-dimensional when the pressure gradient simply counteracts the electrical force on the mixture. This is the same situation as occurs in one dimensional fluid flow where incompressibility dictates that the motion is determined by the concentrations without reference to any forces and is discussed in the context of other such multiphase systems in Drew (1983).
In more general cases the force balances, described above in (2.25), are too naive and need to be supplemented by multiphase viscous dissipation terms as described in Drew (1983). Such issues, however, are beyond the scope of this work but they are addressed in this context in Liu & Monroe (2014).
A final comment about the general moderately-concentrated problem is that the electrochemical potential of the i'th species µ i has essentially the same form as that written down for a dilute 1:1 solute in (2.16). The only modification is that the species mole fraction c i /c T is replaced by its activity a i to reflect the fact that the solution is non-ideal. Hence we find where z i is again the charge state of the i'th species.

Stefan Maxwell equations for a binary 1:1 electrolyte
We now take the general theory of §2.2.1 and restrict attention to the case of a 1:1 electrolyte comprised of a solution of Li + ions and a generic negative counter ion species dissolved in a single solvent species. Although battery electrolytes are often based on rather complex solvent mixtures, which are usually closely guarded industrial secrets, this approach provides a reasonable description of many battery electrolytes. In Figure  2, we parameterize the model against experimental data for the most common lithium ion electrolyte LiPF 6 in 1:1 EC:DMC from Valøen & Reimers (2005), treating the two component solvent (EC:DMC) as if it they were a single solvent.
In line with Newman & Thomas-Alyea (2012) we denote the three species making up the electrolyte, namely the solvent, the Li + ions and the generic counterions, by the subscripts i = 0, p, n, respectively. We have therefore z 0 = 0, z p = 1 and z n = −1. Combining (2.23)-(2.25) and expanding in component form yields where, by using (2.23) (i.e. c T = c 0 + c p + c n ) and (2.24)-(2.25), the drag coefficients can be expressed in the form Henceforth we can omit (2.34), noting that the choice of physically realistic functions µ k (with k = 0, n, p) ensures this is satisfied. Using (2.31) the electrochemical potentials of Charge transport modelling of lithium ion batteries 13 the ion species are As in dilute theory, charge neutrality can be assumed so that we can write (2.37) As discussed above, where the electrolyte is charge neutral (i.e. c n = c p ), the solution to the pressure equation (2.30) is such that p is constant and hence the pressure gradient terms in (2.32)-(2.34) vanish. We now seek to find a constitutive equation for the current density j. We will broadly follow the derivation given in Newman & Thomas-Alyea (2012) but attempt to clarify their argument. As in the dilute case, the current density is given by (2.12), which can be rewritten as (2.38) Substitution of (2.38) into (2.32) and (2.33), taking account of (2.37) and the fact that the pressure gradient terms vanish, yields On substituting for v p − v n (in terms of j) from (2.38), and for µ p and µ n from (2.36), and rearranging the resulting expression we obtain the following expression for j: are the transference number of the (positive) lithium ions with respect to the solvent velocity, and the conductivity of the electrolyte as a function of the concentration. Note how this equation for j, and the definition of transference number and conductivity compare to the dilute version given in (2.13)-(2.15). We now rewrite equations (2.42)-(2.43) in a form that avoids the use of a chemical potential for each ion species and instead uses a chemical potential for the entire electrolyte µ e . Without any loss of generality we can express v p and v n in the form for any function α. Our procedure is to substitute v p − v n = j/(F c) in the final terms of these expressions and then to replace v n and v p , everywhere else on the right-hand side of these expressions, from (2.42) and (2.43). Choosing α = t 0 + , as defined in (2.46), allows the electric potential φ to be eliminated from the term (1 − α)v p + αv n . The expressions for v p and v n , in (2.47), can then be rewritten in the form These equations can be further simplified by introducing the electrolyte chemical potential µ e (c) (a function of electrolyte concentration only) and the chemical diffusion coefficient D, as we did in (2.11), defined by and noting that, from (2.36), µ e = 1 2 (µ p + µ n ). We then find that (2.48)-(2.49) can be rewritten in the form Notably, D n0 and D p0 can vary with concentration independently, without affecting the preceding analysis. It follows that D and t 0 + may also vary independently as functions of concentration.
The remaining equations governing the behaviour come from considering conservation of the ions and solvent as well as the volume they occupy. The mass conservation equations for the two ion species are the same as (2.27), namely Taking the differences of these two equations, and substituting for v p − v n from (2.38), yields an equation for current conservation and this can be compared to its dilute theory counterpart given in (2.11)). These equations couple to the mass conservation equation for the solvent which, from (2.27), is given by Finally we require an equation of state. On using the charge neutrality condition c n = c p = c in (2.22) we find that this can be written in the form In one dimension we now have sufficient equations to specify the problem; these are composed of the five equations (2.45), (2.54), (2.55), (2.56) and (2.57) for the five variables c, c 0 , φ, j and v 0 . Note that, in more than one dimension, these equations are not sufficient, at least not in general, as can be seen by multiplying (2.53) by H c and adding to (2.53) multiplied by H 0 . On using (2.57) to eliminate the time-derivative from the resulting equation this yields the scalar PDE which in multiple dimensions is insufficient to determine the vector v 0 . As has been described in Drew (1983) this conservation equation needs to be supplemented by a momentum equation, but this is beyond the scope of this work.
Having posed the governing equations, an additional simplifying assumption is commonly made, which appears to be adequate for most solvent based battery electrolytes. This consists of assuming that the electrolyte is sufficiently dilute so that c 0 ≈ c T in which case (2.57) can be approximated by c 0 ≈ 1/H 0 and the solvent velocity is small (i.e. |v 0 | |v p | and |v 0 | |v n |). This limit is equivalent to the approximation This assumption also avoids the complication, that arises in multiple dimensions, of requiring to be supplemented by a momentum equation. Note however that, even in this small concentration limit, we still include interphase drag between negative and positive ions. This is because interphase drag is very often significant even at quite at low concentrations (often as low as 0.1 molar which is usually a small mole fraction). It arises because positive and negative ions interact via a strong long-range force (the Coulomb force).

2.2.3
The potential measured with respect to a lithium electrode.
As in the dilute case it is useful to reformulate the model in terms of ϕ, the potential measured with respect to a lithium electrode. Once again we assume that the reaction (2.19) occurring on the surface of the reference electrode is in quasi-equilibrium, so that sum of the electrochemical potentials of compounds on each side of (2.19) are identical (i.e. µ p + µ e − = µ Li ). However, owing to the slightly different definition of µ p in the moderately-concentrated case (compare (2.36) with (2.16)), this leads to a modified version of the relation between φ and ϕ, namely which, on referring to the definition of the electrolyte chemical potential in (2.50), can be re-expressed as This is a clearer way to write Ohm's law than (2.45) because it allows j to expressed solely in terms of the measurable quantities ϕ and µ e (notably the activities of the individual ion species a n and a p are not directly measurable). This expression differs slightly from that typically presented by Newman and co-authors because they tend to use the dilute limit of (2.61), namely (2.21).

Summary: model for a moderately-concentrated electrolyte
The moderately-concentrated theory is relevant in those situations where the salt concentration is small compared to the solvent concentration so that c c T , the interactions forces are large (D n0 + D p0 ) D pn , and the solvent velocity can be neglected v 0 ≈0. In this limit the moderately-concentrated charge transport model (2.54), (2.55) and (2.61) takes the form ( 2.65) Note that if we consider the limit of a dilute electrolyte with c c T , and the interaction forces are not large (D n0 + D p0 )/D pn = O(1), then we revert to the dilute solution model (2.11) and (2.21) provided that we identify D p0 with D p and D n0 with D n . It is worth noting that (2.62)-(2.64) can be rewritten in the form where a e (c) = (a n a p ) 1/2 is the activity coefficient of the electrolyte (such that µ e (c) = µ 0 e + RT log(a e (c))) and the effective diffusivity is given by An implicit assumption in Newman's formulation of the equations is that the salt solution behaves as an ideal solution so that a e (c)/a e (c) = 1/c. To fit the model it is necessary to determine the lithium diffusivity D eff (c), the electrolyte conductivity κ(c) and the transference number t 0 + (c). Doing this function fitting leads to a relatively robust way of modelling the electrolyte for engineering applications. An example of fitting the phenomenological model functions D eff (c) and κ(c) to data is shown in Figure 2 for the electrolyte LiPF 6 in 1:1 EC:DMC at T = 293K (this data comes from Valøen & Reimers (2005)). For this electrolyte the transference number of lithium ions, with respect to the solvent velocity, is found to be approximately constant with t 0 + = 0.38. Lines represent the fit to to the experimental data (circles) taken from Valøen & Reimers (2005). The fitted function for the diffusivity and conductivity are given by D eff (c) = 5.3×10 −10 exp(−7.1×10 −4 c) and κ(c) = 10 −4 c(5.2−0.002c+2.3×10 −7 c 2 ) 2 respectively.

Lithium transport and electric current flow through the electrode particles
Here we review the modelling of lithium transport in individual electrode particles and current transport through the solid matrix formed by the agglomeration of electrode particles, polymer binder material and conductivity enhancers (such as carbon black). Transport of lithium through the microscopic electrode particles is typically slow; the diffusion timescale for a lithium ion traversing a microscopic electrode particle frequently being comparable to, or even longer, than that required for a lithium ion to traverse the whole cell in the electrolyte. Since the concentration of lithium ions at the surface of the electrode particles strongly influences the rate at which lithium ions are intercalated into the electrode particles from the electrolyte (or vice-versa) a battery charge transport model must treat both microscopic transport of lithium through the electrode particles and its macroscopic transport thorough the electrolyte. The coupling between these micro-and macro-scale transport processes occurs via a reaction rate condition which specifies the rate of lithium intercalation (or de-intercalation) at the particle surface in terms of (a) the lithium concentration c s on the surface, (b) the potential difference φ s −φ between the particle surface and the adjacent electrolyte 3 , and (c) the lithium concentration c in the adjacent electrolyte. The condition that is typically used in practice is the Butler-Volmer relation (see, for example, Bockris & Reddy (1970)), which accounts for both the reaction rates and the effects of the double layer, and is based on the quantum mechanical Marcus Theory described in Marcus (1965). Usually the potential in the electrode matrix φ s and the corresponding current density flow j s through the matrix is modelled by the macroscopic Ohm's Law where κ s is the effective conductivity of the electrode matrix (formed of electrode particles, binder and conductivity enhancer).

The standard approach
Here we set out the standard approach to modelling lithium transport within the electrode particles of a lithium ion cell and the current transfer process between the electrode particles and the surrounding electrolyte. This approach is widely adopted in the modelling literature (e.g. Dargaville & Farrell (2010), Doyle et al. (1993Doyle et al. ( , 1996, Fuller et al. (1994b), Ma et al. (1995), Newman et al. (2003), Srinivasan & Newman (2004b)) although it has recently been challenged by an alternative approach which is discussed in §3.4. The central idea of the model is to consider a one dimensional problem between the two current collectors describing behaviour of the electrolyte and the charge transport in the solid electrodes. The lithium motion within the particles of the electrodes is on a much smaller scale and movement of this is described using a separate dimension representing position within each particle. Hence we will use x to represent the position between the current collectors and r to represent the position within any particle. The problem is therefore often referred to as a pseudo two-dimensional model. More precisely, one might describe the structure of the model as being multiscale, where both the micro-and macroscopic models are one-dimensional. Here we present the basic model but later we describe how it can be systematically derived.

Microscopic lithium transport models in individual electrode particles
At its very simplest, transport of intercalated lithium within electrode particles is modelled by linear diffusion in an array of uniformly sized (radius a) spheres (see e.g. Doyle et al. (1996)). Hence the lithium concentration in an electrode particle at position x within the electrode, c s (r, x, t), evolves according to where r measures distance from the particle's centre and D s is the solid phase Li diffusion coefficient. This model for lithium transport in the electrode particles is typically coupled to the processes taking place in the adjacent electrolyte through a Butler-Volmer relation, which gives the transfer current density j tr flowing out through the surface of the particle, in terms of the lithium concentrations on the particle's surface c s , the adjacent electrolyte concentration c and the potential difference between the electrolyte and the electrode particle ϕ−φ s . Following Faraday's laws of electrolysis, the flux of lithium on the particle surface is proportional to the current density, leading to the following boundary condition on the electrode particle surface −D s ∂c s ∂r = 1 F j tr at r = a.
Lithium transport within electrode particles is often better described by nonlinear diffusion (with D s = D s (c s )), an approach that is used in Karthikeyan et al. (2008) (2015), a positive electrode material discussed in detail in §3.3.2. We will also discuss in §3.3 and §3.4 more complicated models of lithium transport applicable to cases where phase separation occurs or the behaviour is highly anisotropic.

The Butler-Volmer relation
The transfer current density j tr is the normal component of the current density on the particle surface, from the electrode particle into the surrounding electrolyte, and can be expressed in terms of a Butler-Volmer relation of the form 3) with U eq depending usually only on c s . Here i 0 (c s , c) is the exchange current density and U eq (c s , c) is the open circuit potential (OCP). Both i o and U eq depend on the electrode material while the dimensionless constants α a and α c are anodic and cathodic transfer coefficients lying between 0 and 1, and are conventionally both taken to be 1/2. Note that when measurements of the OCP of a material are taken, they are done so with the use of a reference electrode which is typically lithium. Hence, the potential difference across the double layer is difference between the electric potential φ s in the solid and ϕ in the electrolyte (the potential measured relative to a lithium electrode). Note this this has been assumed in writing (3.3) and as in line with the discussion given below (3.2).
The OCP is found by considering equilibrium, when j tr = 0, and measuring the potential difference between the electrolyte and the electrode particle, φ s − ϕ, for various different levels of the concentrations. For most materials there is a maximum concentration of Li that can occur, denoted by c s,max , and the OCP is typically plotted as a function of y = c s /c s,max , which is the Lithium stoichiometry. The OCP varies widely depending on the electrode material and to illustrate this in Figure 3(a) it is plotted for the lithiated graphite Li x C 6 while in Figure 3(b) it is plotted for lithiated iron phosphate Li y FePO 4 . Note the OCP for Li y FePO 4 has a large almost flat plateau symptomatic of a two-phase state within the material.   Srinivasan & Newman (2004b)). Each is shown as a function of Lithium stoichiometry.
The exchange current i 0 is much less well quantified in experiments and commonly is taken to be given by the formula where k is a kinetic rate constant.
Charge transport modelling of lithium ion batteries 21

Appropriate forms for the Butler-Volmer relation
The Butler-Volmer relation is widely used in the literature (see, for example, Doyle et al. (1993), Fuller et al. (1994c), Gomadam et al. (2002), Smith & Wang (2006a,b), Smith et al. (2007), Chaturvedi et al. (2010), Kim et al. (2012)), and here we focus on why it is necessary to take care in selecting the functions used to describe U eq and i 0 . In the literature (see, for instance, Weng et al. (2014)) U eq is typically fit to some polynomial, or other family of functions, such as exponentials. Furthermore, as explained in Weng et al. (2014), most of the fitting is done to data points which lie in the middle of the stoichiometry range, namely 0.1 -0.9 where many batteries operate, thereby avoiding extreme cases which are prone to triggering failure and/or accelerated degradation, see Wang et al. (2012), Birkl et al. (2017). Arguably more important than the behaviour for intermediate charge states is the what happens when the electrodes get close to being fully intercalated c s = c s,max , or fully depleted c s = 0, or where the electrolyte approaches depletion c = 0. To discuss the allowable behaviours in these cases we consider the scenario in which c s → 0 and then indicate how the same ideas can be extended to apply to the other limiting cases.
The basic difficulty is that as c s → 0 we need the Li flux out of the particle to go to zero to prevent predicting nonphysical negative concentrations in the solid. However, if we start with a depleted particle and do not allow Li to enter, then it will never charge. We need a model that avoids both physically unrealistic situations. It is of course always possible to simply chose i 0 and U eq and then impose some switching logic to turn the flux on or off as required to avoid such problems. However, it is preferable to have a Butler-Volmer relation that incorporates mechanisms that automatically ensures such physically unrealistic situations cannot occur.
It is common to consider the two parts of a Butler-Volmer relation in (3.3) as representing reverse and forward reactions, which are given by the first and second term of the right hand side of that relation, respectively. As c s → 0, the reverse (anodic) reaction should be first order, while the forward (cathodic) should be bounded and positive (like a zeroth order reaction). Such behaviour can be readily achieved by taking i 0 and U eq to have the local form, when c s → 0, where k 1 and k 2 are positive constants (which can be different for each electrode and can also change inside non-homogeneous electrodes). Such a formulation near this extreme of the concentration will automatically ensure the flux cannot become positive as the concentration reduces, thereby avoiding negative concentrations, and the flux can be finite and positive so the particle can be charged from a completely depleted state. Note that a reverse reaction of order greater than one can be considered but this is not usually used. Such a local behaviour of the Butler-Volmer condition is included in very few formulations, see Doyle et al. (1993) and West et al. (1982), and, without explanation of why this choice is taken. A much greater number of articles in the literature do not include this local behaviour and, also, use a formulation of U eq independent of c. These are not wellsuited for extremely low electrolyte concentrations. Nevertheless, it is useful to model 22 G. W. Richardson et al. these cases since, for example, during high discharge rates, the electrolyte can be almost depleted of lithium in some parts of the cell (see Bizeray et al. (2015)).
One reasonable way to rewrite the Butler-Volmer relationship that emphasises the forward and reverse parts of the reactions is to introduce two strictly positive bounded functions i a (c s ), and i c (c s , c) where where u eq is a bounded function, so that the Butler-Volmer equation (3.3) for the transfer current can be written in the form ( 3.8) The functions i a and i c must be strictly positive except at the extremes of the concentrations and locally these functions must have the following behaviours: (1) The function U eq , the OCP, also depends on the electrolyte concentration c (instead of being independent of it).

13)
Charge transport modelling of lithium ion batteries

23
(4) If c s → c s,max , with c = 0, then i 0 → 0 and U eq → −∞ and (3.14) These conditions ensure that the lithium flux is constrained to prevent concentrations in the solid being taken into nonphysical regimes and also ensure that, if the solid is nearly fully depleted (or filled) with lithium, that the transfer current is not artificially forced to zero. Thus a Butler-Volmer relation of this form allows the flux to move the system away from these depleted (and filled) states. This gives a model which is capable of describing battery charge from a completely depleted state or discharge from a fully charged state. The simplest functions i 0 and U eq that we can use in the Butler-Volmer Equation (3.3) that satisfy the above conditions, have the form 16) and in this instance, according to (3.6) and (3.7),  2012))) take a constant value for i 0 and some of them (cf. Smith & Wang (2006b)) claim that it exhibits modest dependency on electrolyte and solid surface concentration. Although this can be valid for appropriate particular constrained cases, in a general situation this is inappropriate since i 0 should vary near the extreme cases when the battery is either close to being fully charged or fully discharged.

Practical methods of fitting data
When fitting the various functions characterising the Butler-Volmer equation (3.3) to the data it is necessary to include the singular behaviours as c → 0, c s → 0 and c s → c s,max as given previously. A sensible approach is to make non-singular modifications to the simple model (3.17) as follows. There is usually very little accurate data on the exchange current i 0 and so (3.17)-(3.18) is typically the form that is taken. For the OCP one G. W. Richardson et al. approach is to take where f (·; a 1 , a 2 , · · · , a N ) is some family of suitably bounded functions, and g is a similar function allowing for changes in c, while k c , k a , a 1 , · · · a N , b 1 , · · · b M is a set of parameters to be chosen by fitting experimental data, with k, k c , k a satisfying (3.18). Note this could also be expressed in terms of the non-dimensional lithium stoichiometry y = c s /c s,max For example, using the function f proposed in Weng et al. (2014) and a polynomial function for changes in c, the OCP can be fitted to (3.20) For practical particular cases, because there is usually very little data related to variations with c, it may be appropriate to take b i ≡ 0, so that the only c dependency is in the logarithm and this will still ensure the extreme behaviour near c = 0 is physical.

Properties of common electrode materials
The depictions of lithium transport in the active material of the electrode (electrode particles) in §3.1.1, and of the charge transfer reaction between the active material and the electrolyte in §3.1.2, although commonly employed, are an oversimplification of the true behaviour of these materials. In order to highlight some of the nuances of modelling these materials, we briefly review a small part of the copious literature, focusing on some commonly used materials.
The standard negative electrode (anode) material used in commercial lithium-ion batteries is graphitic carbon which alloys with lithium to form Li y C 6 , see, Persson et al. (2010). Other materials (such as silicon) are currently being developed with the aim of supercedeing graphite in this role in the future but none has reached the stage of commercialisation. In contrast to the situation for negative electrodes, there are a wide range of positive electrode materials currently used in commercial batteries; these are reviewed in Julien et al. (2014). Julien et al. (2014) notes that most electrode materials fall into three categories, based on the lithium ion diffusion pathways in the material. In spinel materials lithium transport is three-dimensional, in layered materials it is predominantly two-dimensional and in olivines it is predominantly one-dimensional. Commonly used commercial positive electrode materials include LiMn 2 O 4 (often referred to as LMO) which is a spinel, NMC (mentioned above) which is a layered material and LiFePO 4 (often referred to as LFP) which is an olivine. Li y C 6 (the standard negative electrode material) has a layered structure. The dimensionality of lithium transport in layered (2-d) and olivine (1-d) materials suggests that use of an isotropic model of lithium transport, such as (3.2), is inadequate and should, at the very least, be generalised to an anisotropic diffusion model capable of capturing, for example, the dependence of transport properties on alignment with the crystal axes. Nevertheless an isotropic transport model is frequently used even when it might seem inappropriate. For example, Srinivasan & Newman (2004b) and Arora et al. (2000), Doyle et al. (1996) use isotropic diffusion models for lithium transport in LFP (an olivine) electrode particles and graphitic carbon (a layered material) electrode particles, respectively. There may, however, be good reasons for doing this; for example electrode particles are rarely formed from single crystals and lithium transport in conglomerate particles, formed from randomly oriented crystals, might reasonably be expected to appear isotropic on the lengthscales of interest.
3.3.1 Lithium transport in Li y C 6 (negative electrode material).
Much of the early modelling work on lithium-ion batteries, such as Arora et al. (2000), Doyle et al. (1996), Fuller et al. (1994b, modelled lithium transport within Li y C 6 electrode particles by a linear diffusion equation (3.2). However both Takami et al. (1995), Verbrugge & Koch (2003) and Krachkovskiy et al. (2018) suggest a very strong dependence of solid-state lithium diffusion coefficient D s (c s ) with lithium concentration c s in Li y C 6 particles, with a range of variation of up to about two orders of magnitude, depending upon the exact form of carbon used. In both sets of experiments the size of the carbon particles used was around 10µm and diffusion decreased markedly as lithium concentration c s was increased. To complicate matters further Li y C 6 is known to exhibit (at least) three phases as lithium stoichiometry y is increased. The presence of these phases can be seen inferred from colour changes to the electrode particles (dark blue-low lithium, red-intermediate lithium and gold-high lithium). In Harris et al. (2010) optical microscopy measurements are used to characterise the phase transitions occurring (with increasing lithiation) within a Li y C 6 half-cell anode subject to uniform charging. This shows different phases co-existing (in distinct graphite electrode particles) at different positions in the anode. However, the relevance of these results to commercial devices should perhaps not be overstated, because the width of the negative electrode used in Harris et al. (2010) is particularly large (around 800µm) compare to the standard electrode size in commercial devices (around 100µm). For this reason, the charging process observed in Harris et al. (2010) is likely to be limited by lithium diffusion within the electrolyte, as the electrode particles deplete the surrounding electrolyte of lithium ions. Thomas-Alyea et al. (2017) have also applied optical microscopy to graphite electrodes and observed considerable spatial nonuniformity even after the electrode was left quiescent for an extended period. They were able to predict such states by employing a Cahn-Hilliard phase field model which will be discussed further in §3.4.
Graphite reacts with the electrolyte, consuming lithium ions, to form a thin layer of solid material on the graphite which is referred to as the solid-electrolyte interphase (SEI) layer and, as noted in Bruce et al. (2008), is essential for maintaining the structural integrity of the electrode particles. However, the consumption of electrolyte by this reaction means that graphite electrode particle size cannot be reduced to the nanoscale (in an attempt to improve the charge/discharge rate of the battery) without severely compromising battery capacity Bruce et al. (2008); some lithium makes up the SEI layer and can no longer participate in the useful reactions that store charge. This SEI layer also forms a barrier to lithium ion (and current) transfer between the electrolyte and electrode which may have a significant bearing on electrode performance and has thus been incorporated into some models, for example Srinivasan & Newman (2004a). Diffusion of lithium along the graphene sheets of pure single crystal graphite is extremely fast Persson et al. (2010) (diffusion coefficient of the order of 10 −7 − 10 −6 cm 2 s −1 ), so that, even in electrodes comprised of quite large electrode particles (∼ 100µm) discharged (or charged) at very high rates, it should not significantly affect cell performance. However, Persson et al. (2010) demonstrates that diffusion perpendicular to the graphene sheets and along grain boundaries is many orders of magnitude slower (diffusion coefficient of the order of 10 −11 cm 2 s −1 ) and uses this to infer that this high degree of anisotropy can be used to explain the widely disparate measurements of diffusivity reported in polycrystalline graphite.

Lithium transport in NMC and LMO (positive electrode materials).
Lithium diffusion in NMC Wu et al. (2012) is highly nonlinear so that D s (c s ) decreasing by about two orders of magnitude as lithium concentration within the material increases. Furthermore NMC can be charged and discharged at high rates and the OCP is smooth without the stepped plateau features that usually characterise phase transitions. In contrast, Julien et al. (2014) notes that LMO undergoes a number of phase transitions as it charges and discharges, which are associated with plateaus in its OCP curve. Diffusivity of lithium in single crystals of LMO is about an order of magnitude lower than in multicrystalline particles Das et al. (2005) suggesting that grain boundaries form an easy pathway for lithium diffusion. Furthermore LMO has the disadvantage of capacity loss and fade after repeated cell cycling. This capacity fade has been ascribed, by Das et al. (2005), to the formation of a SEI layer, and consequent loss of lithium mobility.

Bruce et al. (2008) point out that intercalation in LFP involves a phase transition between
FePO 4 and LiFePO 4 , which is reflected in its flat OCP curve (see Figure 3(b)). Kang & Ceder (2009) note that the transport of lithium is dominated by transport along channels in particular crystalline directions, the b-direction, and in single crystal nanoparticles is extremely rapid, so fast indeed that it is doubtful that lithium intercalation in LFP single crystal nanoparticles will ever limit battery performance. This point is clearly made by Johns et al. (2009), who demonstrate that discharge in a half cell nanoparticulate LFP cathode is limited by conduction and transport in the electrolyte. In larger LFP electrode particles, Jugović & Uskoković (2009) point out that performance is significantly impaired because lithium ion transport along the b-direction channels is easily obstructed by grain boundaries and crystal defects. This gives rise to an apparent lithium diffusivity in LFP that decreases sharply as the size of the particle increases, see Malik et al. (2010). Standard models of this material include the so-called shrinking core model, presented in Srinivasan & Newman (2004b), which attempts to capture the phase transition by using a one-dimensional free-boundary model of the phase transition (akin to a Stefan model, see e.g. Rubinšteȋn (2000)) in a spherically symmetric electrode particle. This approach is used in a 3D-scale model, that captures agglomeration of LFP nanocrystals into agglomerate particles, by Dargaville & Farrell (2010) who show good agreement with experimental discharge curves for a wide range of discharge rates. However the shrinking core model is known to predict distributions of the two-phases that are not observed in practice and it is also not easy to implement in a form that allows numerous charge-discharge cycle because of the appearance of multiple free-boundaries, as noted by Farkhondeh & Delacourt (2011). A simpler alternative, suggested in Farkhondeh & Delacourt (2011), is to model lithium transport within LFP particles by a phenomenological nonlinear diffusivity D s (c s ) and this appear to fit discharge data well.

An approach based on Cahn-Hilliard equations of phase separation
The lack of an entirely satisfactory theory of lithium transport in electrode materials that exhibit phase transitions (such as graphite and LFP) has recently led to an alternative, and more fundamental approach in which phase separation with the electrode material is modelled using a Cahn-Hilliard equation. The first use of this approach in this context was by Han et al. (2004), who used it to simulate a generic two-phase material. Subsequently Bai et al. (2011), Cogswell & Bazant (2013, Zeng & Bazant (2014) and Singh et al. (2008) applied this method to LFP using it to study phase separation (and its suppression) in LFP nanoparticles. Both Ferguson & Bazant (2012) and Dargaville & Farrell (2013) have incorporated a Cahn-Hilliard based phase-field description of lithium transport within LFP electrode particles into a porous electrode model. Notably Dargaville & Farrell (2013) compare their results to experimental discharge curves over a wide range of discharge rates, but are unable to obtain a particularly good match to data. In Zeng & Bazant (2014) it is observed that the model is sufficient to capture transitions from solid-solution radial diffusion to two-phase shrinking-core dynamics. In Ferguson & Bazant (2014) fit to data from both LFP and graphite half cells, at very slow discharge rates, with some degree of success (particularly in predicting the positions of the phase transitions across the graphite electrode). One remarkable feature of this work is that it predicts the observed steps in the OCP curves, as a consequence of the phase transitions rather than having to fit a stepped OCP to a potential function as in the standard Newman type model. However these Cahn-Hilliard type models are probably only directly applicable to small single crystal electrode particles, because of the extra physics required to model, for example, obstruction of lithium transport by grain boundaries and defects in larger particles. As mentioned previously in practical applications where the crystals are very small they usually do not limit battery discharge and so accurately capturing their internal transport may be of secondar importance in predicting cell-level features.

A coupled device scale model
The aim of this section is to discuss how macroscopic device-scale equations can be systematically derived from a model of the electrolyte surrounding the electrode particles, the geometry of the electrode particles and a description of the electrolyte reactions G. W. Richardson et al. taking place on the surface of the electrode particles. In this we follow the work of Richardson et al. (2012) who derived the device scale equations for an ideal (dilute) electrolyte and Ciucci & Lai (2011) who derived the device scale equations from a model of a moderately-concentrated electrolyte. We remark that the moderately-concentrated electrolyte model used by Ciucci & Lai (2011) predicts that electrolyte conductivity κ(c) is proportional to the ionic concentration c (as it would for an ideal solution) and so is incapable of adequately describing electrolytes at the concentrations typically occurring in a commercial lithium ion cell. Here we shall extend the work in Richardson et al. (2012) to the moderately-concentrated solution model described in §2.2 and which is applicable to most battery electrolytes.

The microscopic model
The purpose of this section is to set out the equations and boundary conditions of a detailed microscopic model of the battery electrode, including both lithium transport and current flow through the electrode particles and the electrolyte. A portion of a typical electrode geometry is illustrated in Figure 4(a), in which a periodic array of electrode particles occupying regionΩ per (here they have ellipsoidal shape as a possible example) is surrounded by the electrolyte, which occupies the regionV per , and the interface is ∂Ω per . As will be discussed further in §4.2.1 we will assume that the volume occupied by the binder and conductive filler is negligibly small so that the electrode particles and electrolyte completely fill the electrode.
Charge transport in the electrolyte is described by the equations (2.66)-(2.68). At the interface with an electrode particle the transfer current density is given by the Butler-Volmer relation (3.3) and this can be equated to the current flowing into the electrolyte via the boundary condition j · N | ∂Ωper = j tr (c, ϕ, c s , φ s )| ∂Ωper , (4.1) where N is the unit outward normal to the interface (it points into the electrolyte region). A boundary condition for (2.66), the equation for conservation of lithium ions, is provided by noting that all charge transfer across this surface takes place via the motion of lithium ions. Hence, we have the following condition on q p , the flux of positively charged lithium ions, Note that (4.1) and (4.2) imply that the flux of negative ions is zero with q n ·N | ∂Ωper = 0, as required physically.
In each individual electrode particleΩ per , a diffusion equation is solved for lithium concentration in the active material. Generalising (3.2) to an arbitrary shaped particle and allowing for nonlinear diffusion gives

Homogenising the equations in a porous electrode
Here we homogenise the moderately-concentrated electrolyte equations (2.66)-(2.68), with boundary conditions (4.1)-(4.2), over a porous electrode formed by an array of electrode particles permeated by the electrolyte. In order to do this we assume that the electrode can be subdivided into an array of cells over which the electrode structure is locally periodic; that is, the structure inside neighbouring cells is virtually identical but may differ significantly between cells separated on the macroscopic lengthscale. 4 An example of cell microstructure, around an array of ellipsoidal electrode particles, is illustrated in Figure 4. To average the problem using homogenisation we introduce a variablê x to indicate position in the microscopic cell and another variable x to indicate macroscopic position in the entire electrode. Since a very similar analysis has been conducted in Richardson et al. (2012) for a dilute electrolyte we omit the details of the analysis here and merely write down the results (in dimensional form). These consist of macroscopic equations for the lithium ion concentration c and the electrolyte potential (measured with respect to a lithium electrode) ϕ that are formulated in terms of the microscopically volume averaged lithium ion flux q p , the microscopically volume averaged current density j and the microscopically surface averaged transfer current densityj tr . These averaged quantities are formally defined in terms of integrals over the microscopic cells as follows: where the integrals are over the microscopic variable and the averaged functions depend only on the macroscopic space variable x and time t. The quantities |V per |(x), |Ω per |(x) and |S ∂Ωper |(x) are defined by and respectively give the volume of electrolyte, the volume of electrode particle and the surface area of electrode particle in the microscopic cell. The macroscopic homogenised electrolyte equations then take the form v ∂c ∂t where the final term in the current equation (4.9) takes the standard form 2(1−t 0 + )(RT /F c)∇c when the activity is that for an ideal solution (i.e. a e (c) = c/c T ). It is interesting to note that the equations for the electrolyte after homogenisation, (4.8)-(4.9), are very similar in nature to the original equations for the pure electrolyte, (2.66)-(2.68), except that there are now source terms in the conservation equations corresponding to the transfer current from the electrodes. Here v is the volume fraction of the electrolyte defined by and B is the dimensionless permeability tensor whose nine components are defined by the relations for i = 1, 2, 3 and j = 1, 2, 3, (4.10) in which the three characteristic functions χ (j) (j = 1, 2, 3) are solutions to the local cell problems∇ 2 χ (j) = 0 inV per , ∇χ (j) · n| ∂Ωper = e j · n| ∂Ωper , χ (j) periodic inx onV per , where e j is a basis vector in thex j -direction and n is the unit outward normal (pointing fromΩ per intoV per ) to the surface ∂Ω per . We note also here the possibility of using this type of homogenisation technique in conjunction with microscale three-dimensional image data obtained from real battery electrodes in order to obtain more realistic representations of the geometric parameters v , B ij and b et , as discussed for example in Gully et al. (2014) and Foster et al. (2015). In many papers the tensor B is taken to be a constant times the unit tensor, which correspond to a highly symmetric set of particles, such as spherical particles on a regular lattice, and proves to be quite a reasonable model.
The homogenised electrolyte equations (4.8)-(4.9) must be solved in conjunction with the macroscopic equations for the current flow through the solid part of the electrode. These can be obtained by using a constitutive law for current flow in the electrode matrix (3.1) with a current conservation equation that accounts for transfer of charge from the electrode matrix into the electrolyte (4.12) The system of macroscopic equations (4.8)-(4.9) and (4.12) require to be solved along with the microscopic lithium transport equations (4.3) and (4.4), at each point in macroscopic space, in order to determine c s | ∂Ωper , which is required to obtain the transfer current j tr (φ s − ϕ, c s , c) (and hencej tr as given in (4.6)). Note that where the electrode particles are spherical (and isotropic) c s | ∂Ωper is uniform over the particle surface and is thus just a function of the macroscopic variables. It follows therefore that j tr also just a function of the macroscopic variables and, as a consequence of the averaging equation (4.6) it follows that for j tr = j tr for spherically symmetric electrode particles.
(4.13) 4.2.1 Remarks on the role of binder.
In the discussion above we have assumed that the electrode was formed solely from electrode particles bathed in electrolyte. While this is often a reasonable description of research cells, many commercial devices also incorporate a significant volume fraction of polymer binder material that acts both to enhance the structural integrity of the device and, in combination with a conductivity enhancer (such as carbon black), maintain good electrical contact between electrode particles (these are often poor conductors). Three dimensional images of typical commercial electrodes using focused ion beam in combination with scanning electron microscopy (FIB-SEM) can be found in Gully et al. (2014), Foster et al. (2015 and Liu et al. (2016). These show a porous binder material filling almost all the space between electrode particles with the exception of some linear features around the electrode particles where it appears that the binder has become delaminated from the electrode particles. 5 The porosity and pore size of the binder materials varies significantly between different electrode types. Typical pore sizes are usually in the range 10-500nm, much smaller than typical electrode particle sizes which are usually at least micron sized. To account for these effects within an homogenisation approach the analysis could be modified in two possible ways. Firstly it could be performed directly on a microstructure in which all three constituents (electrode particle, electrolyte and binder) 32 G. W. Richardson et al. are resolved by the cell problem and would follow a very similar pattern to that described above in which the binder were treated as part ofΩ per with an interface with the electrolyte on which the transfer current density j tr ≡ 0. Secondly, because obtaining a good representation of the pore geometry in the binder is challenging, even using modern high-performance microscopy, see Liu et al. (2016), it is probably better to treat the electrolyte permeated nanoporous binder as a single electrolyte material (albeit it one with reduced electrolyte volume fraction, electrolyte diffusivity and conductivity) and homogenise over the electrode particles and this composite material. Indeed this second approach is a standard way of treating such two phase (electrolyte/binder) materials in the literature which are often termed porous solid polymer electrolytes (see, e.g. Miao et al. (2008)).
As mentioned previously a common approach in the literature to modelling practical batteries is to use the so-called pseudo-2d model where the behaviour on the macroscale is one-dimensional (with spatial position denoted by x) and transport of lithium on the microscale takes place within spherical particles, and is thus also one-dimensional taking place in the particles' radial direction (with radial position denoted by r). These assumptions are tantamount to assuming that the averaged macroscopic vector-valued currents and fluxes from §4.1 and 4.2 are only non-zero in the x-direction. Henceforth we will replace these vector quantities with scalar counterparts. Here we describe the system of equations that arise for such a situation exploiting the homogenisation results described previously in (4.8)-(4.12) and including the typical units. As alluded to above, x denotes position across the cell in the direction perpendicular to the current collectors which are positioned at x = L 1 and x = L 4 respectively, so that x ∈ (L 1 , L 4 ). The cell is subdivided into three regions (as illustrated in figure 1) with the negative electrode occupying the region x ∈ (L 1 , L 2 ), the separator occupying the region x ∈ (L 2 , L 3 ) and the positive electrode occupying the region x ∈ (L 3 , L 4 ). In the electrolyte, which permeates the whole cell (i.e. x ∈ (L 1 , L 4 )), we seek to determine the Li concentration c(x, t) (mol m −3 ), the electric potential (measured with respect to a reference lithium electrode) ϕ(x, t) (V) and the ionic current density j (x, t) (A m −2 ). In the negative electrode (x ∈ (L 1 , L 2 )) we seek solutions for the solid phase potential φ (a) s (x, t) and the solid phase current j (a) s (x, t). In addition we seek to determine the lithium distribution c (a) s (r, x, t) within the (negative) electrode particles (of radius R (a) (x)) as a function of position r ∈ [0, R (a) (x)] within the particle and the position x of the particle within the electrode. Similarly in the positive electrode (x ∈ (L 3 , L 4 )) we seek solutions for the solid phase potential φ  s (x, t) and also the lithium distribution c (c) s (r, x, t) within the (positive) electrode particles (of radius R (c) (x)) as a function both of position r ∈ [0, R (c) (x)] within the particle and the position x of the particle within the electrode. A sketch of the device geometry as well as an illustration of the domains of definition of the dependent variables is shown in figure 1. Throughout what follows we assume that the transport of lithium within both negative and positive electrode particles occurs through nonlinear isotropic diffusion, though as discussed in §3.3 this is not the only possibility. Furthermore since the electrode particles are spherical we can make use of the simplification (4.13) in order to writej tr = j tr . The assumptions that the particles are spherical and that their radii vary slowly, i.e. that their size depends on macroscopic position but neighbouring particles are almost the same size, gives rise to the following relationships where n(x) is the number density of electrode particles and, owing to our assumption that volume fraction of binder is negligible, 1 − v is the volume fraction of electrode particles.
Here we set out the pseudo 2d-model, which follows from the homogenisation described in §4.2 and describes the performance of a cell at constant temperature. In this model we denote variables and parameters relating to the negative electrode (or anode) by the superscript (a) and variables and parameters relating to the positive electrode (or cathode) by the superscript (c). The transport equations for the electrolyte, obtained from (4.8)-(4.9), are where the volumetric current source term S(x, t) is given by (4.16) and where the electrolyte volume fraction v and the B 11 component of the permeability tensor are evaluated appropriately in each of the three regions. As discussed in §4.2, B 11 can be computed by solving the appropriate cell problems. However, a common approach is to instead estimate its value using B 11 = p v where p is the Bruggeman porosity exponent (a nondimensional constant), which is commonly taken to be p = 1.5, see Bruggeman (1935), Gully et al. (2014) and Gupta et al. (2011). We note also the work of Shen & Chen (2007) who discuss some alternative estimation methods beyond the Bruggeman approximation. Boundary conditions on the electrolyte equations are enforced by the requirements that there is no flow of electrolyte current or flux of lithium ions into the current collectors at x = L 1 and x = L 4 and are (4.17) In the negative electrode matrix conservation of current and Ohm's Law, as given by (4.12), are described by Here I is the current flowing into the cell and A is the cell's area. In the same region, lithium transport (as described in (4.3) The transfer currents densities j s | r=R (c) (x) , c) that describe the flow of current out of the anode and cathode particles, respectively, and which act to couple together the lithium transport problems in the electrolyte to those in the electrode particles are typically determined by the Butler-Volmer condition as discussed in §3.2. The existence and uniqueness of solution of the system (4.15)-(4.27) has been proved in Díaz et al. (2019). We point out that, in some works in the literature (see, for instance, Gomadam et al. (2002), Smith & Wang (2006a), Smith & Wang (2006b), Smith et al. (2007), and Kim et al. (2012)), authors present the following (incorrect) boundary conditions at x = L 1 and x = L 4 . We note in passing that if these incorrect conditions are used then it can be proved (see Ramos (2016)) that the corresponding system of boundary value problems does not have any solution unless I(t) ≡ 0.
Once the model described above has been solved, we can estimate the state of the charge of the negative electrode SOC (a) (t) and of the positive electrode SOC (c) (t), and the cell voltage V (t), at time t, by computing where there is a constant resistance R f which accounts for the potential dropped in the current collectors, but we remark that this is often negligible in practice. We note some ambiguity about how the state of charge of an electrode or cell should be defined. The definition given here might be referred to as the state of charge measured with respect to the theoretical maximum capacity. In the experimental literature it is common to define the state of charge with reference to the capacity measured from a cell under a low-rate (dis)charge, see e.g. Johns et al. (2009), Barai et al. (2015).  (2015) where a series of experiments were conducted on a high energy pouch cell produced by Kokam; the values used here are summarised in Table 1. Figure 5 shows the discharge curve, and internal concentration and potential profiles during a relatively low-rate usage where a current demand of 0.13A is applied for 4000s and the cell is the immediately recharged at the same rate until it reaches a cut-off voltage of 4.2V. We observe that under these conditions concentration and potential gradients in the electrolyte are relatively modest. Likewise, the concentration is through the radius of the electrode particles is almost uniform; a consequence of the relatively large diffusivity in LNC. The largest gradients are observed internal to the anode particles and although these are not large enough to hamper the initial discharge, it is the inability of the graphite to transport intercalated Li from its surface into its interior that ultimately causes the recharging process to be interupted. At the final snapshot in time in panel (d) we observe that the concentration on the surface of the graphite particle has reached its maximum and therefore intercalation cannot proceed further at this location despite their being available space to accomodate Li in the particle's interior. Figure 6 shows the same undergoing a similar discharging and subsequent recharging protocol, but at a more aggresive demand of 1.3A for a shorter time of 400s, followed by a subsequent aggresive recharging again at 1.3A. Note that this faster discharge supplies the same amount of charge to the external circuit as the slower protocol but in a time window 10 times smaller. At the increased rate we observe that gradients in the electrolyte are much more pronounced, and in fact they are sufficiently large that the deep regions of the cathode approach depletion during discharge. The same is true in the anode during recharging. This larger polarisation contributes to a diminished cell voltage during discharge and we can observe that at the deepest discharge state the voltage has dropped to 2.5V, which is markedly lower than the 3.5V attained during the slower protocol despite the devices supplying the same amount of charge. There are now also noticeable concentration gradients within the LNC electrode particles and gradients in the graphite particles are very high. Once again, it is the graphite which ultimately causes recharging to terminate because the surface of the graphite particles becomes saturated and the recharging can only for around 140s.

Numerical solutions to the pseudo 2d-model
The solutions shown in Figures 5 and 6 were determined numerically using an in-house ultra-fast and robust solver called DandeLiion, the details of which will be described in a forthcoming paper, see Korotkin et al. (2020). Whilst is is beyond the scope of this work to reiterate, in detail, the workings of the numerical methods used it is pertinent to outline the approach and highlight some of difficulties in solving (4.15)-(4.27) numerically. First, a spatial mesh is defined. We introduce N (a) grid points across the anode for x ∈ (L 1 , L 2 ), N (s) across the separator for x ∈ (L 2 , L 3 ) and N (c) across the cathode for x ∈ (L 3 , L 4 ). At each point in the anode and cathode a microscopic transport problem must be solved and so at each of the N (a) + N (c) grid points a further M grid points need to be introduced to on which to discretise the microscopic transport equations within the electrode particles. Consequently, the complete discrete geometry is comprised of (N (a) + N (c) ) × M + N (a) + N (s) + N (c) grid points. Second, a suitable approximation (e.g. finite volumes or finite elements) can be used to remove the spatial derivatives and reduce the problem to a large system of coupled differential-algebraic equations (DAEs). The algebraic equations arise largely from the elliptic PDEs, e.g. those for the electron conduction in the solid (4.18), whereas the ordinary differential equations arise from the parabolic PDEs, e.g. those for transport in the electrode particles (4.20)-(4.22). We note the importance of using a spatial discretisation method which is conservative; if such a method is not used, on repeated cell cycling, the total amount of lithium within the system changes markedly and introduces significant errors. It is for this reason that many approaches based on finite difference approximations are not recommended. Third, a scheme for timestepping the system of DAEs must be found. The DandeLiion software uses uses a selection of implicit backward differentiation formula methods, of orders 1-6, and also offers adaptive time stepping. The choices of timestepping methods is restricted, in comparison to those that can be used for pure ODE systems, because of the additional constraints imposed by the algebraic equations.
Implementing the steps outlined above and implementing in C++ gives rise to a numerical scheme for solving the pseudo 2d-model in a very short time on a standard Table 1. The parameter values used to carry out the simulations shown in §5. These are largely based on the work of Ecker, Tran, Dechent, Käbitz, Warnecke & Sauer (2015) and Ecker, Käbitz, Laresgoiti & Sauer (2015). The functions used for the electrode conductivities were fitted to data in Ecker, Tran, Dechent, Käbitz, Warnecke & Sauer (2015) and Ecker, Käbitz, Laresgoiti & Sauer (2015) and the functions themselves are given in Korotkin et al. (2020).   Figures 5 and 6 each took around 5 seconds to run and were performed on a discretised geometry comprised of 20,300 grid points with 100 spatial points inside the anode, separator, cathode and each of the electrode particles.
We have reviewed the existing modelling of charge transport models of lithium ion batteries at the cell scale. This includes a description of ionic motion in the electrolyte in both the dilute and the more practically relevant moderately-concentrated regimes. We resolve a common source of confusion in electrolyte modelling in the literature which arises from the definition of the electric potential. In the electrochemical literature this is usually chosen to be the potential measured with respect to metallic lithium electrode, in contrast the standard definition used in the physics community it is with respect to a vacuum at infinity. Crucially this choice of potential affects the coefficients in the electrolyte transport equations. We have also examined the Butler-Volmer relation describing reaction at the interface between the electrolyte and the solid electrode particles and demonstrated that these should have a particular functional form in order to avoid nonphysical predictions. The dependency in the Butler-Volmer relation should not only account for the solid electrode particle becoming completely intercalated or deintercalated but allow for cases where the electrolyte concentration gets very low. The specific behaviour of various common solid materials used in electrodes have been considered including the dominant mechanisms for lithium transport and the possible modelling approaches that can be used. The problem of up-scaling the models from the microscopic (a single electrode particle) to the macroscale (the whole cell) has been considered and the appropriate approximations discussed that allow the models to account for moderatelyconcentrated electrolyte behaviour reviewed. In addition numerical solution to the Newman model is discussed and some representative realistic solutions to the resulting pseudo 2-dimensional model have been presented and discussed. It is our hope that this work will prove a useful guide to people who are new to this topic allowing them to develop an appreciation of this highly fertile and technologically important area of research.     Table 1. Thick curves indiciate profiles at the beginning of the (dis)charge stages and the different snapshots are taken every 500s. Panels (a)-(c) show the cell potential, electrolyte potential and electrolyte concentrations respectively, and panels (d) and (e) indicate profiles within an anode and cathode particle both of which are located half way through the thickness of their respective electrodes.  Table 1. Thick curves indiciate profiles at the beginning of the (dis)charge stages and the different snapshots are taken every 50s. Panels (a)-(c) show the cell potential, electrolyte potential and electrolyte concentrations respectively, and panels (d) and (e) indicate profiles within an anode and cathode particle both of which are located half way through the thickness of their respective electrodes.