Modeling and Simulation of Lithium-Ion Batteries from a Systems Engineering Perspective

Empirical models employ past experimental data to predict the future behavior of lithium-ion batteries without consideration of physicochemical principles. Polynomial, exponential, power law, logarithmic, and trigonometric functions are commonly used as empirical models. The computational simplicity of empirical models enables very fast computations, but since these models are based on fitting experimental data for a specific set of operating conditions, predictions can be very poor for other battery operating conditions. Such battery models are also useless for the design of new battery chemistries or materials.

The electrochemical engineering field has long employed continuum models that incorporate chemical/ electrochemical kinetics and transport phenomena to produce more accurate predictions than empirical models. Electrochemical engineering models of lithium-ion batteries have appeared in the literature for more than twenty years.¹⁴ Below is a summary of electrochemical engineering models, presented in order of increasing complexity.

Single-particle model

The single-particle model (SPM) incorporates the effects of transport phenomena in a simple manner. Zhang et al.¹⁵ developed a model of diffusion and intercalation within a single electrode particle, which was expanded to a sandwich by considering the anode and cathode each as a single particle with the same surface area as the electrode.¹⁶ In this model, diffusion and intercalation are considered within the particle, but the concentration and potential effects in the solution phase between the particles are neglected.^16,17 The following typical reactions are considered in each of the particle in the SPM (MO refers to metal oxide):

Due to these simplifications, this model can be quickly simulated, but is only valid for limited conditions, such as low rates and thin electrodes.¹⁷ Greater efficiency can be obtained by including a parabolic profile approximation for the lithium concentration within the particle.^16,18

Ohmic porous-electrode models

Pseudo-two-dimensional models

The pseudo-two-dimensional (P2D) model expands on the ohmic porous-electrode model by including diffusion in the electrolyte and solid phases, as well as Butler-Volmer kinetics (see Fig. 4). Doyle et al.¹⁴ developed a P2D model based on concentrated solution theory to describe the internal behavior of a lithium-ion sandwich consisting of positive and negative porous electrodes, a separator, and a current collector. This model was generic enough to incorporate further advancements in battery systems understanding, leading to the development of a number of similar models.^14,22–32 This physics-based model is by far the most used by battery researchers, and solves for the electrolyte concentration, electrolyte potential, solid-state potential, and solid-state concentration within the porous electrodes and the electrolyte concentration and electrolyte potential within the separator. This model based on the principles of transport phenomena, electrochemistry, and thermodynamics is represented by coupled nonlinear partial differential equations (PDEs) in x, r, and t that can take seconds to minutes to simulate. The inclusion of many internal variables allow for improved predictive capability, although at a greater computational cost than the aforementioned models.

Zoom In Zoom Out Reset image size

Multiscale, multidimensional, and multiphysics electrochemical-thermal coupled models are necessary to accurately describe all of the important phenomena that occur during the operation of lithium-ion batteries for high power/energy applications such as in electric/hybrid vehicles.

Thermal models

Stress-strain and particle size/shape distributions

Intercalation of lithium causes an expansion of the active material, such as graphite or manganese oxide, while lithium extraction leads to contraction. The diffusion of lithium in graphite is not well understood, but some work has been done to model the diffusion and intercalation of lithium into the electrode material.^27,56,57 Since lithium diffuses within the particle, the expansion and contraction of the material will not happen uniformly across the particle (i.e., the outer regions of the particle will expand first due to lithium intercalating there first). This spatial nonuniformity causes stress to be induced in the particle and may lead to fracture and loss of active material,^58,59 which is one of the mechanisms for capacity fade. Various models have been developed to examine the volume change and stress induced by lithium-ion intercalation for single particles.^60–62 A two-dimensional microstructure model was developed⁶³ to extend the stress-strain analysis from single particles and was eventually incorporated into the full P2D model.⁶⁴ These models show that high rates of charging result in increased stress and increased chance of fracture, which can be somewhat mitigated by using smaller particles, or ellipsoidal particles. Additionally during battery cycling, some particles are lost or agglomerate to form larger sized particles, which results in performance degradation. In addition, porous materials rarely have uniform particle size and shape. Some continuum models have accounted for the distribution of particle sizes and its effect on the battery performance,^65,66 for example, through the equation⁶⁵

where is the fraction of total current flowing in solution, N (r) is the number of particles per unit volume of composite electrode with a radius between size r and r + dr in the porous electrode, Y(r) is a function that relates the outward normal current density per unit surface area of a particle to the potential difference, and is the potential difference between the solid particle and the adjacent solution. A promising future direction would be to extend such models to include variations in particle size and shape distribution by (1) writing N in terms of the multiple independent particle coordinates that define the particle shape (typically 3), and (2) replacing the single integral with a more complicated volume integral. The time-dependent change in the particle size distribution due to breakage and agglomeration can be modeled by a spatially-varying multi-coordinate population balance equation:

where f (l, x, t) is the particle size and shape distribution function, x is the spatial coordinate, l_i is the i^th independent size coordinate, l is the vector with elements l_i (typically of dimension three), G_i (l,t) = dl_i/dt is the growth rate along the i^th independent size coordinate (which is negative for shrinkage), h(l, x, t, f) is the generation/disappearance rate of particle formation (e.g., due to breakage and agglomeration), and t is time.^67–70 The expression for h(l, x, t, f) for breakage and agglomeration contains integrals over the f (l, x, t), and the h and G_i have dependencies on additional states such as local lithium-ion concentrations. This model to capture the effects of morphology within a material, called a mesoscale model,^71,72 would enable the material degradation due to spatially-varying and time-varying changes in the particle size and shape distribution to be explicitly addressed.

Stack models

The Kinetic Monte Carlo (KMC) method is a stochastic approach that has been used to model the discharge behavior of lithium ions during intercalation. Such models^76–79 have been used to simulate diffusion of lithium from site to site within an active particle to aid in understanding on how different crystal structures affect lithium mobility⁷⁵ and how the activation barrier varies with lithium-ion concentration.^78,79 Additionally, Monte Carlo methods have been used to predict thermodynamic properties.⁸⁰ KMC has also been applied to simulate the growth of the passive SEI-layer across the surface of the electrode particle, to simulate one of the mechanisms for capacity fade.⁸¹

Molecular dynamics

Density functional theory

Multiple numerical methods are available for the simulation of any particular battery model. For empirical models, analytical solutions are usually possible and can be easily solved in Microsoft Excel or Matlab.⁸⁸ Analytical solutions can be implemented in a symbolic language such as Mathematica,⁸⁹ or Maple,⁹⁰ or Mathcad,⁹¹ or in a compiled language such as FORTRAN or C++. Analytical solutions based on linear model equations often involve eigenvalues, which might have to be determined numerically. For nonlinear model equations, sometimes analytical series solutions using perturbation methods⁹² or other symbolic techniques⁹³ can be derived. Numerical simulation methods are more flexible, with multiple methods available for any particular battery model. The best numerical methods tend to be more sophisticated when moving toward the upper right of the battery models shown in Fig. 3.

For SPMs for a single electrode, analytical solutions have been derived for constant-current operation and cannot be obtained directly for the constant-potential operation, due to the fact that the boundary flux is implicitly determined by the nonlinear Butler-Volmer equation particularly when the open circuit voltage changes with state of charge. At this scale, especially for AC impedance data, analytical solutions are easily obtained and have been heavily used even for estimating unknown diffusion coefficients. A numeric symbolic solution was also derived for the AC impedance response that showed similar results to the analytical solution.^94–96

When two electrodes are included in an SPM, an analytical solution is available for constant-current operation but not for constant-potential operation, for reasons as stated above, or when film formation for the SEI layer is modeled. Beyond SPM and porous electrode ohmic resistance models, analytical solutions are not possible for simulating charge-discharge curves. A SPM with two electrodes consists of a single partial differential equation for each electrode. Conversely, a finite-difference scheme discretized with 50 node points in the radial direction generates 50 × 2 + 50 × 2 = 200 differential algebraic equations (DAEs). Recall that the SPM is computationally efficient but is not accurate, especially for high rates. For P2D models¹⁴ typically the finite-difference approach has been used. A P2D model with polynomial approximation¹⁸ for the solid phase, when discretized with 50 node points in the spatial direction for each variable, results in a system of 250 DAEs for each electrode and 100 DAEs for the separator. Thus, the total number of DAEs to be solved for the P2D model across the entire cell is 250 + 250 + 100 = 600 DAEs. The addition of temperature effects to this model results in 750 DAEs to be solved simultaneously. Stack models are much more computationally expensive, as the number of DAEs is equal to the number of cells in the stack (N) times the number of equations coming from each sandwich. Using the finite-difference discretization of spatial variables in x, y, and r with 50 node points along each direction in a pseudo-3D thermal-electrochemical coupled model would generate 15,000 + 7500 + 15,000 = 37,500 DAEs to be solved simultaneously for a single sandwich.

The speed and accuracy of a numerical method depends upon the complexity of the model equations, including operating and boundary conditions, and the numerical algorithm. The most common numerical methods for simulation of lithium-ion batteries are the finite-difference method (FDM), finite-volume method (FVM, or sometimes called the control volume formulation), and finite-element method (FEM). The main continuum simulation methods reported in the literature for the simulation of batteries can be classified as

• (1)  
  DUALFOIL.²⁶ This software employs Newman's BAND subroutine,⁹⁷ which is a finite-difference method used to simulate electrochemical systems for more than four decades. Symbolic software such as Mathematica⁸⁹ and Maple⁹⁰ can be used for determining analytical expressions for the Jacobians and for generating the associated FORTRAN code for use with the BANDJ subroutine.²³
• (2)  
  FVM with various time-discretization schemes,⁹⁸ which has been applied to P2D models.
• (3)  
  COMSOL⁹⁹/BATTERY DESIGN STUDIO,¹⁰⁰ which implements the FEM/FDM in a user-friendly interface and includes a module that implements the P2D battery model.
• (4)  
  Finite-difference or reformulation schemes in spatial coordinates with adaptive solvers such as DASSL in time.²³

Each approach has its advantages and disadvantages. DUALFOIL is a freely available FORTRAN code. The FDM has been used extensively in battery simulation²³ as it is easy to implement and modify. The FVM is closely related to the FDM but more easily handles irregular geometries. The FEM handles both irregular geometries and heterogeneous compositions, but is much harder to implement by hand, and so is usually only applied to batteries using commercial FEM software such as COMSOL. An advantage of commercial software like COMSOL is ease of use and that the numerical implementation is invisible to the user and results from COMSOL can be directly integrated to MATLAB environment, which is a widely used tool for control and optimization. However, a disadvantage is that COMSOL's numerical implementations cannot be modified by the user to (1) increase computational efficiency by exploiting additional mathematical structure in the model equations or (2) integrate such efficient simulation results into advanced systems engineering algorithms for optimal design, operation, or control in a computationally efficient manner.

When optimization fails while using COMSOL-like codes, detective work is required to determine whether the numerical simulation was robust enough to provide accurate numerical Jacobians. Also, as of today, global optimization methods are readily available only for algebraic equations. Algebraic optimizations can be formulated by discretization of all the variables and parameters including the control variables,^101,102 but these optimization schemes typically have too high complexity to be solvable using existing global optimization software. Many groups are working on the development of optimization software that is more computationally efficient at computing local optima for dynamic optimizations or on ensuring convergence to a global optimum.^103,104

BATTERY DESIGN STUDIO¹⁰⁰ has a module for the simulation of P2D lithium-ion battery models. Adaptive solvers provide advantages in speed compared to fixed time-discretization schemes. Researchers have used DASSL for solving battery models.²³ DASSL/DASPK use backward differentiation schemes in time, which are numerically stable and efficient. For the same set of equations, these adaptive schemes can provide an order of magnitude savings in time. Battery models more advanced than the P2D model are usually solved offline in the literature (an exception is the P2D thermal model from Gu et al.^44,48 and the stress-strain model from Renganathan et al.⁶³).

To understand the importance of capacity fade in a lithium-ion secondary battery system, significant efforts have been devoted to the development of mathematical models that describe the discharge behavior and formation of the active and passive SEI layers. The majority of these models are empirical or semi-empirical.^105,106 Other works have attempted to simulate capacity fade by considering the lithium deposition as a side reaction and the resulting increased resistance.^31,107–111 Others have simulated capacity fade by modeling the active material loss, or change of internal parameter with cycling.^31,108–113 Other researchers have used KMC methods to examine the SEI layer formation at the microscale level.⁸¹ Such a model, however, is computationally expensive, which makes online simulation difficult. Further work is needed to couple such fundamental models to the popular continuum models in use.

Several researchers have applied optimization to design more efficient electrochemical power sources. Newman and co-workers obtained optimal values of battery design parameters such as electrode thickness and porosity.^{21,24,26,114–117} To simplify the optimization, many of these papers employed models with analytical solutions, which are only available in limiting cases. Battery design optimization using a full order model has been demonstrated by several researchers.^{11,24,26,115,116} Newman and co-workers report the use of Ragone plots for studies regarding the optimization of design parameters, changing one design parameter at a time, such as electrode thickness, while keeping other parameters constant, Ragone plots for different configurations can be obtained. Hundreds of simulations are required when applied current is varied to generate a single curve in a Ragone plot, which is tedious and computationally expensive. An alternative is to simultaneously optimize the battery design parameters and operating conditions such as the charging profile.¹¹ Parameters have been simultaneously optimized for different models and goodness of fits compared based on statistical analysis.¹¹⁸ Parameter estimation has also been used in a discrete approach to analyze and predict capacity fade using a full-order P2D model.^110,111 Golmon et al.¹¹⁹ attempted a multiscale design optimization for improving electrochemical and mechanical performance of the battery by manipulating both micro- and macro-scale design variables such as local porosities, particle radii, and electrode thickness to maximize the capacity of the battery. A surrogate-based framework using global sensitivity analysis has been used to optimize electrode properties.¹²⁰ Simulation results from P2D models have been used to generate approximate reduced-order models for use in global sensitivity analysis and optimization. Rahimian et al.¹² used a single-particle model when computing the optimum charging profile for maximizing the life of battery during cycling. The following section describes the systems engineering tasks of (1) parameter estimation, (2) model-based optimal design, and (3) state estimation that have been applied to lithium-ion batteries.

Parameter estimation is typically formulated as the minimization of the sum-of-squared differences between the model outputs and their experimentally measured values for each cycle i, for example,^121–123

where y_i(t_j) is the measured voltage at time t_j for cycle i, y_model, i(t_j; θ_i) is the voltage computed from the battery model at time t_j for cycle i for the vector of model parameters θ_i (the parameters being estimated from the experimental data), and n_i is the number of time points in cycle i. Solving the optimization [3] is known in the literature as least-squares estimation.^121–123 Many numerical algorithms are available for solving the nonlinear optimization [3], such as the steepest descent, Gauss-Newton, and Levenberg-Marquardt methods.¹²² These iterative methods reduce the sum-of-squared differences between the model outputs and the experimental data points until the error is no longer significantly reduced. More sophisticated Bayesian estimation methods employ the same numerical algorithms but use optimization objectives that take into account prior information on the model parameters.¹²⁴

Battery design parameters such as cell thickness and electrode porosity and operating profiles can be optimized using the same numerical algorithms, for objectives such as maximization of performance (e.g., energy density, life) or minimization of capacity fade and mechanical degradation. These optimizations are solved subject to the model equations and any physical constraints. The optimal estimation of unmeasured states in lithium-ion batteries can also be formulated in terms of a constrained model-based optimization. The optimization objectives, models, and constraints differ for different systems engineering tasks, but can all be written in terms of one general formulation:¹²⁵

where Ψ is the optimization objective,¹²⁶ z(x) is the vector of differential state variables, y(x) is the vector of algebraic variables, u(x) is the vector of control variables, and p is the vector of design parameters. Although there are many numerical methods for solving constrained optimization problems,^127–129 this paper summarizes only control vector parameterization (CVP) as this is the method that is easiest to implement and most commonly used in industrial applications. The CVP method parameterizes the optimization variables, by employing basis functions or discretization, in terms of a finite number of parameters to produce a nonlinear program that can be solved using standard software. First-principles models for lithium-ion batteries tend to be highly stiff, requiring adaptive time-stepping for reasonable computational efficiency.¹⁰⁴ CVP is well suited for optimizations over such models, as CVP incorporates the model equations by calling a user-specified subroutine for simulating the model equations. Any speedup obtained by an adaptive time-stepping for the model equations directly translates into a speedup on the CVP calculations.

More specifically, the control variable u(x) in CVP is parameterized by a finite number of parameters, typically as a polynomial or piecewise-linear function or by partitioning its values over space, and the resulting nonlinear program is solved numerically. Most numerical optimization algorithms utilize an analytically or numerically determined gradient of the optimization objective and constraints to march toward improved values for the optimization variables in the search space. In CVP, as the number of intervals increases, the number of equations increases and makes optimization more computationally expensive. Hence the fastest and most efficient battery model and code for the desired level of accuracy is recommended when applying CVP or any alternative optimization methods.

A discussion of simulating lithium-ion batteries at the systems-level is incomplete without addressing issues pertaining to the estimation of state-of-charge and health of the battery. Designing a tool to predict the life or performance of a battery is an interesting optimization problem with implications on material modifications during the initial battery formulation for a particular application, allowance for making a specific maintenance plan during the course of the life of the battery, and, most importantly, on the cost of the battery. Precise estimations of SOC and SOH are also essential to ensure the safe operation of batteries, that is, preventing the battery from overcharging and thermal runaway.

Some commonly used methods in the industry to monitor the SOC of the battery include monitoring of the cell impedance,^130–133 equivalent circuit analyses,^134,135 techniques based on fuzzy logic,^136,137 or pattern recognition.¹³⁸ Optical and eddy current methods^139,140 are being devised to monitor available capacity in battery systems with flat response surfaces. Based on the algorithm used for estimation, the models used to estimate SOC and SOH can be classified broadly into two categories. Some utilities such as the battery packs used in on-board satellites during the lack of solar energy or cells used in watches follow a routine or pre-programmed load. In such instances, it is possible to develop a degradation model based on a priori testing, knowing the operating conditions and the design parameters of the cell. Such a model does not require frequent updates for the parameters, unless there is a significant change in the operating conditions. In some other applications, such as battery packs used in vehicles, the battery is subjected to a dynamic load that changes as frequently as every few milliseconds. In these cases, the degradation mechanism and hence state of charge or the state of health of the power system depends on the load conditions imposed in the immediate past and it is necessary to monitor the cell on a regular basis. There are some differences between the algorithms used to make life-estimates for the case with the known operating parameters compared to the dynamic-load case. The latter situation is less forgiving in terms of the calculation time, for example. SOC and SOH estimators have been an integral part of battery controllers; however, the estimations have been primarily based on empirical circuit-based models that can fail under abusive or non-ideal operating conditions. Precise estimations of SOC and SOH are very essential for the safe operation of the batteries, in order to prevent them from overcharging and thermal runaway. Santhanagopalan et al.¹⁴¹ reviewed past efforts on the monitoring and estimation of SOC in the literature, and reported an online Kalman filter-based SOC estimation for lithium-ion batteries based on a single-particle model. Klein et al.⁹ recently published state estimation using a reduced order model for a lithium-ion battery. Smith et al.'s¹⁰ analysis of a 1D electrochemical model for a lithium-ion battery indicated that the electrode surface concentration was more easily estimated from the real-time measurements than the electrode bulk concentration. Domenico et al.¹⁴² designed an extended Kalman filter for SOC estimation based on an electrochemical model coupling the average solid active material concentration with the average values of the chemical potentials, electrolyte concentration, and the current density.

Once the battery is manufactured and closed in a sealed case, the battery is discharged (used) according to the requirements of the application. The only variables that can be manipulated during battery operation to make best use of the battery is the charging current profile and operating temperature, which can affect transport and electrochemical rates resulting in modified performance.

Before the battery is sealed, the design variables such as the electrode dimensions, the type of materials, and materials properties such as porosity, active surface area, and microstructure can be selected so as to provide the best possible performance. The resulting battery design can be verified at small scale (e.g., few milli- or micro-Ah batteries) relatively easily in the laboratory, but scaling up to the large-scale batteries required for some industrial applications is challenging.

The physicochemical understanding is incomplete for much of the phenomena that occur inside a battery, such as capacity fade, stress-strain effects, mechanical degradation, and mechanisms for failure due to shocks, defects, and shorts. Much progress has been made in the last twenty years on failure mechanisms, stress-strain models, capacity fade mechanisms involving side reactions, SEI-layer formation, and other phenomena, and studies have been published with the objective of understanding battery operation at the molecular scale, using Kinetic Monte Carlo simulation, molecular dynamics, and density functional theory calculations, and at the mesoscale using population balance models. The molecular-scale models are simulated off-line (that is, not in real-time) and their predictions have been fed to continuum-scale models. A potential future application of molecular- and mesoscale models would be in the real-time prediction of the states of the battery at the small length scales for use in more accurate prediction of the whole battery performance in real time.

Battery models result in multiple DAEs to be simulated with unknown initial conditions while operating for multiple cycles of charge and discharge. For these models adaptive time steppers are usually more than an order of magnitude faster than uniform time-discretization. Several adaptive solvers are available for solutions of DAE models.^143–146 Recently, many easy-to-use ODE solvers have been made available (ode15s, ode15i, etc.) from MATLAB,⁸⁸ "NDsolve" from MATHEMATICA,⁸⁹ and "dsolve" from MAPLE⁹⁰ to solve non-stiff, stiff and moderately stiff DAE models of index-1.

In spite of recent advancements, many of these DAE solvers and initialization routines can fail due to numerical convergence problems during Newton iteration to solve nonlinear equations and singular/ill-conditioned Jacobian matrices resulting from small integration steps. The complexity in battery model simulation is increased by steep variations of the dependent variables (concentrations and potentials) between charging and discharging.

Battery simulations for extended operations, such as during switching from constant-current to constant-potential operations, typically require some form of event detection. The DAEs for battery models increase in complexity and also in number as the accuracy and predictability of models increase. Simulation times for battery models range from milliseconds for empirical circuit-based models to minutes for P2D/P3D models and even days for a multiscale model such as a P2D model coupled with KMC simulation, limiting the options for real-time simulations.

Although much literature exists for capacity fade, SEI-layer formation, and other phenomena, no existing model simulates all of the mechanisms for capacity fade or battery failure. More detailed information is required to sufficiently specify a hypothesized mechanism for a phenomenon before it can implemented in a simulation model, such as

• Which chemical species are formed and consumed in each phase and at the interface between phases?
• What is the physical configuration of each chemical species at the interface between phases (e.g., is a molecule on an electrode surface sticking out into the electrolyte or flat against the surface)?
• How many sites does each molecule on a solid surface cover?

Substantial experimental design efforts are required to answer such questions so the answers can be incorporated into first-principles lithium-ion battery models. Also, most applications using batteries for long-term requirements depend on projections made from model predictions coupled with limited test data; however, the relationship between failure modes during the test conditions and those during actual operating scenarios have not been clearly established – necessitating the tools used in SOC and SOH predictions to be independent of the operating or manufacturing conditions. Quite often in such scenarios, the use of look-up tables limits the confidence in the predictive capabilities of the models.

Conventional degradation models based on extensive testing of batteries under various operating conditions and loads have in general attributed the degradation of battery performance to loss of the active material and loss of lithium that can be cycled. Several detailed models to quantify the signature of these parameters on the aging profile of lithium ion batteries have been presented.^31,147 Other approaches include the use of arbitrary empirical parameters obtained by regressing test data. These models usually interpolate the SOC and the health of the battery based on pre-stored database of information. Such models are widely employed in the industry when sufficient information on the physics of the materials in the batteries is not available – this problem is commonplace among module and pack manufacturers, who assemble the units from cells manufactured by a third-party. It is standard industrial practice to calibrate such models^148,149 since monitoring the evolution of all of the physical parameters such as transport coefficients and the reaction rates within each cell inside the pack is expensive, if not impossible. Network models have also been used to address non-uniform degradation in large format cells.¹⁵⁰

Currently, batteries are charged at constant current until a cutoff potential is reached or a time limit followed by charging at constant potential. However, these charging protocols may result in thermal runaway, leading to under-utilization and possibly even explosions. Given the limited variables that are available for manipulation, it is especially important to make the best utilization of these variables during battery operations. A first-principles battery model can be employed in a dynamic optimization framework to compute a time-varying charging profile that maximizes life, minimizes capacity fade, and improves battery performance.

The determination of an optimized charging profile requires a first-principles model that has high predictive accuracy for a wide range of operating conditions, since charge transfer, reaction kinetics, and diffusion rates may be quite different than in the experiments used in the model development. A first-principles model that describes the battery behavior at the meso- and microscale models would be able to take these effects into account during the dynamic optimization. The application of dynamic optimization to compute an optimal charging profile is illustrated here for a P2D model¹¹ for lithium-ion batteries. The dynamic optimization for a cell was formulated as:

where the optimization objective E is the total energy stored in the cell, V is the voltage obtained from the cell as computed from the first-principles model, i_applied is the applied current to the cell, the charging time t_f was restricted to 1 hr, the maximum allowed voltage was 4.05 V, and the value for V as a function of time. The implementation of dynamic optimization is facilitated by the use of a reformulated model⁶ to compute the optimization objective. The time profiles for the electrolyte concentration at the cathode/current collector interface in Fig. 6 are for three different charging scenarios: (1) conventional charging at constant-current followed by constant-potential charging, (2) constant-current charging at an optimized value obtained by solving the dynamic optimization for a fixed value, and (3) the time-varying charging profile given by Eq. 5. The electrolyte concentration at X = 0 (the cathode/current collector interface) has the highest peak value during dynamically optimized charging, due to its higher initial current. For the chosen chemistry, mass transfer limitations in the electrolyte occur at higher currents. This protocol indicates that to increase the energy density, more energy should be stored at shorter time, albeit causing mass transfer limitations in the electrolyte, and allow the concentration to equilibrate at longer times to ensure longer operability of the battery. During dynamically optimized charging, the electrolyte concentration decreases in the latter part of charging, as lithium-ion transfer slows while more lithium ions are packed into the carbon matrix in the negative electrode. In contrast, after the first 10 minutes the electrolyte concentration is nearly constant during optimized constant-current charging. When a meaningful global objective function was chosen at the system level and robust optimization tool and meaningful models are used, improvements in 'local' battery behavior are observed.

Zoom In Zoom Out Reset image size

The above approach can be considered as a top-down approach, where operating conditions or charging protocols are determined at the system level (battery as a whole), and the system-level behavior is affected by the local mass/charge transfer and reaction effects (Fig. 1) and indirectly manipulates non-measurable internal variables such as the electrolyte concentration or potential or also the solid-phase concentrations as shown schematically in Fig. 6. Physics-based models are required in the dynamic optimization to correctly relate the local effects to the system-level behavior as quantified by the optimization objective. The more detailed and accurate the model, the more optimal 'local' behavior can be determined using the few manipulated variables at the system level.

Note that the SPM model lacks sufficient information on the behavior in the cell to be of much usefulness in the above optimizations. If the first-principles model employed in the optimization includes a high fidelity thermal model, then the localized temperatures in the cell can be included as a constraint in the optimization. A more detailed multiscale model that includes more of the physicochemical phenomena would be needed for optimization of battery operations for very quick charging generally involving rates of 2C or higher.

Another approach that can be used to address the sparsity of manipulated variables is to have the limited number of material properties (manipulated variables) vary spatially. If the electrode architecture is designed to minimize and address every possible local nonideality at the sandwich level, then the system level performance will improve. This can be viewed as the bottom-up approach, where the material properties or electrode architecture, etc. are determined at the electrode level (micro-scale), to produce improved performance at the system level (Fig. 1). Physics-based models are required in the optimization framework to correctly relate the local effects to the system-level behavior as quantified by the optimization objective. For example, consider the minimization of the ohmic resistance at the sandwich level (Fig. 1). Optimization of spatially-uniform porosity reduced the ohmic resistance by 20%, whereas optimization for a spatially-varying profile results in a reduction of 33% (Fig. 7).¹³ Physics-based models are required in the optimization framework to correctly relate the local effects to the system-level behavior as quantified by the optimization objective. Note that improved performance for both solid-phase potential and current are obtained locally, which leads to reduced ohmic resistance across the sandwich, which then relates to improved performance for charge-discharge curves at the system level.

Zoom In Zoom Out Reset image size

To address all the issues in Fig. 1, a more detailed model is required (i.e., moving right along the diagonal in Fig. 3). Possible material properties that can be varied as a function of distance are given in Fig. 2. Note that for particle radius, optimization with the P2D model would yield only the smallest possible radius, but stress-strain models would suggest a different size for mechanical stability.¹¹⁹

The more sophisticated the battery model, more computationally intensive the simulations and optimization. While the value of adding more physicochemical phenomena into battery models is clear, and discussed in more detail below, there is also a need to improve the computational efficiency in the simulation of these models by reformulation or order reduction.

Different simulation methods are effective at different scales (see Fig. 5), which has motivated efforts to combine multiple methods to simulate multiscale systems. Battery models that dynamically couple the molecular- through macro-scale phenomena could have a big impact in understanding and designing lithium-ion batteries. The above continuum models could be coupled with stress-strain models and population balance models to describe the time evolution of the size and shape distribution of particles. Probably the first step would be to couple molecular models with P2D models, to thoroughly validate the coupled simulation algorithms before moving to more computationally expensive 3D continuum models. KMC methods could be combined with P2D models to analyze surface phenomena such as growth of the SEI layer in a detailed manner, similarly as has been applied to other electrochemical systems.^72,151–160 For a 125 × 125 mesh, 2D KMC coupled with P2D model with time steps ranging from nanoseconds to seconds would require simulation times ranging from minutes to hours and even days for a single cycle. Another multiscale coupling that could be useful is to occasionally employ molecular dynamics to update transport parameters in a P2D or 3D model. Molecular dynamics can provide information that cannot be predicted using a P2D or 3D continuum model, but long times cannot be simulated using molecular dynamics, so the combination of the two approaches has the potential to increase fidelity while being computationally feasible.

Zoom In Zoom Out Reset image size

The current literature review suggests that typically researchers have expertise and skills in one or two of the models/methods reported in Fig. 6. If researchers with expertise in different fields collaborate, the task of multiscale model development becomes easier and faster progress can be expected. While black-box approaches are available for some of the methods in Fig. 5, it is strongly recommended that, at least for case studies, hard-coded direct numerical simulation is carried out to enable better understanding of coupling between models at different length and time scales.

The complexities of battery systems have made efficient simulation challenging. The most popular model, the P2D model, is often used because it is derived from well understood kinetic and transport phenomena, but the model results in a large number of highly nonlinear partial differential equations that must be solved numerically. For this reason, researchers have worked to simplify the model though reformulation or reduced order methods to facilitate effective simulation. One method of simplification is to eliminate the radial dependence of the solid phase concentration using a polynomial profile approximation,¹⁸ by representing it using the particle surface concentration and the particle average concentration, both of which are functions of the linear spatial coordinate and time only. This type of volume-averaging^161,162 combined with the polynomial approximation^163,164 has been shown to be accurate for low to medium rates of discharge.^18,165–168 At larger discharge rates, other approaches have been developed to eliminate the radial dependence while maintaining accuracy.^{106,165–168} Approximate solution methods have also been developed for phase-change electrodes, for solid phase diffusion.¹⁶⁹ Recently, discretization in space alone has been used by researchers to reduce the model to a system of DAEs with time as the sole independent variable in order to take advantage of the speed gained by time-adaptive solvers such as DASSL/DASPK.^5,6,144 Such solvers also have the advantage of being capable of detecting events, such as a specific potential cutoff, and running the simulation only to that point.

Complications arise when determining consistent initial conditions for the algebraic equations. Consequently, many good solvers fail to solve DAE models resulting from simulation of battery models.¹⁷⁰ As a result, it is necessary to develop initialization techniques to simulate battery models. Many such methods can be found in the literature for a large number of engineering problems. Methods and solvers specifically focusing on initialization of battery models are also available in the literature.^170,171 Recently, a perturbation approach has been used to efficiently solve for consistent initial conditions for battery models.¹⁷² An alternative continuum representation of the discrete events in the charge/discharge cycle of a battery that does not require initialization between the discrete events of a given cycle or between any two cycles was also proposed.¹⁷³

Proper orthogonal decomposition (POD) has been used to reduce the computational cost in various sets of model equations, by fitting a reduced set of eigenvalues and nodes to obtain a reduced number of equations.⁵ Alternatively, model reformulation techniques have been used to analytically eliminate a number of equations before solving the system.⁶ Other researchers have used orthogonal collocation and finite elements, rather than finite differences, in order to streamline simulations.^75,174,175

For stack and/or thermal modeling of certain battery systems, many attempts have decoupled equations within the developed model.^33–42 This approach breaks up a single large system into multiple, more manageable systems that can be solved independently. This allows the model to be solved quickly, but at the expense of accuracy. For this reason, efficient models that maintain the dynamic online coupling between the thermal and electrochemical behavior, as well as between individual cells in the stack are preferred.

Numerical algorithms for optimization can get stuck in local optima, which can be nontrivial to troubleshoot when the number of optimization parameters is large. This problem can at least be partly addressed using a sequential step-by-step approach (see Fig. 8). For illustration purposes, consider the maximization of the energy density with l_p, l_n, l_s, ɛ_p, and ɛ_n, where l is the thickness of each region and ɛ the porosity (p – positive electrode, s – separator, and n – negative electrode).

• (1)  
  Choose a battery model that can predict the optimization objective and is sensitive to the manipulated variables (e.g., a P2D model).
• (2)  
  Reformulate or reduce the order of the model for efficient simulation. This step has to be judiciously made to ensure that the reduced order model is valid in the range of manipulated variables for optimization.
• (3)  
  Maximize energy density with l_p,
• (4)  
  Using the solution from Step 3 as an initial guess, find optimal values for the two parameters (l_p, ɛ_p).
• (5)  
  Add parameters one by one, in the same manner as in Step 4.
• (6)  
  Arrive at optimal performance with multiple parameters.
• (7)  
  If needed before Step 3, find results with a simpler and less accurate model for a good initial guess.
• (8)  
  After convergence, feed in more sophisticated models (for example, including stress effects) to make sure mechanical stability is not compromised.

Zoom In Zoom Out Reset image size

A similar approach can be used for CVP for dynamic optimization with the total time interval divided as 2, 4, 8, etc. for subsequent optimizations until convergence.

The above algorithm will tend to have better convergence if the parameters in Steps 3–5 are rank ordered from having the largest to the lowest effect on the optimization objective. While advances have been made in the computation of global optima for dynamic optimizations,^104,176 it will be at least a decade before such methods are computationally efficient enough for application to the optimal design of lithium-ion batteries using nontrivial physics-based models.

Fig. 9 shows improved performance at each step of an optimization while successively adding manipulated variables. Capacity matching was used a constraint for the thickness of the negative electrode.

Zoom In Zoom Out Reset image size

Uncertainty quantification methods have been applied to hundreds of different kinds of systems to assess the progress of the development of first-principles models and to assess the confidence in model predictions.^124,177,178 The Monte Carlo method and its many variants for uncertainty quantification are computationally expensive and have become less used over time compared to power series and polynomial chaos expansions. These expansion-based approaches avoid the high computational cost associated with applying the Monte Carlo method or parameter gridding by first computing an approximation to the full simulation model, followed by application of robustness analysis to the approximate model. These expansion-based methods are computationally efficient enough for application to lithium-ion batteries.

For example, consider the discrete estimation of model parameters as a way to track the effects of capacity fade. As of today, capacity fade is attributed to many reasons. This depends upon the chemistry, mode of operation, and size. A wide range of reasons can be linked to transport and kinetic parameters as published elsewhere.^110,111,179 Five effective transport and kinetic parameters were estimated by applying least-squares estimation to the 250 mAh button cells experimental voltage-discharge data. The estimated parameters were the effective diffusion coefficient of lithium ion in the solution phase (D₁), effective diffusion coefficient of lithium in the solid phase for the negative and positive electrodes (D_sn and D_sp), and electrochemical reaction rate constants for the negative and positive electrodes (k_n and k_p).

The effective negative-electrode solid-phase diffusion coefficient and reaction rate constant (D_sn and k_n) were found to decrease monotonically with cycle #, whereas the other three parameters did not follow any particular trend. This suggested that the voltage-discharge curves may not contain sufficient information to accurately estimate the effective values of D₁, D_sp, and k_p, and that the change in the voltage-discharge curves with cycle # could be captured by estimation of only the effective solid-phase diffusion coefficient D_sn and reaction rate constant k_n for the negative electrode. A more detailed analysis suggested that the voltage-discharge curves were very sensitive to the value of the effective solid-phase diffusion coefficient D_sn but weakly sensitive to deviations in the model parameters D₁, D_sp, k_p, and k_n from their nominal values, resulting in large uncertainties in their values when fit to experimental voltage-discharge curves. That the voltage-discharge curves were much sensitive to a negative-electrode parameter (D_sn) suggests that mechanisms for capacity fade in the negative electrode, rather than the electrolyte or positive electrode, were the most important for this battery under these operating conditions.¹¹¹

The overall trend in the variation of model parameters is more reliably assessed by plotting nominal estimates over many cycles. A discrete approach was adopted for the prediction of capacity fade by tracking the change in effective transport and kinetic parameters with cycle number (N). The model parameters D_sn and k_n fit to the experimental data for cycles 25, 100, 200, 300, 400, and 500 were used to predict the remaining battery life based on voltage-discharge curves measured in past cycles. To characterize the degradation in the model parameters, a power law was fit to the estimated parameter values from cycles 25 to 500 as shown in Fig. 10a. By implicitly assuming that the changes in the parameter values are the result of the same mechanism in later cycles, the parameter values for the subsequent cycles were predicted using the power-law expressions. The voltage-discharge curve predicted by this model was in very good agreement with the experimental data at cycle 1000, indicating that the model was able to predict capacity fade as shown in Fig. 10b. Each red dot is an experimental data point and the blue line is the model prediction. It is likely that when more detailed multiscale models become available, there will not be a need to perform fitting and tracking of transport and kinetic parameters with cycles.

Zoom In Zoom Out Reset image size

A rapid update of the parameters usually involves some form of a moving horizon algorithm that estimates the parameters used in the model using an initial set of data points (for example between from the start of the experiment to some interval of time t). These values for the parameters (θ_t) are then used to predict the cell performance for the next few data points (e.g., between times t and t + Δt). The error between the model predictions and the actual data points collected between t and t + Δt is then used to calculate the updated set of parameters θ_t+Δt. This process is repeated at periodic intervals of time or the load. Filtering techniques are commonly employed for on-line estimation;^141,180,181 in most of these algorithms, the measured change in the response is divided between the actual battery response and system noise, based on pre-determined weights assigned to the functions calculating the noise and the battery models. The predicted response for the voltage is compared at the next time step and a correction term is introduced to the weights. More elaborate moving horizon estimates include the influence from several sets of parameters from the past on the current estimates. One example is the use of exponential forgetting functions.^182,183 In this example, the effect of the parameter values θ_t, θ_t+Δt, θ_t+2Δt, etc. on the current estimate θ_t+kΔt is assumed to decay exponentially. The steps are summarized below:

• Step 1: Choose a subset of data points N₀ that end when the parameters need to be updated. Calculate the initial value for the SOC.
• Step 2: Calculate the value of the exponential forgetting function at the end of N₀.
• Step 3: Use the next set of data points N₀+1 to N₁, to calculate the updated values for the parameters in the model equations.
• Step 4: Update the SOC for the next set of data points using the parameter values from the previous step.
• Step 5: Update the exponential forgetting function, based on the data points N₀+1 to N₁, new values for the parameters, and the current value of the SOC.
• Step 6: Repeat Steps 1 to 5 until the end of the data set. This procedure produces a set of values for the SOC updated whenever the error between the model and the experimental data is significant.

The use of such online-tracking algorithms, together with reliable models requiring modest computational effort, greatly reduces the uncertainty associated with assessing the failure mode of the batteries, and can be implemented for a variety of operating conditions. The states of interest are tracked as and when the system operates; the advantages offered by this approach are two-fold: (1) any mitigation scheme can be implemented fairly quickly since the operator does not wait until performing the scheduled capacity checks on the batteries and (2) the proposed methodology does not rely solely on a characterization chart made under lab-scale testing environment, and captures the wear-and-tear imposed by the system on the batteries.