Properly Lumped Lithium-ion Battery Models: A Tanks-in-Series Approach

In the absence of complexities such as phase-separation or concentrated solution effects,^38,39 solid phase transport is modeled by Fick's second law in spherical coordinates. For the positive electrode particle, we have

Where the r coordinate denotes the radial distance within the particle and is thus the second 'pseudo-dimension'. The subscript 1 denotes variables and parameters pertinent to the positive electrode (region 1).

In porous electrode theory, Equation 1 is valid at each point along the electrode thickness x. The superscript s is used to denote the solid phase. D^s₁ is the diffusion coefficient in the positive particle. The second-order problem in r requires the specification of two boundary conditions. At the surface of the solid particle, the diffusive flux is related to the local rate of electrode reaction, or the pore-wall flux as

A symmetry boundary condition is applied at the center of the positive particle

The analogous set of equations are written for the negative electrode (region 3)

With the boundary conditions given by

The Butler-Volmer equation is a common constitutive relation for the pore-wall flux in each electrode as follows

Where k₁ is the rate constant for the positive electrode reaction, and c^{s, max}₁ denotes the maximum concentration of Li in the positive electrode particle. F and R denote Faraday's constant and the gas constant respectively. α's are the charge transfer coefficients for each electrode reaction. The quantity η = ϕ_{s, 1} − ϕ_{l, 1} − U₁(c^s.surf₁) is the surface overpotential, expressed as the difference between the solid and liquid phase potentials minus the electrode open circuit potential U₁(c^{s, surf}₁) vs. a Li/Li⁺ reference electrode. c^{s, surf}₁ is the solid particle concentration evaluated at the surface of the particle, i.e.

The dependence of open circuit potential (OCP) on the surface concentration is indicated accordingly. The equivalent expression for the negative electrode is given by

Equations 1 and 4 can be volume-averaged over their respective electrode volumes. For the positive electrode, we obtain

The volume-averaged form becomes

Similarly, for the negative electrode, we have

The corresponding boundary conditions can also be expressed in volume-averaged form.

Numerical solution of these equations entails spatial discretization in the spherical dimension. Discretization of Equations 11 and 12 results in a system of Differential Algebraic Equations (DAEs), with a convenient linear form for constant D^s_i. For discretization, we employ an efficient three-parameter model based on a biquadratic profile for the radial dependence of c^s_i.^14,40 This approximation is expected to ensure higher accuracy than a two-parameter parabolic profile even at relatively high rates of discharge. The discretized system of equations is therefore

Where the three-parameter model has been expressed in terms of the particle-averaged solid phase concentration , the particle-averaged concentration gradient , and the particle surface concentration . The particle average concentrations are related to the State of Charge (SoC) at the cell level and is directly obtained from the simulation results in the above formulation.

As mentioned earlier, the focus of this article is the development of efficient equations for the electrolyte phase, and therefore the most commonly adopted solid-phase approximation is used in this article. While other reformulation and approximation techniques may be more suitable at higher discharge rates and parameter combinations, we chose this approximation to (a) explain the concepts with a simpler approximation for brevity and easier adoption of current work (b) to alert users the importance of more detailed and relevant approximations published elsewhere.¹⁵ Importantly, the accuracy of the quartic profile approximation was quantitatively verified against nearly error-free numerical methods (collocation, finite difference) for the cases considered in this article.

For charge transport in the solid phase, the governing equation may be written as a conservation law for charge as follows^37,41

The time-derivative for charge density is ignored due to electroneutrality.

Similarly, we have, for the negative electrode

Here, are i_{s, 1} and i_{s, 3} denote the solid phase current densities. The constitutive equation for the solid phase current density is an Ohm's law expression based on the effective electronic conductivity and local potential gradient as follows

Where 1 − ε_i − ε_{f, i} is the fraction of solid phase in electrode i, after subtracting the liquid and inert volume fractions. This factor corrects for the actual conduction pathways in the electrode material. More detailed correction factors may also be applied. The solid phase volume fraction also appears in the specific interfacial area, which for perfectly spherical particles is given by .

Volume-averaging gives us

Using the boundary conditions that impose the solid phase current density at the interfaces, one can simplify Equation 20 as

Volume averaging thus connects the applied cell current density i_app and average pore-wall flux. The sign convention for the model is so adopted that i_app is negative during discharge.

Similarly, for the negative electrode, we have

Equations 21 and 22 can be used in conjunction with the volume - averaged forms of Equations 1–6 to determine the temporal evolution of average solid-phase concentrations for a given i_app. In effect, the active material in each electrode is now modeled by a single representative particle,¹⁸ the Li concentration profiles through which will be influenced by factors such as the applied current density i_app, solid phase diffusion coefficients D^s_i and the characteristic particle radius R_i. This set of equations also determines the evolution of the averaged surface particle concentration, which in turn affects the surface overpotentials and thus the cell voltage response through constitutive Equations 7 and 9. These equations are volume-averaged in order to obtain a relationship between the average pore-wall fluxes and the average potentials and . This is illustrated for the positive electrode in Equation 24 below

Unlike the equations for electrolyte concentration, the highly non-linear nature of the constitutive equation renders evaluation of Equation 24 cumbersome and likely impossible without the use of a full-order solution that gives the actual spatial dependence of the variables. To this end, averages of non-linear quantities are approximated by their value at the average values of the variables on which they depend (i.e. ).

Mathematically, this can be stated as

The classic Single Particle Model (SPM) employs an additional simplification by ignoring the dynamics of the electrolyte phase. Therefore, , implying that the electrolyte concentration is always equal to its initial value. Neglecting liquid phase variations also means that is often set to a constant reference, e.g. for all i.

Neglect of ohmic and electrolyte concentration effects restricts the accuracy of SPM to operating regimes characterized by low ohmic losses, low currents, and kinetically limited electrodes, which usually result in spatially uniform pore-wall flux distributions.^16,30 To this end, the Tanks-in-Series descriptions of electrolyte dynamics are expected to augment and improve the practical applicability of SPM.

We begin with the governing equation for electrolyte concentration for an isothermal model in one spatial dimension. The equations for c based on porous electrode theory may be expressed in the form of conservation laws³⁷

In the positive electrode,

Due to the absence of solid active material, the conservation equation for c in the separator is characterized by a lack of a source term as

The subscript 2 is used to denote the variables in the separator domain.

The equation for the negative electrode is identical in form to that of the positive electrode

N₁, N₂, N₃ may regarded as electrolyte fluxes, which need to be related to local concentration gradients. Noting the similarity of the governing equations to Fick's second law, we have the following constitutive equations⁴¹

In the above equations, D(c) denotes the concentration-dependent electrolyte diffusion coefficient, corrected by a Bruggemann-type factor to account for porous medium tortuosity.

The governing equations for electrolyte concentration are second-order in space, which entails the specification of two boundary conditions for c₁, c₂ and c₃. The boundary conditions are defined at the extremities of each domain. The positive and negative current collectors are physical barriers to the transport of Li⁺ ions, and thus the electrolyte flux at these locations is set to zero. These boundary conditions are thus given by

And,

In addition, electrolyte concentrations and their fluxes must be continuous at the interface between the separator and electrodes. At the positive electrode-separator interface, this is expressed as

In general N_ij is used to denote the flux at the interface of regions i and j.

Similarly, at the interface between the separator and negative electrode, we have

Now, Equation 26 is integrated over the volume of the positive electrode V₁ as

For the one-dimensional model in cartesian coordinates, the differential volume dV is given by dV = Adx, where A is a constant that may be considered a cross-sectional area. In addition, we express the integrals in terms of average quantities as , and , with denoting the volume average of variable v in a given 'tank'. Substituting these relations converts the volume integral into a one-dimensional integral over electrode thickness, i.e. from x = 0 to x = l₁.

Equation 34 thus becomes

Here, we make the reasonable assumption that the electrode porosity ε₁, specific interfacial area a₁ and Li⁺ transport number t⁰₊ in the electrolyte phase are constant in both space and time.⁴²

Using Equation 30, Equation 35 reduces to

The same sequence of operations gives us the volume-averaged equations for the separator and negative electrode as follows

It is worth noting that the steps applied so far represent a rigorous volume-averaging of the equations in each porous domain, followed by the use of the boundary conditions to eliminate interfacial flux terms where possible. No approximations have been made up to this point.

In order to track the average concentrations in each 'tank', we begin with the volume-averaged concentration Equations 36–38. Inspection of these equations reveals the presence of the unknown interfacial flux terms that require suitable approximations to achieve closure. In doing so, we can exploit the flux boundary conditions 32 and 33, which establish the mass flow coupling between adjacent tanks in series.

A simple approximation for the interfacial diffusive flux is in terms of a 'driving force' Δc and a 'length scale' approximation δ_{i, ij} within domain i for the interface between domains i and j. The concept is depicted in Figure 3. This is analogous to the 'diffusion-length' approach attempted previously, but we dispense with assumptions on the spatial profiles for c.^32,43 The 'driving force' Δc is expressed in terms of a difference between the average concentration and the unknown interfacial concentration , and we use as a first approximation. We therefore have, using the constitutive Equations 29

On the separator side of the interface, we assume , which we use for the flux approximation

An equivalent interpretation of the above flux approximations is that we have assumed that the concentration at the midpoint of the porous domain equal to the volume average. This approximation is mathematically equivalent to Gaussian integration with one point, accurate to l²₁.

Zoom In Zoom Out Reset image size

Substitution of the above approximations into the continuity conditions of Equations 32 gives us

The interfacial concentration is now expressed in terms of tank averages as

In general, we use the notation v_ij to denote the value of a given variable v at the interface of domains i and j. An identical sequence of steps yields the concentration at the separator-negative interface as

Substituting the values of interfacial concentrations back into the flux approximations allows us to express the interfacial fluxes in terms of tank-average variables. Thus, we have

Equations 44 and 45 can thus be inserted into the volume-averaged forms 36–38 to obtain a system of Ordinary Differential Equations (ODEs). It is important to equate the approximations for the total flux, and not just the driving forces.⁴⁴ This ensures true mass conservation.

In the three 'tanks', after substituting the known values of average pore-wall fluxes, we have

As specified previously, a slight difference between this approach and others in literature^18,30,36,45 is that it avoids the assumption of uniform reaction rate (given by the pore-wall flux) to solve the PDEs for concentration, but instead deals in average quantities. This allows the prediction of average concentration trends even when the constant pore-wall flux assumption is not applicable, with the flux approximations as the sole source of error. Inspection of Equations 46–48 also suggests their decoupling from those for other electrochemical variables, indicating that they may be solved independently as an ODE system, the solutions of which may be used to compute other relevant quantities during post-processing. While all model equations are simulated simultaneously throughout this article, such a segregated approach may be computationally efficient in real-time control or resource-constrained environments and is enabled by the Tanks-in-Series model. Leaving the model in this form allows for the incorporation of nonlinear diffusivities.

The governing equation for electrolyte current is related to charge conservation in the electrolyte phase. Thus, we have

Where the constitutive equation for electrolyte current is given by a modified Ohm's law based on concentrated solution theory. Thus

The concentration-dependent ionic conductivity κ(c) is corrected by a tortuosity factor specific to each region.

Now, volume-averaging Equation 49 is redundant, ultimately resulting in Equation 21 due to the overall charge balance imposed by porous electrode theory. However, the interfacial boundary conditions for i_{2, i} provide for the estimation of liquid phase ohmic effects. At the interface between the electrodes and separator, the entire current i_app is carried by the liquid phase, and the solid-phase current density is zero. Using the constitutive Equation 50, we therefore have, at the positive electrode-separator interface

Where we have defined the thermodynamic factor as . The tank-averaged equations for the electrolyte potential are now written as

Where an approximation and , analogous to Equations 39 and 40 is made for the gradients of ϕ_{l, i}, in terms of tank-average and interfacial values. Equating the interfacial current density results in Equation 53, which can be solved in conjunction with continuity conditions to obtain the interfacial electrolyte potential.

Solving for ϕ_{l, 12} gives us

A similar form is obtained at the separator-negative electrode interface

Equations 39–43 can now be combined with Equations 51–55 to give us the algebraic equations governing electrolyte potential

The final step in the model formulation is the specification of a reference potential. A convenient reference is the electrolyte potential at the interface between the separator and positive electrode. We thus set ϕ_{l, 12} = 0. This modifies Equation 54, and completes the DAE system for the Tank Model.

We can now assemble the complete Tanks-in-Series Model in Table IV.

Table IV. Governing Equations of the Tanks-in-Series Model.

Positive Electrode (Region 1)	Separator (Region 2)	Negative Electrode (Region 3)

Solving the system of Table IV now allows the prediction of the cell voltage as

In using Equation 58 to calculate cell voltage for the Tanks-in-Series Model, it is implicitly assumed that the solid phase potentials at the cell termini may be approximated by their respective electrode averages. In reality, the terminal potentials must be determined via interfacial approximations for the constitutive Equations 19 analogous to Equations 52 for ϕ_{l, i}. Accounting for this potential drop may be necessitated at high current densities, or for electrodes with poor electronic conductivity. In practice, σ_eff ∼ 50 S/m and the ϕ_{s, i} gradients within the electrode are negligible. This can be seen by the following rough estimation for an aggressive |i_app| = 100 A/m² and σ_eff = 1 S/m, and electrode thickness l₁ = 100μm. We therefore have

Thus, the upper limit of this solid phase ohmic drop is ∼5 mV per electrode. This justifies the assumption of uniform ϕ_{s, i} for most situations of practical salience.

In this section, the Tanks-in-Series model is simulated for galvanostatic discharge at C-rates of 1C, 2C, and 5C. In comparing the model predictions to the full p2D model, we seek to ascertain the accuracy of the Tanks-in-Series approximations, and to quantify its deviation from p2D as a function of applied current density. The results are also compared against curves from SPM, chiefly to illustrate improvement over commonly used fast physics-based models.

Galvanostatic discharge curves

Figure 4 compares the cell voltage-time predictions of our Tanks-in-Series model (hereafter termed the 'Tank Model') with SPM and the full p2D model. The agreement of all three models at 1C and 2C indicates relatively low liquid phase polarizations compared to other overpotentials, possibly due to the values of the Base Case parameters, which are those of a power-cell rated for higher discharge rates. Even so, there is a discernible improvement with the Tank Model. The improved accuracy of the tank model at 5C discharge is substantially more pronounced, particularly during the intermediate phase of the discharge process. This is further illustrated in Figure 5, which depicts the instantaneous error for V_cell. The qualitative trends are identical for both the Tank Model and SPM, but the curves for the Tank Model appear shifted by nearly a constant value compared to SPM. This difference is due to the Tank Model's estimates of liquid phase ohmic drops and concentration overpotentials, which reduces error. The magnitude of these overpotentials increases with current density as characterized by the C-rate.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The increased accuracy of the Tank Model at 5C is illustrated in Figures 5c, and 5d, with the absolute error only slightly exceeding 20 mV, in contrast to >60 mV for SPM. A more than threefold reduction in RMSE is obtained at 5C due to the incorporation of the lumped equations for the electrolyte.

The oscillatory instantaneous error profiles in Figure 5 can be observed in other reduced-order models.^18,25,26 The maximum amplitude of the oscillations, corresponding to the maximum absolute error, expectedly increases with discharge rates. Sources of this error may be in the liquid phase approximations of the Tank Model, and in the estimation of electrode averaged overpotential from the electrode averaged current as in Equation 25. In this case, the contribution of errors from the solid phase concentration is expected to be negligible given the choice of particle size R_i, which ensures rapid diffusion dynamics and attainment of steady state in ∼50 s. This can also be seen in Figure 6, in which the averaged particle surface concentrations across the three models at 1C are compared. The close agreement of the Tank Model to the finely discretized SPM and p2D models indicates the accuracy of the biquadratic profile approximations 14–16. Even at a higher current density of 5C, agreement is ensured in both electrode particles, as depicted in Figure 7. However, the average values may still obscure variations across the electrode thickness due to a non-uniform reaction distribution. In particular, the particle surface concentration at the termini will also deviate from its electrode average due to variations in local pore-wall flux j_i. Such deviation from the average values of Equations 22 and 23 results in errors in reaction overpotentials, which contributes to errors in V_cell.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 8 compares the spatial distribution of pore-wall fluxes at 1C. The positive electrodej₁ achieves rapid uniformity, even though the reaction is initially confined near the separator interface at t = 0. Minor oscillations and deviations from the average are observed toward the end of discharge, but the relative change is negligible compared to that for the negative electrode j₃. The local value near the separator is six times higher than the Tank Model average at the beginning of the discharge, whereas j₃ ∼ 0 near the negative collector. The maxima shift toward the negative collector as the discharge progresses. The reaction 'front' propagates through the electrode, and successive peaks are observed at the two points. Throughout the discharge, the local intercalation rate near the terminal deviates substantially from the Tank Model, sometimes by a factor of two, approaching the average only toward the end of discharge. The reaction distribution in a porous electrode is the result of the relative balance between reaction kinetics, ohmic resistances, and the functional dependence of the open circuit potential U(c^{s, surf}_i). Non-uniform reaction distributions are observed in electrodes with faster kinetics with respect to ionic or electronic transport, and flat OCP curves.^17,49 The U(c^{s, surf}_i) curve for graphite has a substantially flat section corresponding to multiple phase transitions at different degrees of lithiation.⁵⁰ In contrast, U(c^{s, surf}_i) for the NCM positive electrode exhibits a more monotonic dependence with a larger (in magnitude) average slope over the full lithiation range. The non-uniform reaction distribution thus contributes to the error in the Tank Model, in addition to the approximations for liquid phase.¹⁹ These effects are even more pronounced in Figure 9, which examines these trends at 5C. The larger relative deviation close to the negative collector is expected at the higher overall current density. Additionally, oscillatory trends are also observed for the positive electrode beyond t ∼ 300 s, with the deviations about the average increasing to ∼20%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Despite the chemistry-specific considerations and limitations of the Tank Model equations, the Tank Model results in an RMSE of 14.3 mV even at 5C discharge rate, making it competitive in terms of error metrics for online applications. This error is expected to reduce even further based on the specific battery chemistries being considered, such as in the case of a negative electrode with slower kinetics and a more monotonic U(c^{s, surf}_i) curve,¹⁷ which is expected to result in an inherently more uniform reaction profile. This would further bolster the case for the Tank Model with its substantially improved computational efficiency.

Given the approximate nature of the Tank Model equations, it is also expected that errors in estimating the liquid phase overpotentials will contribute to the error. In order to characterize the extent of this error, we now examine the liquid phase variable predictions. The next section evaluates the Tank Model predictions of liquid phase quantities.

Figures 10a and 10b compares the average and interfacial concentrations from the Tank Model and the p2D models, at 1C discharge rate. While there is a disagreement in the steady state values of attained, the deviation from the p2D model is ∼2% in the two electrodes. On the other hand, there is near perfect agreement for the average electrolyte concentration in the separator . The interfacial concentrations also exhibit agreement, despite the somewhat naïve approximations used for interfacial flux. This is likely due to the smaller thickness of the separator relative to the electrodes (25 μm vs. ∼40 μm) which improves the validity of our two-point flux approximations in Equations 44 and 45. The Tank Model predicts different dynamics for compared to the p2D model. An approach to a constant steady state value compared is indicated, in contrast to the p2D model, which shows small fluctuations. The establishment of this steady state occurs on the characteristic electrolyte diffusion timescale, given by in the electrodes.

Zoom In Zoom Out Reset image size

The concentration drop across the separator is ∼ 40 mol/m³ given its smaller thickness relative to the electrodes. Therefore, the diffusion coefficient at both interfaces is comparable in magnitude, since the concentration change is within 4%. Coupled with the comparable thicknesses of both electrodes (Table III), this results in similar magnitudes of interfacial concentration drops in both electrodes due to the imposition of flux continuity by the Tank Model as in Equation 41.

Figures 11a and 11b compares the concentrations at 5C discharge rate. The close agreement between the concentrations is notable, indicating the ability of the Tank Model equations to capture the concentration variations. The average concentrations in the electrodes agree within 5%, with near-exact agreement for the values in the separator, and at the interfaces. The transient dynamics and concentration drop in the electrodes for both models are qualitatively similar to the 1C case, though the magnitude of the concentration drop has increased to sustain the higher flux at 5C. The close agreement between the concentrations profiles suggests that the contribution of concentration overpotential errors to V_cell is not significant up to 5C for the cell parameters considered here. It also indicates the accuracy of the approximation, likely due to the electrode thicknesses considered. It may be possible to refine these flux approximations by altering the diffusion lengths δ_{i, ij}. This can be used to in order to match the concentration profiles with attendant benefits for the accuracy of prediction of V_cell through the concentration overpotentials. For example, it is common in literature to assume a parabolic or cubic spatial dependence for c_i in the electrodes.^{24,25,27,34,45} Applying volume-averaging on a parabolic profile in a given domain returns a value of , the use of which may result in more accurate predictions.^32,34 It is worth noting however that the δ_{i, ij}value can be treated as an adjustable parameter that can be estimated, providing a representative measure of the diffusion length scales for a given operating situation. This allows more flexibility, in that different values may be adjusted for different domains without altering the model formulation, in contrast to assuming spatial profiles in a somewhat heuristic fashion, which entails rederiving the DAE system each time a spatial profile is changed.

Zoom In Zoom Out Reset image size

Figures 10c and 10d depict the average and interfacial values of the electrolyte potential ϕ_l over a 1C discharge. Of note, the average values rapidly reach a steady state for the Tank Model. We observe excellent agreement in the separator . This indicates the validity of the chosen flux approximations in the relatively thin separator, which also ensure agreement between the interfacial values ϕ_{l, ij}. This is a likely consequence of the charge conservation imposed by the Tank Model through Equation 56, which adjusts the values of interfacial potential to maintain the total liquid phase current i_app.

More significant differences are observed in the electrodes. In particular, the value predicted by the p2D model exhibits substantial fluctuations, whereas the Tank Model predicts a steady state within ∼30 s of the discharge process. This qualitative difference is likely due to the errors in the prediction of average potential from the average pore-wall flux, owing to the non-uniform reaction distribution in the negative electrode, as discussed previously. The errors due to approximating the terminal overpotential based on the average in the Tank Model is thus manifested in , with the instantaneous deviation from the p2D model tracking the time-varying local reaction rates. In addition, substantial errors are observed at shorter timescales, where the Tank Model predicts the attainment of steady state. The maximum absolute value of this deviation is ∼3 mV, corresponding to approximately 40%. The choice of flux approximation is also likely to affect the deviation of the average value from that at the terminal, in addition to producing mismatches in ohmic drop.

In contrast, the uniform reaction distribution in the positive electrode results in a more spatially uniform overpotential profile over the discharge process, which is also reflected in the nature of the profiles. The average values reach quasi-steady values differing by ∼1.5 mV. The sum of the average errors in both electrodes value is ∼3 mV, which is comparable in magnitude to the voltage errors in Figures 5a and 5d. The main contributions to the error V_cell are thus the liquid-phase potentials in the electrodes.

These effects are further accentuated at 5C, as illustrated in Figures 11c 11d. In particular, the initial deviation in has increased beyond 10 mV, with a commensurate increase in average errors in both electrodes. The deviations in both electrodes are ∼ 7 mV, indicating comparable errors in ohmic drops. Interestingly, the profile for from the p2D model also shows a slight fluctuation at t ∼ 300 s, corresponding to the non-uniformities in j₁ at 5C as seen in Figure 9a. The total average error in electrolyte potentials is comparable to the values of Figures 5c and 5d. This indicates a contribution of errors in to V_cell through . The separator and interfacial values ϕ_{l, ij} are in close agreement even at 5C, except for relatively small errors in ϕ_{l, 23} during the approach to steady state.

Refinement of the gradient approximation for may help reduce errors by matching the ohmic drops across the electrodes, and by improving predictions of reaction overpotential. As with electrolyte concentration, it is common to assume polynomial profiles for the spatial variation of ϕ_l.^26,32 Similar adjustments of the analogous length-scale approximations may help improve estimates of the migration contribution to ionic current, helping reduce errors in the estimation of both ohmic drops and reaction overpotentials. Preliminary studies suggest that altering gradient approximations marginally increases the errors for and , and thus for the positive overpotential , but this error is more than balanced by the reduction in errors for , and the overall electrolyte ohmic drop.

Interfacial fluxes

Figure 12a depicts the comparison of the Tank Model approximations N_{i, j} with the actual interfacial values from the p2D model. The interfacial fluxes exhibit excellent agreement even with the somewhat simple approximations for the Tank Model. Discrepancies in the temporal profiles are observed at relatively short timescales, where the full model predicts a different trajectory for the approach to steady state and subsequent minor fluctuations (not shown). However, this error is less than 1% even for 5C discharge, as in Figure 13a. The close match between fluxes is observed once the steady state gradients have been established at approximately 30 s, which, expectedly, is of the order of the diffusion timescale for the electrolyte. For the Tank Model, the match between the fluxes at the two interfaces is also consistent with the observation of nearly equal interfacial concentration drops and comparable electrode thicknesses, as discussed in the section on electrolyte concentration. There is a rapid establishment of steady state. We expect the approximations to reduce in validity as the current density is increased further. Errors will also arise if the characteristic timescales for different processes are increased, which will result in significant spatial gradients (analogous to the reaction distribution). This in turn is expected to result in errors of concentration overpotentials, with its manifestation in errors in V_cell.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 12b compares the approximations for the ohmic contribution to the electrolyte current density at 1C, to evaluate the approximations for the electrolyte potential gradients. Here too we observe substantial agreement with respect to the p2D model. This is expected given the match in the interfacial mass flux predictions above. The Tank Model ensures charge conservation at the interface through Equations 52. Agreement between the diffusional contributions forces the equality of the ohmic contributions to maintain the current density at the interfaces. The relative error is comparable even at 5C, as depicted in Figure 13b. There is a slight discrepancy in values during the approach to the quasi-steady value and in subsequent fluctuations, which tracks the errors in electrolyte potentials discussed previously.

Given the close agreements in predicted concentrations and concentration overpotentials, the flux approximations for may be modified to reduce errors in and . One easy way to achieve this could be to use a length scale approximation corresponding to a higher order polynomial profile. For the cell and current densities considered in this article both concentration overpotentials and liquid phase ohmic drops are small relative to the 'thermal voltage' RT/F.^5,18 However, the errors in ohmic potential drop and concentrations could cause errors in V_cell at higher current densities, which will be reflected in systematic errors in both |N_{i, j}|and .

The tanks-in-series model thus needs to be examined in light of different limiting performance scenarios, especially for situations in which significant spatial gradients can arise across the electrodes, giving rise to polarizations that are underestimated – e.g. extremely high rates of discharge, or very thick electrodes. This is related to the key assumptions in the Tank Model, namely the interfacial fluxes, which determines the accuracy of the corrections for liquid phase effects.

The impacts of these approximations are now studied in a special case, by the comparison of the Tank Model for the case of a cell with relatively thick electrodes, which may be of salience for high-energy density batteries.^51,52 For this case study, the electrode thicknesses l₁ and l₃ were increased by a factor of 6 while holding all other parameters constant. Thus, the individual capacities of each electrode increase by a factor of 6, while their ratio is held constant. For a consistent comparison with the Base Case, the model was compared at a 5C/6 C-rate, corresponding to the current density i_app equivalent to the 5C discharge for the Base Case. The modified parameter values are listed in Table V.

For this '6x case', significant electrolyte diffusion limitations and ohmic drops are expected to arise compared to the Base Case. A modified parameter S_e based on Doyle et al. can be used as a measure of expected diffusion limitation.^5,53

Thus

This factor is thus the ratio of the characteristic time for electrolyte diffusion to the total discharge time, higher values suggesting exacerbated diffusion resistances.

The voltage-time discharge curve for this case is depicted in Figure 14. Here we can observe the large liquid phase polarization and spatial non-uniformities that SPM is clearly unable to capture, and instead predicts near-complete cell utilization. The p2D model predicts termination of discharge at an intermediate time point, likely owing to severe electrolyte depletion in the positive electrode. What is noteworthy however is the prediction of the Tank Model, which predicts a premature termination but within 5 minutes of the discharge process. This is despite substantially improved agreement compared to SPM over the small portion of its operation. This suggests substantial error in the prediction of c_i and ϕ_{l, i}. This can now be seen in Figure 15, which compares the average concentration predictions from the two models. In particular, the Tank Model predicts the being driven to zero at t ∼5 minutes, which coincides with the termination of discharge. The values disagree by more than 20%, while the Tank Model predicts at t ∼ 300 s, in contrast to the actual value of ∼500 mol/m³. This discrepancy thus suggests the limitation of the diffusive flux approximations, which substantially overestimate the concentration gradients. One possible reason for this error is that the diffusion layers have not built up to their steady state value, as also indicated by the diffusivity parameter S_e. The errors due to the resulting flux approximations may be circumvented by refining the same, or by the use of a time-dependent exponential correction for the length approximation as suggested by some workers.^32,54 The time constants for this exponential change may be defined based on characteristic diffusion timescales.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This example thus demonstrates a combination of design and operating parameters in which the Tank Model has reduced performance and necessitates improved approximations. In particular, the large polarizations considered in this example may be encountered for the Base Case cell as well, such as during high-current dynamic and pulsed operation. A more comprehensive analysis is useful to establish the parameter combinations at the cell level that deem the Tank Model useful. The limits of a BMS that uses such a model can also be defined through such an analysis.

In this article, the Tank Model is compared against SPM because of its substantial ubiquity in advanced BMS applications, as it continues to be the model of choice when physical detail is desired. In introducing the Tank Model, we seek to propose an alternative for various applications where SPM is substantially common. However, for additional perspective, this section discusses the error metrics of the Tank Model relative to a recent version of the SPMe, which corrects SPM with electrolyte dynamics as mentioned previously.¹⁸ For this study, the electrode parameters from Ref. 18 were chosen, in addition to electrolyte transport correlations from Ref. 42. The representative values of electrolyte transport properties required by the SPMe formulation were evaluated at the initial electrolyte concentration.

The error metrics up to 3C discharge rate are illustrated in Figure 16. The substantial improvement in error in going from SPM to lumped electrolyte models is evident, particularly at 3C, where the RMSE reduction is nearly three-fold for SPMe. However, the Tank Model exhibits even lower error (13 mV vs. 35 mV). This suggests more accurate estimation of concentration overpotentials and liquid phase ohmic drops for the chosen parameters. The SPMe uses representative values, whereas the Tank Model only requires the evaluation of electrolyte transport properties at the interfaces. Consequently, the concentration-dependent transport properties can be included in an efficient manner. The expressions in SPMe may need to be rederived to account for these variations.

Zoom In Zoom Out Reset image size

Computational times for the Tank Model are listed in Table VI. For the Tank Model, each conservation law in the electrolyte phase is replaced by its volume-averaged form, while the solid phase in each electrode is replaced by 3 linear DAEs. The original ∼ 200 – 1000 DAEs are thus replaced by 14 average conservation equations. This results in a 1C discharge curve being simulated in ∼ 2 ms, in contrast to >1000 ms for a standard Finite Difference implementation. The computational speed of this model is comparable to fast SPM implementations. The Tank Model is also competitive with the state-of-the-art reformulated p2D model from our group, achieving up to an order of magnitude reduction in computation time. The mathematical similarity of the Tank Model to the reformulated model with n = (1,1,1) collocation prompts the question as to why the (1,1,1) reformulated model is not used instead of developing the Tank Model. This is because the proposed model is conservative and exhibited higher accuracy. In addition, it can be rewritten so as to contain fewer adjustable parameters compared to reformulated models. Increasing the accuracy and numerical convergence of reformulated models requires using its ability to guarantee convergence for different chemistries, parameters and operating conditions by increasing the number of collocation points.

Recalling the volume-averaged electrolyte Equations 46–48 of the Tank Model, we have

The above equations may be conveniently written as

Where the governing equations have been expressed in terms of equivalent 'mass transfer' coefficients, volume and source terms. In the above equations, the have units of mass transfer coefficient (m/s), and V_{l, i} denotes the electrolyte volume (per unit area) in region i with units of length. The transport parameters in the electrolyte current balance 56 may also be grouped similarly

Where, the thermodynamic factor v(c_ij) = 1 for simplicity. Here, have units of S/m². The lengths used in the flux approximation, which for this case is , are present in the expressions for K_ij and S_ij. Estimating these quantities is equivalent to deriving a suitable flux or gradient approximation. For constant diffusivity, the electrolyte concentration equations constitute a coupled linear system of Ordinary Differential Equations with constant coefficients. The form of these equations suggests a possible parameter estimation strategy in which the Tank Model is matched against experimental data by direct estimation of the 'consolidated' transfer coefficients. Variations due to concentration-dependent transport properties, or transients in which the 'steady state layer' is attained over a characteristic timescale comparable to discharge time can potentially be handled by estimating average values over the course of the process. The similarity of the mathematical model to conventional Continuous Stirred Tank Reactor (CSTR) systems is noteworthy and suggests scope for exploiting the extensive control theory developed for CSTR systems.