A Discussion of the Cell Voltage during Discharge of an Intercalation Electrode for Various C-Rates Based on Non-Equilibrium Thermodynamics and Numerical Simulations

For the electrolyte we consider exclusively the material model^9–11 of an incompressible liquid electrolyte accounting for solvation effects, i.e.

with mole fraction

molar concentration n_α, and total molar concentration of the mixture (with respect to the number of mixing particles⁹)

In (4) T denotes temperature, k_B  the Boltzmann constant, g^R_α denotes the reference molar Gibbs free energy (or chemical potential of the pure substance), p^R the reference pressure and v^R_α the partial molar volume of constituent α in the mixture. Throughout this manuscript we assume an isothermal temperature of .

Note that denotes the number of free solvent molecules, whereas and are the densities of the solvated ions. This is crucial for various aspects of the thermodynamic model, and we refer to Refs. 9,10,12,13 for details. Overall, the material model for the electrolyte corresponds to an incompressible mixture with solvation effects. We assume further

whereby the incompressibility constraint^9–11 implies also a conservation of mass, i.e.

The molar volume of the solvent is related to the mole density of the pure solvent as

Note further that the partial molar volumes v^R_α and the molar masses m_α of the cation and anion are related to the solvation number and , respectively.

We assume that the partial molar volume of the ionic species is mainly determined by the solvation shell, which seems reasonable for large solvents like DMC in comparison to the small ions like Li⁺. We proceed thus with the assumption

For the active particle, we consider an extension of a classical lattice mixture model^14–21 which accounts for occupation numbers as well as a Redlich–Kister type enthalpy term^22,23 for the intercalation material Li_y(Ni_1/3Mn_1/3Co_1/3)O₂ (NMC). We refer to Ref. 24 for a detailed discussion and derivation based on a free energy . The chemical potential of intercalated lithium is derived as

with

and mole fraction

of intercalated cations in the active phase. The number density of lattice sites is constant, which corresponds to an incompressible lattice, and the enthalpy parameter . Note that entails a phase separation²⁰ and requires an additional term in the chemical potential. However, we assume throughout this work that no phase separation occurs, whereby in diffusional equilibrium of the intercalation phase the concentration is homogeneous. An extension of this discussion toward phase separating materials will given in a subsequent work.

For the electrons we consider^9,25

and for the lattice ions

where is the molar volume of the lattice ions, the partial pressure and the constant molar Gibbs energy. The material functions of the active intercalation electrode is essentially model an incompressible solid with a sub-lattice for the intercalated cations .

The explicit surface chemical potentials

are not required throughout this work since we will assume that the double layer is in equilibrium and that the double layer capacity (and thus also adsorption), is negligible for the sake of this work. However, we refer to Ref. 9 for the explicit functions of and the surface free energy of a surface lattice mixture with solvation effects.

The electroneutrality condition of , and can be obtained by an asymptotic expansion of the balance equations in the electrochemical double layer at the respective surface Σ. We only briefly recapture the central conclusions and refer to Refs. 3–5,9,26 for details on the modeling, validation and the asymptotics. Most importantly, we have that

• the double layer is in thermodynamic equilibrium, i.e. ∇μ_α + e₀z_α∇φ = 0 in and , where e₀ is the elementary charge, z_α the charge number and φ the electrostatic potential
• there exists a potential drop between the active particle surface Σ and the hyper-surface outside of the respective space charge layers which is denoted by
  
  where is the electrostatic potential right outside the space charge layer in the electrolyte or the active particle, respectively, and the (continuous) potential at the surface Σ^a. The whole potential drop across the double layer at is denoted by
  
• the chemical potential at the surface can be pulled back through the double layer, i.e.
• the condition entails that the potential drop is constant (with respect to some applied voltage) and determined by
  
• the charge density in the electrolyte vanishes and that for monovalent electrolytes the cation mole fraction (or number density) is equal to the anion mole fraction, i.e.
  
• in the active phase the electroneutrality entails
  
  whereby we abbreviate
  
  which is basically the Fermi energy of the solid material.

In the electrolyte we have two balance equations determining the concentration (or the mole fraction ) and the electrostatic potential in the electrolyte,^27–32 i.e.

with (dimensionless) thermodynamic factor

where

is the thermodynamic driving force for diffusion,¹¹ and the chemical diffusion coefficient, the cation transfer number, and the molar conductivity.

Note that we assumed and which yields the representation 26. Note further that the total number density in the electrolyte writes as

which is determined from the incompressibility constraint 8

and the electrolyte concentration in terms of as

If we consider a simple Nernst–Planck-flux relation for the cation and anion fluxes,^11,33 respectively, i.e.

with constant diffusion coefficients for the anion and for the cation, we obtain (in the electroneutral electrolyte)

Note, however, for general Maxwell-Stefan type diffusion^29–32,34 or cross-diffusion coefficients^7,24,35 in the cation and anion fluxes lead to more complex representations of the transport parameters . In general, three of the transport parameters are independent, and , and are related to each other via

Further, depend in general non-linearly on the electrolyte concentration . However, it is sufficient for the sake of this work to assume constant values for the transport parameters , together with relation 33.

In the active particle we have two balance equations determining the concentration (or mole fraction ) and the electrostatic potential in the active particle, i.e.

and (dimensionless) thermodynamic factor

Note that in principle can be dependent on the amount of intercalated ions, i.e. .

We want to investigate the non-equilibrium thermodynamic modeling of the intercalation reaction

Surface thermodynamics dictates that the reaction rate of this process can in general be written as^4,5,13,36,37

with α ∈ [0, 1]. Note that a non-negative function in 38 ensures a non-negative entropy production due to reactions on the surface, i.e. .

The quantity can be considered as surface affinity of the Reaction 37. The surface reaction rate vanishes when the affinity vanishes, which is the actually the thermodynamic equilibrium condition of 37, i.e. .

Since the electrochemical double layer is in equilibrium, we can pull back the surface chemical potentials through the double layer to the respective points (in an asymptotic sense) outside of the double layer, whereby we obtain for the surface affinity

With the material models 4 and 11 we can rewrite the surface affinity as

with

and

with according to 12. Note again that denotes the evaluation of at the interface and that the surface affinity 40 is dependent on the chemical potential (or the mole fraction) evaluated at the interface.

We consider the cell voltage in a half cell with metallic lithium as counter electrode, denoted by and position at (see Fig. 2. The cell voltage in such a cell is

where is the potential drop in the bulk active particle due to the electron transport, is the potential drop across the double layer at the interface between the active particle and the electrolyte, and the bulk potential drop due to cation electric current.

Zoom In Zoom Out Reset image size

We assume that the counter electrode is ideally polarizable,²⁸ whereby the reaction

at the the interface positioned at is in thermodynamic equilibrium and . The equilibrium condition of 45 entails

where is the chemical potential of the metallic lithium.

For the surface affinity 40 we obtain the compact typeface

with

and

For the single intercalation reaction we have the following expression⁴

for the current density i flowing out of the electrode , where is the double layer capacity. Note that the reaction rate is

Since g(x) is a strictly monotone function, we can introduce the inverse of g, i.e. g^{− 1}. For we have and . For values α ≠ 0.5 the inverse function g^{− 1} is only implicitly given, however, can easily be calculated numerically. Fig. 3 display the functions g and g^{− 1} for various values of α. We call g(x) the reaction rate function and g^{− 1} the inverse reaction rate function.

Zoom In Zoom Out Reset image size

Note that in the Tafel approximation Eq. 51 yields^b

The term can be considered as the exchange current density.²⁸

The Onsager coefficient (or the exchange current density ) of the surface reaction 37 could in principle be a function of the surface chemical potentials (or surface concentrations), i.e. or or the surface affinity, i.e. , as long as the condition is ensured.^4,8,26 Note, however, that surface thermodynamics dictates the dependency of on the surface chemical potentials and not the bulk chemical potentials μ_α.

For a general relation we can pull back the surface chemical potentials through the double layer to obtain

Note that this necessarily restricts the functional dependency of on the mole fractions at the interface .

Consider, for example a model , where the exchange current density is dependent on the electrolyte concentration at the interface. This would be, however, thermodynamically inconsistent since the general functional dependency of 54 requires for the electrolyte concentration at the interface

Another commonly used model is a functional dependency of on the concentration of intercalated ions at the interface, i.e. . Since the space charge layer in the active particle is essentially constant (because is constant), we can indeed write

We discuss this aspect as well as various models for in Discussion of the exchange current density section. Meanwhile we assume and proceed the following derivation and the discussion based on this assumption since it turns out to be very reasonable.

At this stage, it is illustrative to discuss the explicit value of the parameters.

• For the electrode geometry we consider for a planar surface of area A and a thickness which yields and . The electrolyte is considered with a thickness of . This corresponds to the cell dimensions of the cell MX-6 in Ref. 2.
• Throughout this work we consider DMC as solvent with and assume for the solvation number . The reference electrolyte concentration is and average amount of electrolyte is and a parameter of the model.
• Average concentrations (or mole fractions) are abbreviated as
  
  for the electrolyte species and
  
  for the amount of intercalated ions in the active phase.
• For the active particle phase we consider Li(Ni_1/3Mn_1/3Co_1/3)O₂ (NMC) whereby
  
  which is simply computed from the density and stoichiometry of the bulk material.³⁸ As parameters for the chemical potential we consider an occupation number of and a Redlich–Kister interaction energy of .²⁴
• The differential capacity has a prescribed value (actually is a function of , but we proceed here with a constant approximation for the sake of simplicity⁹ of about
  
• The electrode capacity Q is
  
  This yields the non-dimensional capacity
  
  which is sometimes also called status of charge (SOC) or depth of discharge (DOD).Note that during discharge of a complete battery the cathode is actually filled up with lithium. In a half cell with metallic lithium as counter electrode, discharge thus actually means filling up the intercalation electrode, here the NMC cathode material. Hence corresponds to a fully charged cathode (i.e. no lithium in the intercalation compound, ) while corresponds to a fully discharged cathode (i.e. the intercalation compound is completely filled with lithium, ).
• From the charge balance 35 of the active particle we can deduce
  
  where I is the current flowing into the intercalation electrode during discharge and Q(t = 0) the initial charge state. For a galvanostatic discharge I > 0 we obtain thus
  
• The C-Rate defines (implicitly) the current at which after h-hours the intercalation cathode is completely filled during galvanostatic discharge. C₁ is thus the rate at which the battery is charged within one hour and commonly abbreviated just as C-rate C, i.e.
  
  We can hence express the current I in multiples of the C-rate, i.e.
  
  which yields
  
  The only parameter for the current density i = I/A is thus C_h.
• For the time t we consider the interval of one discharge cycle, i.e. t ∈ [0, t_end] with
  
  We can thus introduce the non-dimensional time
  
  whereby the capacity rewrites as
  
• For the current density i at the planar electrode we have thus
  

Discussion of the scaling

Consider the non-dimensional voltage

and abbreviate

which yields

with

The parameters and yield

and

whereby

The double layer contribution in Eq. 51 is thus almost negligible whereby 51 reduces to

We consider for the exchange current density the rescaling

This is the crucial decomposition throughout this work and the parameter of the surface reaction rate .

For the current density and the inverse function g^{− 1} we obtain thus with Eq. 48 for the surface affinity the general expression

for the cell voltage E.

The open circuit potential (OCP) is obtained from 81 as

for C_h = 0 (infinite slow discharge), which entails also as well as . Hence we have

For Li_y(Ni_1/3Mn_1/3Co_1/3)O₂³⁹ as intercalation electrode, the two parameters of the chemical potential function are the occupation number and the interaction energy of the Redlich–Kister type enthalpy contribution. This yields an absolute ℓ²-error of and a relative error of 1.860% vs. experimental data of P. Bruce et al.,³⁹ and Fig. 4 shows a comparison to two experimental data sets of measured OCP data.

Zoom In Zoom Out Reset image size

Infinite conductivity within the active particle phase as well as within the electrolyte yields

and infinite fast diffusion in the active particle and the electrolyte entails

Hence is directly related to the capacity via

whereby the cell voltage of BV 1 is

It is a simple algebraic relation between the measured cell voltage E, the C-rate C_h, the capacity Q and the (non-dimensional) exchange current density .

In order to compare the cell voltage E computed in the approximation BV 1 with other approximations, we abbreviate the voltage computed from 89 as E⁽¹⁾. Note that cell voltage 89 is actually independent of the electrolyte. For we obtain the voltage/capacity relation given in Fig. 5 for a variation of C_h from 0 (open circuit potential) to C_h = 100 (extremely fast discharge).

Zoom In Zoom Out Reset image size

Reaction overpotential

We define the reaction overpotential as

which is actually independent of the status of charge or capacity. Measured voltage data would thus allow to determine and the parameter α ∈ (0, 1).

Fig. 6 shows computations of the reaction overpotential η^R for various values of α and as function of C_h.

Zoom In Zoom Out Reset image size

Finite conductivity within active particle phase entails from Eq. 35

Employing the scaling 71 of the current density i, i.e. , as well as the decomposition

yields

The quantity is the specific conductivity of the active particle phase at C-rate of one. For the parameters given in Discussion of the model parameters section computes as

The measured cell voltage is then

which is (yet again) a simple algebraic relation between E, the C-rate C_h, and the capacity . E⁽²⁾ is additionally parametrically dependent on the conductivity .

We define the active phase conductivity overpotential as

Constant diffusion coefficient

Reconsider that we have assumed yet with respect to space in the intercalation particle. In general, however, we have to solve a (here 1D) diffusion equation

with

This yields at the interface some solution

which will also impact the cell voltage

In order to discuss this impact systematically, we apply the following scaling

as well as

which leads to

The dimensionless flux

yields the dimensionless diffusion coefficient

and thus

At the interface we have thus

Overall we may write

with

Note that we can analytically compute from Eq. 11 (see also appendix B1) as

Since the problem 109 is non-linear, a classical separation Ansatz is not meaningful. We proceed thus with solving the problem 109 with 110 numerically with MATLAB and the pdepe() function. The syntax for pdepe() of the problem 109 with 110 is given in appendix B2.

Based on the numerical solution we compute then numerically for various values of C_h and . The (global) capacity is yet .

We assume the same parameters as before, now additionally with two values of the diffusion coefficient , i.e. slow diffusion and fast diffusion , and compute as function of the capacity (or time τ). Fig. 7 shows computations of for various discharge rates and diffusion coefficients in the active particle phase as function of the cell capacity. The angle bisection in black corresponds to the open circuit potential situation, where . For increasing discharge rates, the concentration at the interface is larger than the average concentration in since the evacuation of intercalated ions is delayed by the finite diffusion. This effect becomes even stronger for smaller values of , i.e. slow diffusion in the active particle.

Zoom In Zoom Out Reset image size

The cell voltage E is then computed a posteriori from 100 based on the numerical solution of . Fig. 8 displays the cell voltage for various discharge rates as well as slow () and fast () diffusion in the intercalation phase. Finite diffusion in the active particle has an enormous impact on the cell voltage and changes qualitatively the shape due to the non-linear feedback. This effect is also found experimentally, see Fig. 1, and extremely important since it determines the maximum amount of charge that can be withdrawn from an intercalation electrode.

Zoom In Zoom Out Reset image size

Two important measures serve to discuss the impact of the diffusion coefficient ,

• the cell voltage at 50% discharge, i.e. ,
• and the capacity at the cut off voltage E^off, here with .

Fig. 9 shows numerical computations of and for various values of the C-rate C_h and diffusion coefficients in the range of 10^{− 3} − 10². For slow discharge rates, i.e. C_h < 1 a diffusion coefficient of is sufficient to achieve a voltage of at 50% discharge and capacity of 90% at the the cutoff voltage. However, for higher C-rates, e.g. C_h = 50, the impact of the solid state diffusion becomes enormous, requiring a diffusion coefficient of to discharge the electrode to 50%.

Zoom In Zoom Out Reset image size

Overpotential

The overpotential due to finite diffusion in the active particle phase can be defined as

which computes as

Fig. 10 shows computations of for slow and fast diffusion.

Zoom In Zoom Out Reset image size

Concentration dependent diffusion coefficient

The diffusion coefficient in 34 was yet assumed to be a constant^c with respect to the mole fraction of intercalated ions in the solid phase. This assumption might be inappropriate or over-simplified, and we seek to discuss this hypothesis again on the cell voltage E as function the capacity Q and the C-rate C. For this sake, we consider in the following a concentration dependent diffusion coefficient

where the constant is scaled equally than before, yielding again dimensionless diffusion coefficient . The term models that the diffusion on a lattice reduces when the lattice becomes more occupied (i.e. ).

We obtain hence the non-dimensionalized transport equation

with boundary conditions

This equation system is again solved numerically with MATLAB and the pdepe() function yielding numerical solutions whereby we compute then for various values of C_h and . The (global) capacity is yet . The cell voltage E is then computed a posteriori from 100 based on the numerical solution of .

We assume the same parameters as before and two values of the diffusion coefficient , i.e. slow diffusion and fast diffusion . Based on the numerical solution of 115 and 116 we compute as function of the capacity (or time τ), and subsequent the cell voltage E from 100 based on the numerical solution of .

Fig. 11 displays the cell voltage for various discharge rates computed numerically with a constant (dashed line) and the concentration dependent (solid line) diffusion coefficient according to 114.

Zoom In Zoom Out Reset image size

Expectably, the two models behave similar when the intercalation electrode is empty, i.e. close zero. However, when the electrode gets filled with lithium (i.e. ), the -term of the concentration dependent diffusion coefficient reduces significantly the diffusivity, which affects the cell voltage in a surprisingly non-linear behavior.

Most important, Fig. 11 shows that the concentration dependent diffusion coefficient 114 could explain the decrease of the capacity at the end of discharge, e.g. here at , for increasing discharge current densities (c.f. the experimental data in Fig. 1). Wu et al. state in Ref. 2 that "the electrode capacity measured at the end of discharge decreases with increasing rate, as a result of transport limitations in the electrode". This can be confirmed on the basis of concentration dependent diffusion coefficient with

First note that an infinite fast diffusion in the electrolyte yet entails , whereby the (coupled) transport equation system 23–24 of the electrolyte reduces to

which yields

Employing the scaling 71 of the current density i, i.e. yields

which motivates the decomposition

Here is a constant reference electrolyte concentration, e.g. , and is the corresponding reference conductivity. Hence

whereby the cell voltage is

Hence, finite conductivity in the electrolyte linearly decreases the cell voltage and scales also with the ratio of the electrode width to the electrolyte width, i.e. . The quantity accounts for concentration dependence of the electrolyte conductivity.

Correspondingly we define the electrolyte conductivity overpotential as

The final contribution to the surface reaction is the space dependent electrolyte concentration. We have yet assumed with respect to space, however, in general the (coupled) equation system 23–24 has to be solved.

Note that simplifies the (coupled) equation system 23–24 to

Further, entails at the interface the condition

We assume that the average electrolyte concentration is constant in time, i.e.

which yields at the right boundary the condition

The concept of an ideally polarizable counter-electrode , positioned at , delivers (or consumes) hence exactly the amount of ions in the electrolyte which flow in (or out) of at , with keeping the reaction 45 in thermodynamic equilibrium. The initial value is

where is the prescribed average electrolyte concentration.

We introduce the scalings

with and . This yields for 125

The corresponding non-dimensional boundary conditions at ξ = 0 reads

which introduces (implicitly) the scaling

leaving

The balance Equation 125 then reads

with

Note that the 2 in accounts for the charge of cations and anions. The charge capacity of a electrolyte is

whereby

The dimensionless transference number and . The PDE is solved with MATLAB's pdepe function, and details are given in the appendix A3. We denote the numerical solution of with and emphasize that the capacity is yet . The numerical solutions at the respective boundaries and are

We discuss now briefly the concentration distribution in the electrolyte as function of the C-rate C_h and the diffusion coefficient based on numerical solutions of 125 with boundary conditions 127 and 129.

Fig. 12 displays computations of the electrolytic cation concentration at the interface of the intercalation electrode, i.e. , and at the interface of the counter electrode, i.e. for slow electrolytic diffusion (, left) and fast diffusion (, right) for various values of the C-rate.

Zoom In Zoom Out Reset image size

After a short time the concentration yields a stationary state and Fig. 13 displays the stationary concentration in the electrolyte, again for slow and fast diffusion as well as for various C-rates.

Zoom In Zoom Out Reset image size

Note, however, that the concentration variation of has additionally an impact on the voltage drop . First reconsider that 126 rewrites as

since

We abbreviate

and insert the scaling which yields

The overall cell voltage of BV 5 is then

In order to show the impact of the electrolyte concentration variation on the cell voltage E, assume as well infinite fast diffusion in the active particle phase . This yields

and numerical computations of the cell voltage for slow and fast diffusion are shown in Fig. 14.

Zoom In Zoom Out Reset image size

Overpotential

Due to the (stationary) concentration gradients in the electrolyte (c.f. Fig. 13) we have a diffusional overpotential , which can be defined as

Fig. 15 displays numerical computations of the overpotential for slow and fast diffusion in the electrolyte.

Zoom In Zoom Out Reset image size

The recursive definition of the various overpotentials allows us to write

with one overpotential for each non-equilibrium process, measuring the deviation from the equilibrium of the respective process. This decomposition is hence a useful tool to systematically investigate the contribution of each process in broadly conceived experimental or numerical studies of a cell batch with varying parameters.

Internal resistance

Note that we can also compute the internal resistance R^int of the electrochemical cell via the implicit definition

With we obtain for the specific resistance

with

and

• intercalation reaction resistance
  
• active phase diffusional resistance
  
• active phase conduction resistance
  
• electrolyte diffusional resistance
  
• electrolyte conduction resistance
  

Equation 146 for the general relation for the cell voltage E in a simple, non-porous intercalation electrode. Note, however, that the scalings, discussion and parameter study of Discussion section can be straight forward adapted to porous electrodes.

We provide finally a validation study for NMC with the following set of the non-dimensional parameters

and a concentration dependent diffusion coefficient . In absolute values, these translate to:

• exchange current density ,
• NMC electric conductivity of ,
• lithium diffusion coefficient in NMC ,
• electrolyte conductivity ,
• electrolyte diffusion coefficient

Figure 16 displays the numerical computation of the cell voltage. In comparison to experimental data for a cell of the same dimension (however, neglecting porosity), we obtain a good qualitative and quantitative agreement to Fig. 1.

Zoom In Zoom Out Reset image size

This is especially remarkable since we assumed essentially for all non-equilibrium parameters constant values, i.e. no concentration dependence of the diffusion coefficient , the cation transference number , and the conductivities. In particular we assumed that the exchange current density is also constant, yielding the reasonable results of the last section. Note, however, that it is frequently assumed that the exchange current density is dependent on cation concentration at the interface . We discuss this aspect in the next section and emphasize again that a consistent thermodynamic modeling as well as coupling through the surface reaction rate yields the reasonable results of the last sections. The scaling of all non-equilibrium parameters to the C-rate is quite illustrative for the sake of galvanostatic discharge and especially for the systematic search of the parameters of a specific battery.

The preceding discussion of the cell voltage E was based on the model 38 of the surface reaction rate , i.e.

with In Discussion of the model parameters section we showed that double layer charging effects are negligible under galvanostatic conditions, whereby the measurable current density i is directly related to , i.e. . The surface affinity is related to the cell voltage E via 48, i.e.

with

Note that we have introduce the open circuit potential E⁽⁰⁾ in BV 0: Open circuit potential section as

We can also evaluate the open circuit potential function E⁽⁰⁾ with the interface concentration , i.e.

Note that this is a crucially different to 165 when finite diffusion in the active particle phase is considered, see BV 1: Infinite fast diffusion and conductivity in the active particle and the electrolyte section. However, this allows us rewrite the surface affinity as

with

and

This yields

This is the general, thermodynamic consistent version of the Butler–Volmer equation.^4,26 The specific form 171 of the current density i in terms of the surface overpotential¹ , the open circuit potential , and the exchange current density is widely employed in the literature^28,40 and thus feasible to discuss various material models of .

In Refs. 1,41,42 as well as subsequent work we find

with η = Φ₁ − Φ₂, where Φ₂ "is measured with a lithium reference electrode",¹ p.1527 and Φ₁ the electrostatic potential in the active phase. For the exchange current density i^BV₀ we find various models:

• In Ref. 41 we find
  
• In Ref. 41 we find
  
• In Refs. 35,43 we find
  

This model for the Butler–Volmer-reaction rate became a standard in the literature of modeling intercalation batteries^44–48 and is implemented in various software packaged to simulate battery cycles (i.e. COMSOL, Battery Design Studio,⁴⁹ BEST^35,43) as well as a basis for the interpretation of experimental data.⁵⁰

We compare the Butler–Volmer Equation 172 to the surface reaction rate 171 and discuss the thermodynamic consistency of the three models 173–175 for the exchange current density. Latz et al.,⁵¹ Bazant,²¹ and Dreyer et al.²⁶ also point out the importance of thermodynamic consistency of the Butler–Volmer equation to achieve some overall predictive model since it couples the different thermodynamic bulk models. Dreyer et al.²⁶ use the same structure 38 of the reaction rate, in combination with a constant exchange current density, however do not study intercalation electrodes.

First of all we mention again that the Butler–Volmer Equation 171 is derived from surface thermodynamics (see Reaction rate based on surface thermodynamics section) and that the exchange current density is the Onsager coefficient of the surface reaction Onsager coefficient of the intercalation reaction. This yields some necessary constraints on in terms of the functional dependency on the concentrations (or mole fractions) evaluated at the interface , which are discussed in Onsager coefficient of the intercalation reaction.

By comparison of Eqs. 171 and 172 we obtain

For a metallic lithium counter electrode, where the reaction

is in thermodynamic equilibrium we have 45 entails

which somehow justifies the interpretation of Φ₂ as the potential "is measured with a lithium reference electrode",¹ p.1527. However, from a thermodynamic point of view this re-definition of the potential is not necessary and could lead to inconsistencies when not applied in all balance equations (e.g. of the electrolyte transport) and boundary conditions of the intercalation battery model.

For the exchange current density we showed in Onsager coefficient of the intercalation reaction section that if is dependent on the concentrations at the interface , the dependency for the electrolyte species is necessarily

and for the intercalated ions in the active phase

or overall

with . Comparing these constraints with the models of the exchange current densities 173–175 clearly shows that the dependency of i^BV₀ on the mole fraction (or concentration ) of the electrolyte concentration is not compatible with a reaction rate based on non-equilibrium surface thermodynamics. The concentration dependence is already embedded in the term of the surface affinity 167. A dependency of the exchange current density on is in principle compatible with surface thermodynamics. All three models propose

which in terms of the surface Onsager coefficient would be

In order to discuss the validity, predictability and finally the necessity (or non-necessity) of a concentration dependent surface Onsager coefficient (or exchange current density), we pursue the same strategy and scalings as in Discussion section, however, now with the model 183. We compute the cell voltage E as function of the capacity and the C-rate C_h in the hierarchy of approximations BV 1 – BV 5 and compare it to the computations based on the constant Onsager coefficient.

Eq. 51 reduces with negligible double layer contributions to

We consider again the scaling

which yields the cell voltage

Consider the approximation of infinite conductivity in both phases as well as infinite fast diffusion in the electrolyte, i.e. the approximation BV 3. Fig. 17 shows computations of cell voltage with constant exchange current density as well as concentration dependent exchange current density, for slow () and fast () diffusion in the active particle phase.

Zoom In Zoom Out Reset image size

The impact of the model 183 for the Onsager coefficient (or the exchange current density) on the cell voltage is surprisingly small. Quite similar to the assumed concentration independence of the diffusion coefficients and we can conclude that is a rather good approximation for the overall modeling procedure. However, well defined and reproduce experimental data sets to compute absolute and relative model errors are rare throughout the literature and the deviations in 17 within the experimental variability. We conclude hence that the model 183 is in principle thermodynamically consistent, when embedded rigorously as stated in Modeling section, however, a constant exchange current density produces also very reasonable results and is thus the first choice.

In this work we discuss the cell voltage E of a non-porous intercalation half-cell during galvanostatic discharge with a continuum model for the active intercalation phase, the adjacent electrolyte, and boundary conditions coupling the phases. Based on non-equilibrium surface thermodynamics a reaction rate for the intercalation reaction is stated and the measured cell voltage E subsequently derived. We emphasize some necessary restrictions on the exchange current density of the surface reaction rate in terms of concentration dependence to ensure surface thermodynamic consistency.

For the detailed investigation of the non-equilibrium processes, scalings of all non-equilibrium parameters, i.e. the diffusion coefficients and of the active phase and the electrolyte, conductivity and of both phases, and the exchange current density of the intercalation reaction, with respect to the 1-C current density are introduced. The current density i, entering the model via the boundary conditions, is then expressed as multiple of , i.e. i = C_h · i, where C_h is the C-rate. Further we derive an expression for the capacity Q of the intercalation cell, which allows us to compute numerically the cell voltage E as function of the capacity Q for various C-rates C_h. Within a hierarchy of approximations, e.g. open circuit potential, infinite conductivity, infinite fast diffusion, and so forth, we provide simulations of E = E(C_h, Q) for various values of the (non-dimensional) parameters (), scaled with respect to the material constant . This provides an overall view of the processes and scalings within a lithium ion half cell which is validated at experimental data of Li_xNi_1/3Mn_1/3Co_1/3O₂(NMC).

We consider complete dissociation of the electrolyte and can thus express the mole fractions y_α in terms of , i.e.

and according to the electroneutrality condition. Note we assume

whereby

We can also express y_α as function of of , i.e.

Further

whereby

and thus

Finally we have also

We want to solve numerically the problem

with boundary conditions

and

Note that it is ever convenient for the numerical computation of to introduce the variable

which yields

The parameter a can be adjusted for numerical computations.

Correspondingly, we obtain

and

with

This yields

with boundary conditions

and

The initial value is

and transfers as

PDEPE takes the form

with boundary conditions

We have hence

and

We consider for the chemical potential in the active particle phase

with

Hence

with

The thermodynamic factor is then

We seek to solve 109, i.e.

with boundary conditions 110

and

Note that it is ever convenient for the numerical computation of to introduce the variable

which yields

We have hence

and

with

PDEPE takes the form

with boundary conditions

We have hence

and

Note that we introduce the stop-event for the time-integration of pdepe.