Efficient Simulation and Reformulation of Lithium-Ion Battery Models for Enabling Electric Transportation

Equivalent circuit models try to describe the underlying system using a representation that usually employs a combination of capacitors, resistors, voltage sources, and lookup tables.³⁷ Capacity fade is often represented by a capacitor with a decreasing capacity, while temperature dependence is modeled by a resistor-capacitor combination. Current research in this area includes adopting the circuit based models by continuously updating the parameters using the current and voltage data.³⁸ Such models occupy the lower left corner of Figure 3; they can be simulated very quickly but are not accurate outside of the operating conditions for which they were developed or as the battery grows older. The parameters also lack any physical meaning, limiting the insight that can be gained from such models.

Despite these limitations, circuit based models are incredibly popular in the BMS literature because of the very small computational requirements of simulation. In fact, it is considered among many that the use of full order models in a BMS is not feasible.^17,39,40 In order to improve the validity of circuit based models, additional components can be included to account for additional phenomena, such as diffusion resistance,¹⁷ hysteresis,²⁵ temperature effects,^26,29 and self-discharge and current inefficiencies.²⁹ This can be achieved by adding linear or nonlinear terms,¹⁷ or using empirical look up tables which are functions of state of charge.²⁵ Plett²⁶ found that calculating the circuit based parameters at discrete temperatures using experimental data did not extend well to temperatures that were not used to determine the parameters. In other words, the parameters did not correlate well with temperature so that linear interpolation did not provide accurate results. This was partially rectified by assuming the parameters were a fourth-order function of temperature.²⁶

Other approaches have used reduced order models based on physical models to incorporate capacity fade effects, such as lithium plating³⁹ and solid electrolyte interface growth.⁴⁰ These are computationally cheap, but ignore the variation of the concentration and pore wall flux across the electrode.^39,40 However, such approximations likely reduce the validity of such models at high rates of charge/discharge.

Moving up the diagonal of Figure 3, the electrochemical engineering community has long employed continuum models that incorporate chemical/electrochemical kinetics and transport phenomena to generate predictions that are more accurate and meaningful than empirical models,^3,23 which can be used for parameter estimation,^41,42 optimization,^12,43 state estimation, and control³⁶ with more confidence than circuit based models. Including additional physical phenomena in a model increases the computational cost in terms of both solution time and memory. Numerical methods often are required as most battery models cannot be solved analytically. The mathematical method used to solve the system of equations can also have a significant impact on the computational cost of simulation.

The single-particle model (SPM) is a simple model that represents each electrode as a single particle²⁰ and considers diffusion in the solid phase, but neglects solution phase effects.^44–46 The governing equations for the SPM are given by describing the diffusion of lithium in an active particle:

with boundary conditions

where the pore wall flux is given using Butler-Volmer kinetics

Battery models are typically solved efficiently using the method of lines in which discretization of the spatial derivatives results in a system of first-order differential algebraic equations (DAEs)^33,34 that can be solved using optimized solvers for initial value problems. DAEs can be difficult to solve because the initial conditions must be consistent with the algebraic equations, which causes many solvers to fail if inconsistent conditions are provided, especially when nonlinear algebraic equations are considered. Techniques for initialization have been presented elsewhere and will not be discussed here.^47,48 Using the SPM with N = 15 node points in each electrode results in a total of 34 DAEs.

This model can be quickly simulated, and has been used to predict capacity fade due to the growth of the SEI layer,²⁰ which makes the SPM a good choice as an initial attempt for implementation in a microcontroller environment. The single particle model has been validated for rates up to a 1C rate of discharge, but the assumptions are not valid at higher rates or thick electrodes where variations in the electrolyte phase are important.^44,45 The SPM can be further reduced if a parabolic profile approximation in the solid phase.⁴⁹ This only tracks the average and surface lithium concentration in the solid phase, reducing the system to only 4 DAEs that must be solved. This makes the SPM very efficient for use in a BMS for low power applications. However, for applications in which higher rates are experienced, a more comprehensive model is needed to accurately estimate the internal states to develop aggressive control strategies.

The pseudo-two-dimensional (P2D) model is a more detailed physics-based model that considers the electrochemical potentials within the solid phase and electrolyte along with lithium concentration in both the solid- and liquid-phases,²³ and is flexible enough to include additional physical phenomena as understanding improves.^{21,45,50–59} The improved predictive capability of the P2D model has contributed to its popularity among battery researchers but it has two independent spatial variables: x to track the variables across the thickness of the cell sandwich, and r to track the lithium concentration radially in the solid electrode particles.²³ Having multiple spatial variables increases the dimensionality of the problem, which greatly increases the number of equations to be solved (and computational requirements) if a finite difference approach is used to discretize both the x and r directions. If 15 node points are used in the radial direction, 50 node points across each electrode, and 25 node points for the separator, nearly 2000 DAEs must be solved. Therefore, appropriate mathematical techniques are required to reduce computational time and memory requirements in order to allow the model to be implemented in a microcontroller environment.^31,42 This high computational cost of simulation has motivated researchers to develop techniques to simplify the battery models and enable faster simulation. For example, proper orthogonal decomposition has been used to reduce the total number of states simulated.⁶⁰ Quasi-linearization combined with a Padé approximation has also been used to simplify the model and improve simulation.⁶¹

Conversely, many commercial software packages, such as COMSOL,⁶² Fluent,⁶³ etc. use well understood numerical methods to solve ordinary differential equations (ODEs) or partial differential equations (PDEs). However, many node points, control volumes, or elements are required for convergence. These methods are robust approaches for solving the problem, but the resulting set of algebraic or differential-algebraic equations can number into the thousands and is computationally expensive, even for linear problems, and is difficult to implement into a microcontroller or other resource-limited environment. Furthermore, many commercial solvers are over-designed in order to handle a wide variety of problems with minimal input from the user. They do not exploit the structure and unique characteristics of the underlying models, which can be used to improve the computational performance without compromising on the robustness.

In order to reduce the number of DAEs that must be solved for the P2D model, the reformulation methods described previously for the SPM can be implemented for the solid phase diffusion in the P2D model. Using the parabolic profile approximation for the concentration profile in the radial direction, the number of DAEs can be significantly reduced, thereby improving computational efficiency.³⁴ For the case with 50 node points across each electrode and 25 node points for the separator, roughly 500 DAEs must be solved, much less than the 2000 for a full finite difference approach. The parabolic profile approximation is valid for long times and low rates, but has inaccuracies when there is a large gradient in the solid phase particles, which become significant for rates greater than about 4C.³⁴ Ramadesigan et al.³² developed a mixed finite difference (MFD) approach for the solid phase using unequal node spacing across the radius of the particle so that fewer node points are required to achieve convergence at high rates.³² The equivalent MFD solution results in approximately 1000 DAEs. Using higher order approximations for the concentration profiles in the radial direction can also be used for greater accuracy at high rates while minimizing the computational requirements.

Reformulation in the x-direction can also be applied to further reduce the computational demands of simulation. Spectral methods have faster convergence than finite differences so that fewer equations are required, but require more up-front work for implementation and the resulting system of equations is not sparse, unlike for the finite difference method. In spectral methods, the unknowns are approximated as a series solution of trial functions, such as cosines, with time-dependent coefficients. The coefficients are determined by minimizing the residual of the governing equations across the domain typically by using the Galerkin or orthogonal collocation (OC) methods,^31,64 though OC can better handle non-linear parameters as the integrations required for the Galerkin approach are computationally prohibitive.

Since each dependent variable is approximated as a series solution, it may be possible to solve some equations analytically a priori. Symbolic math tools such as Maple⁶⁵ or Mathematica⁶⁶ can play an important role in solving for unknown variable to reduce the number of equations that the solver must compute. However, this can increase the complexity of the remaining equations, so testing is often required to determine if this approach is indeed advantageous. Using a single non-constant term in the series approximation coupled with the parabolic profile approximation in the solid phase can provide reasonable results (see Figure 4 and Table II) with only 21 DAEs.

Zoom In Zoom Out Reset image size

Table II shows the simulation time and numerical error for several simulation schemes using the FORTRAN solver DASKR⁶⁷ run on a 3.33 GHz, 24 GB RAM machine for a 1C discharge. The use of the SPM provides faster simulation, but the average error is large compared to any of the P2D simulations due to the limitations of the model itself. Note that the same physical parameters and values used in the P2D model were used in the SPM. It is possible that including a correction factor to account for the electrolyte resistance would provide a better accurate fit. Graphically, Figure 4 shows the voltage-time curve as predicted using the lowest order models considered for implementation in a microcontroller. It is readily apparent that the reformulation approach provides the most accurate results. Conversely, the low order finite difference solutions deviate significantly from the full order solution. The SPM solution predicts a significantly higher voltage for the entire discharge due to its neglect of electrolyte phase resistance, but all the models considered do predict the total capacity reasonably accurately.

Importantly, model reformulation reduces the computational time of the P2D model to be comparable to SPM and circuit-based models which makes the implantation of such models into a BMS practical. Furthermore, using reformation allows for the possibility of using more detailed models (moving further up the diagonal in Figure 3) in a BMS. For example, thermal effects can be included, and/or a P3D model could be used to account for the spatial variation parallel to the electrodes. Additionally, for state of health estimates and life modeling, capacity fade mechanisms must be included in the model at increased computational cost, making reformulation even more useful.

Table I shows the relationship between the choice of models and possible functionalities of the BMS to provide more clarity about the advantages of using a detailed physics based models in BMS. Empirical and circuit based models can optimize conditions which describe the cell as a whole, for example, the cell voltage, or total SOC. Using the single particle model allows for constraints to be applied to electrode averaged variables, such as anode or cathode SOC. The porous electrode model allows for objectives to be set on local values of SOC or potential, as well as other variables. Applying constraints to local variables can be very important, especially at high rates which can cause significant variations across the electrode. Under such conditions, the average values may suggest that there is nothing to be concerned about, but there may be areas within the electrode that experience conditions which are detrimental to performance and/or life.

Inclusion of the stress and other effects into the single particle framework allows for constraints to be implemented to reduce capacity fade in the cell,⁶⁸ while inclusion of the same phenomena into the P2D framework allows from local variation to be accounted for, so that the maximum stress development can be minimized. Under conditions with high spatial variation, the maximum stress (as predicted by the P2D model) may be much larger than the average stress (as predicted by the single particle model). As fracture and capacity fade occur at any point which exceeds the yield stress, the maximal stress is a more important metric to predict internal damage than average stress. In general, more detailed models allow for aggressive, safe operation of batteries by utilizing physics based local constraints. However, this comes at a computational cost and efficient simulation helps bring more physics models to real-time simulation and control for the BMS.

The use of detailed physics based models in a BMS enables the ability to provide additional charge and discharge constraints on the system to maximize performance. For example, circuit based models can enforce a terminal voltage conditions on the battery on charge and discharge. The specific cycling window varies with the specific chemistry used, but is typically constrained to 3 < V < 4. Using the single particle model, however, allows for voltage constraints to be set for each individual electrode. This can be critical at the anode as lithium plating occurs when the overpotential at the anode when η_anode < 0 V vs. Li/Li⁺ which can result in severe capacity fade and possible dendrite formation. Conservative charging protocols can alleviate this problem if using circuit based models, but using a model which can directly estimate the overpotentials can provide the confidence needed to be more aggressive. The importance of considering lithium plating has led to the development of reduced order models which account for lithium plating from a control perspective.³⁹ Although this can be seen as improvement over circuit based models, such a reduced order approach makes several assumptions, such as neglecting variation across the electrode,³⁹ which makes it invalid at higher rates of charge.

As an example, the reformulated P2D model was used to develop an optimal charging profile to maximize the total amount of charge stored subject to a maximum charging rate and a limited charging time. To ensure that lithium plating does not occur, the additional constraint of η_anode > 0 V vs. Li/Li⁺ at all points in the electrode throughout the charging time is enforced. Furthermore, constraints are applied on the solid phase concentration in the anode and cathode to maintain reasonable limits. By applying constraints directly to the overpotential, it is not necessary to apply a constraint on the overall cell voltage, allowing for higher voltages to be used. The optimal charging profiles are given in Figure 5 with maximum charging rates of 2C, 2.5C, and 3C and compared to traditional constant-current constant-voltage charging (CC-CV). Clearly the optimal charging profile allows for greater charge storage by allowing for the maximum current to be applied for a longer time. Although the optimal charge curves are qualitatively similar to the constant-current constant-voltage charge curves, the voltage-time curves given in Figure 6 show that a constant potential condition is never applied, but are allowed to continually increase. Rather, a constant overpotential is applied, as shown in Figure 7, and the anodic overpotential is never < 0 V, thus safe operation is maintained. The overall cutoff voltage of 4.2 V is overly conservative. Furthermore, overpotential increases during the constant voltage portion of the CC-CV charge, suggesting that more aggressive charging could be applied. This allows for a higher charging voltage to be used. Importantly, such an aggressive charging profile begins to deplete the lithium available in the cathode. Thus, the inflection point in the current at the end of discharge is due to the bounds enforced on the solid-phase concentration in the cathode. This results in a sharp drop in applied current to ensure that the cathode is not overcharged (i.e. by removing too much lithium from the metal oxide lattice, such that the process is not reversible). This allows the cell to relax to an extent, as shown by the increasing overpotential and decreasing cell voltage, while still applying some small amount of charge.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Being able to accurately estimate the overpotential in a lithium-ion cell is critical for estimating the amount of power that can be delivered or received by the cell. We showed that applying a cutoff to the cell voltage limits the amount of charge stored. Thus, during a regenerative breaking event, for example, constraining the overpotential rather than cell voltage can allow more energy to be supplied to the battery, especially in conditions in which the cell is nearly fully charged already. Allowing more power to be sent to the battery directly increases the overall efficiency and utilization and reduces the braking energy which must be dissipated as heat. Importantly, the overpotential cannot be directly measured, so an accurate model with effective state-estimation techniques is required if a BMS is to perform such aggressive charging.

Furthermore, Figures 8, 9, and 10 show the current-time, voltage-time, and overpotential curves for a 3.5C rate of charge for optimal charging and for CC-CV charging. In this case, the CC-CV charge results in a negative overpotential at the anode. This can result in lithium-plating as an undesirable side reaction, and ultimately reduce the lithium available for cycling. By using a model which can account for the overpotential, the occurrences of such events can be minimized. Although the optimal charging profile does not apply the maximum charge for as long as the CC-CV approach, the rate of decrease is less, so that more charge is being supplied by the end of charge in the optimized case.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is possible that such a constant-current constant-overpotential charge may cause increased capacity fade, and thus may not be ideal for all conditions, but can improve the amount of charge that can be stored in a short time.

The reformulated P2D model allows for variation to be considered across the electrode in real-time, without the need to make assumptions reducing the physics significantly. This model can be further extended to include additional physical phenomena, including SEI layer growth. For example, the electrolyte concentration difference across the battery can be constrained to be less than a specified value. This can improve energy efficiency by reducing the diffusive resistance that must be overcome on during charging and discharging. Of course, maintaining the concentration gradient to be within a certain amount can reduce the total amount of power supplied to the battery, which highlights the tradeoff between total power and efficiency and safety. The energy efficiency of battery systems can be much lower at high rates due to the large internal resistances which cause a large amount of heat to be generated in the system. Thus, it may be advantageous to charge a cell at a lower rate, if the downsides of high-rate charging outweigh the benefits. A BMS which accounts for all these phenomena can be used to improve efficiency and performance, and ensure safety and longevity.