Faster-Than-Real-Time Simulation of Lithium Ion Batteries with Full Spatial and Temporal Resolution

A one-dimensional coupled electrochemical-thermal model of a lithium ion battery with full temporal and normal-to-electrode spatial resolution is presented. Only a single pair of electrodes is considered in the model. It is shown that simulation of a lithium ion battery with the inclusion of detailed transport phenomena and electrochemistry is possible with faster-than-real-time compute times.The governing conservation equations of mass, charge, and energy are discretized using the finite volumemethod and solved using an iterative procedure. The model is first successfully validated against experimental data for both charge and discharge processes in a Li x C6-LiyMn2O4 battery. Finally, it is demonstrated for an arbitrary rapidly changing transient load typical of a hybrid electric vehicle drive cycle. The model is able to predict the cell voltage of a 15-minute drive cycle in less than 12 seconds of compute time on a laptop with a 2.33GHz Intel Pentium 4 processor.


Introduction
Lithium ion batteries, by virtue of having high power density and high open circuit voltage, have evolved as the frontrunner among energy storage systems.Today, they are used in virtually all electronic devices, as well as in hybrid electric vehicles, among other applications.Modeling of such batteries provides a means to better understand the coupled physical and chemical phenomena that occur within the battery.Simulation results can shed light on phenomena such as thermal runaway, transient response of the battery, and effect of external thermal conditions on the performance.Simulations can also lead to design of better thermal management strategies.
The literature is full of several well-established, validated models for the prediction of the performance lithium ion batteries with a single electrode pair.Some of these models account for the detailed electrochemistry and transport phenomena within the battery [1][2][3][4][5][6][7][8][9][10][11][12], while others do not [13,14].Similarly, some models account for spatial nonuniformity of temperature within the battery [3,7,11], while others do not [1, 2, 4-6, 8-10, 12].To the best of the authors' knowledge, the earliest model was developed by Newman and coworkers [1], and later extensively used and improved by White and coworkers [4][5][6][7][8].Although the governing equations proposed by Newman and coworkers are general conservation equations of mass, charge, and energy and are applicable to any geometry, until the late 1990s, most of the modeling effort focused on simple one-dimensional (1D) geometries, in which variations along the electrodes are neglected, and only variations across the electrodes and separator are considered.Although somewhat simplistic, such models have found widespread use because they are still able to capture the coupling of the electrochemistry and transport phenomena within the electrodes.In the past decade or so, with the advent of faster computers, researchers have embarked upon multidimensional modeling of lithium ion batteries [3,15].In particular, it has been brought to light by Wang and coworkers [3] that to fully understand the thermal response of a battery, it is necessary to account for multidimensional effects.
The vast majority of the transport-electrochemistry models presented in the literature have been validated and demonstrated for constant load charge and discharge processes with relatively low charge or discharge rates (up to about 4C).For relatively low charge/discharge rates, diffusion of Li ions International Journal of Electrochemistry within the active particles can keep up with the reactions occurring at the surface of the particle, and therefore, the diffusion resistance to Li ions within the active particles does not play a significant role in the performance of the battery.In such scenarios, the active particles can be treated as lumped masses, and a simplistic subgrid scale model can be used for the solid phase concentration field.Smith and Wang [2] pointed out that while such an approximation is adequate for batteries used in electronic devices, it is inadequate for batteries used in hybrid electric vehicles because they often encounter much higher charge/discharge rates than 4C.
With growing interest in using lithium ion batteries for hybrid electric vehicles, control of these batteries has become an important issue.When the load (current drawn) on a battery changes suddenly, as in an automotive drive cycle, the voltage does not respond instantaneously.The delay in response is caused by the combination of a variety of transport and electrochemical processes within the battery.Ideally, a physics-based mathematical model of the types discussed earlier is able to predict such delays in response.Unfortunately, it is widely believed that such detailed physicsbased models that resolve the battery both temporally and spatially are not efficient enough to be used for control applications.Consequently, so-called equivalent circuit models of batteries that are used either in time or frequency domain have found widespread development and use [16][17][18][19], particularly within the hybrid electric vehicle community [20][21][22].While these models are very efficient, they provide little or no insight into the fundamental processes that occur within the battery.In recent years, efforts have been made to develop reduced-order models that include some of the key physical/chemical aspects while being more amenable to use for control purposes [14,23].
In this paper, it is demonstrated that it is possible to simulate in real time (as warranted by control applications) the performance of a lithium ion battery for realistic hybrid vehicle drive cycles with the inclusion of detailed transport phenomena and electrochemistry and full temporal and normal-to-electrode spatial resolution.

Mathematical Model
A typical lithium ion battery with a single pair of electrodes is comprised of three regions: a negative electrode, a positive electrode, and a separator.A one-dimensional schematic representation is shown in Figure 1.The negative electrode is a porous structure comprised of lithium carbide particles tightly packed together, while the positive electrode is a porous structure comprised of lithium metal oxide particles tightly packed together.The most commonly used metal oxides are LiMn 2 O 4 , LiCoO 2 , or LiFePO 4 .The two electrodes are separated by a separator (or ionic conductor) that allow lithium ions to migrate across, but prevents any electron transport across it.The entire cell is filled with an electrolyte that occupies any space not occupied by the active particles or the binder and filler (not shown in Figure 1).The electrolyte is comprised of a lithium salt dissolved in an organic solvent.
Under normal operating conditions, the following reactions occur at the two electrodes (assuming a manganese oxide positive electrode): During the discharge process, lithium ions that are generated at the negative electrode migrate across the separator to the positive electrode, where they recombine with the electrons that flow through the external circuit to form the metal oxide.The reverse occurs during the charging process.
In the above equations,   is the molar concentration of lithium ions in the electrolyte, while   is the volume fraction of the electrolyte containing part of the porous electrodes and/or separator. eff  is the effective diffusion coefficient of lithium ions within the porous electrodes/separator.It is related to the free-stream diffusion coefficient of the lithium ions through the so-called Bruggeman relationship: where   is the free-stream diffusion coefficient of lithium ions in the electrolyte. + is the so-called transference number and physically represents the charge carried by the lithium ions relative to the solvent due to drift. is the Faraday constant.The electrolyte (ionic) and solid (electronic) phase potentials are denoted by   and   , respectively, while the effective ionic and diffusional conductivities are denoted by  eff and  eff  , respectively.They are also related to their freestream values through the Bruggeman relationship (5).The diffusional conductivity (electrical conductivity in the electrolyte phase due to diffusion of lithium ions) is related to the ionic conductivity (electrical conductivity in the electrolyte phase due to drift or electromigration of lithium ions) by the relationship where  is the universal gas constant and  the absolute temperature.The mean molar activity coefficient of the mixture is denoted by , which is assumed to be constant in this study due to lack of better information.Due to electrochemical reactions at the surface of the active particles, electrons are transferred to the solid (or electronic) phase, while positively charged lithium ions are transferred to the electrolyte phase.This process is represented by a current known as the transfer current and is denoted here by  Li .The transfer current depends on a number of factors, such as the lithium ion concentrations both in the solid and electrolyte phase, the temperature, and the surface overpotential.It is customary to express the rate of generation of this current using the Butler-Volmer kinetics: where   is the active surface area to volume ratio of the electrode and  0 is the so-called exchange current density and is given by where  ,max is the maximum concentration of lithium ions in the solid phase and  ,surf is the concentration of lithium ions at the surface of the active particle, that is, at the solid-electrolyte interface.  and   are kinetic constants that depend on the chemical composition of the electrodes and the electrochemical reactions that occur on the active particles.The exchange current also depends on the kinetic constant  0 .In (7), the quantity  SEI represents resistance due to irreversible film formation at the solid-electrolyte interface.It is generally believed [4-6, 9, 10] that these films form in all batteries during the initial assembly process but grow only if the battery is subjected either to extremely high rates of charge/discharge and/or deep discharge.Aging of the battery is often attributed to growth of this film [10] at the anode, which causes increase in the internal resistance of the battery.In the present model, although there is scope to include the effect of this film (as indicated by ( 7)), due to lack of understanding of the exact mechanism of film formation and growth, this resistance is neglected in this preliminary work.
The overpotential, , drives the electrochemical reactions, and is given by where  is the open circuit (or equilibrium) potential.
The open circuit potential is dependent on the chemical composition of the electrode, its state of charge/discharge, and temperature.Details on how this quantity is calculated are presented later.Equation ( 4) represents the conservation of energy, and conduction is assumed to be the only mode of heat transfer.The average (including both phases) density, specific heat capacity, and thermal conductivity are denoted by ,   , and , respectively.The volumetric heat generation rate within the battery is a net result of three effects: Joule heating, irreversible heat generation, and reversible heat generation and is written as [3] The first term on the right hand side of (10) represents irreversible heating, and the second term represents reversible heating.The third term represents Joule heating due to the transport of electrons within the two electrodes, while the last two terms represent Joule heating due to the transport of ions within the entire battery.The boundary conditions for the governing conservation equations [( 1)-( 4)] are well known [1][2][3] and are summarized in Table 1.

Subgrid Scale Active Particle Model.
Equations ( 1)- (10), along with the boundary conditions shown in Table 1, represent the macroscopic model that can predict the performance (voltage-current) characteristics of a typical lithium ion battery for both charge and discharge cycles.However, the model is not complete in terms of closure.This is because the solid-phase concentration at the surface of the active particles,  ,surf , (as appearing in (7)) is an unknown and requires additional equations to be determined.While the dependent variables in the macroscopic model, namely,   ,   ,   , and , are defined within the entire battery, the solid phase concentration is defined only within the active particles.Since the active particles are much smaller than the characteristic dimensions of the electrodes (active particles have a diameter of less than 1 m, while the electrodes are several tens of ms thick), they cannot be resolved by the same computational mesh that is employed to solve the macroscopic model equations ( 1)- (10).The active particle is essentially at a length scale much smaller than the grid scale of the macroscopic model.Therefore, any equation (or model) used to describe the lithium ion concentration field within the active particles is at the so-called subgrid scale, and is henceforth referred to as the subgrid scale model.In reality, the active particles are bound together within the electrodes using filler and binder materials.As a result, even though they start as spherical nanoparticles, once bound together, they form a complex network better represented by sets of overlapping spheres of various diameters.Doyle et al. [1] proposed treating the active particles as nointeracting perfect spheres of constant radius.In this work we adopt the same assumption, although, with significantly larger computational effort, it is possible to treat more complex shapes, as demonstrated by Kamarajugadda and Mazumder for fuel cell electrodes [24] and by Wang and Sastry [15] for battery electrodes.Under the isolated sphere assumption, and assuming that there is no variation of any quantity along the surface of the sphere (polar and azimuthal uniformity), the concentration of lithium ions within a single active particle can be described by an unsteady diffusion equation written in spherical coordinates with only radial dependence: where   is the diffusion coefficient of lithium ions within the active particle (or in the solid phase).Equation ( 11) requires two boundary conditions, which are written as Equation (12a) represents symmetry at the center of the active particle, while (12b) represents a diffusion-reaction flux balance at the surface of the active particle ( =   ).The active surface area to volume ratio is dependent on the radius of the active particles and their overall volume fraction: where   is the volume fraction of the active particles in the electrode.
Since the transfer current,  Li , is a nonlinear function of the solid phase concentration   , as indicated by ( 6) and ( 7), (11) represents a partial differential equation whose boundary condition (12b) is non-linear.Thus, closed-form analytical solution of ( 11) is not possible.Under the assumption that the mass transport by diffusion within the active particle is much faster than the reaction occurring at its surface (i.e., small mass transport Biot number), the active particle can be treated as a lumped mass, and the solid phase concentration of lithium ions at the surface of the active particle is given by [2] where  ,avg is the average solid-phase concentration of lithium ions within the active particle.Unfortunately, the lumped mass approximation breaks down in situations where the load is large, as demonstrated by Smith and Wang [2].Therefore, for hybrid vehicle applications, it is desirable to solve (11) in its original form, which is done in the present work.Solution of (11) in its original form can be computationally expensive.This is because of two reasons.First, as discussed earlier, its boundary condition is non-linear.Secondly, the equation has to be solved at each computational node (or cell) of the macroscopic model.Thus, if 100 nodes are used in the macroscopic battery model, (11) has to be solved 100 times because for each cell the transfer current, which goes in as an input to the model via the surface boundary condition (12b), may be different.Smith and Wang [2] solve this equation using a finite-element method.In this work, (11) is solved using the finite-difference method with nominally 11 grid points.

Numerical Procedure
The basic procedure used in the present work to discretize the governing partial differential equations is the finite-volume  method.This method is chosen because it guarantees both local and global conservation irrespective of mesh size, while other methods are inherently nonconservative.Furthermore, it is amenable to handling discontinuities in material properties because the method is basically an integral method.For temporal discretization, the backward Euler method is used because it is unconditionally stable.Since all of the governing PDEs have nonlinear source terms, appropriate source term linearization techniques were employed to improve diagonal dominance of the resulting discrete equations sets.Each governing PDE, when discretized and linearized, resulted in a set of tridiagonal equations that were solved using the Thomas algorithm (TDMA solver).The PDEs were solved sequentially.Coupling between the PDEs was addressed by using an outer iteration loop over the five governing PDEs.This iteration loop was also instrumental in addressing nonlinearities in the system.Within each time step, convergence was deemed to have been achieved when the residuals for each conservation equation decreased by at least 4 orders of magnitude.The residual of each PDE was defined as the l2norm (inner product) of the discretized equations without source term linearization.The overall algorithm is depicted in Figure 2. At each time step, once the solution reaches convergence, the results are postprocessed to extract quantities of engineering interest.One of these is the cell voltage, which is given by where  contact is the contact resistance between the electrodes and the current collectors and  collector is the collector plate surface area.The applied (drawn) current or load is denoted by  app .

Model Validation.
A validation study was first undertaken.In this study, the entire charge and discharge cycle of a Li  C 6 -Li  Mn 2 O 4 battery was simulated.The battery considered here has a nominal rating of 6 Ah, and a charge/discharge rate of 1C was used.Experimental data is available for the same set of conditions.Other necessary geometric parameters, material properties, and operating conditions are summarized in Table 2.The data presented in Table 2 are extracted from Smith and Wang [2], who also provide the experimental data used for the validation study.The only data that were not directly extracted out of Smith and Wang [2] are the two kinetic constants.The values reported in Table 2 for these constants were estimated to obtain an exchange current density of 36 A/m 2 in the negative electrode and 26 A/m 2 in the positive electrode, as reported by Smith and Wang [2].
In addition to the data summarized in Table 2, Smith and Wang [2] also provided curve-fits to experimental measurements for the ionic conductivity and the open circuit potential.The free stream ionic conductivity (electrical conductivity of the electrolyte phase due to drift of ions) in S/m is given by [2]  = 15.8 × 10 −4   exp [0.85(   1000 ) where  and  are the stoichiometries in the negative and positive electrodes, respectively.The open circuit potentials computed using (17) are in Volts.For the external contact resistance,  contact , a value of 20 Ω cm 2 was used.Although the present model has the ability to predict the temperature field, for the validation study, a constant temperature of 288 K was used.
In order to compute the voltage characteristics of the entire charge/discharge cycle, the calculations were performed for a net duration of 3800 seconds, with a time step of 1 second.Within each time step, it took approximately 50 iterations to achieve convergence by 4 orders of magnitude.A typical convergence plot is shown in Figure 3.
Figure 4 shows the comparison of the computed cell voltage with measured voltage for the entire charge-discharge cycle with a charge/discharge rate of 1C.The computed results agree quite well with the experimental data except toward the extreme end of the discharge cycle, where the predicted voltage drop occurs slightly later than the measured voltage drop and at a much sharper rate.This discrepancy could be due to a number of factors including uncertainties in material properties and kinetic constants, neglecting the formation and growth of the SEI film, and thermal effects.
Since the present model resolves the battery both temporally and spatially, it is instructive to carefully examine the evolution of some of the dependent variables.Figure 5 shows spatiotemporal evolution of lithium ion concentration in the electrolyte phase, the electrolyte phase electric potential, and the overpotential.It is seen that with time the lithium ion concentration distribution becomes more lopsided, while the electric potential in the electrolyte phase is fairly uniform.This is because the ionic conductivity of the entire cell is quite high, resulting in a small potential variation.The overpotential, though quite small at 1C discharge rate, shows significant spatial variation.At the negative electrode, towards the end of the discharge cycle, it rises rapidly.This is consistent with the sharp drop in the electrolyte phase concentration near the negative-separator interface.The simulation of an entire charge cycle or an entire discharge cycle (3800 time steps with a time step of 1 s) required about 42 seconds on a laptop with a 2.3 GHz Intel Pentium 4 processor.Of these 42 seconds, about 15 seconds after used to write data to files for postprocessing, implying that a 3800-second cycle was simulated in less than 30 seconds of real time.

Demonstration for HEV Driving Cycle.
In order to demonstrate the model for real-time control applications, a driving cycle, typical of that of a hybrid electric vehicle (HEV) was simulated.Current (load) pulses as large as 10C and lasting over only 10 seconds was used.The load cycle was commenced with the battery at 100% state of charge.Figure 6 shows the computed voltage response of the battery subjected to the prescribed pulse load.The predicted voltage shows that the numerical algorithm is able to cope up with the rapidly changing load.The jumps in cell voltage when the current is reversed (i.e., when the operation mode is changed from charge to discharge or vice versa) are caused by the differences in the output voltage for the same current during charge versus discharge, as observed in Figure 4.The 900second drive cycle simulation required about 12 seconds of CPU time (which is almost the same as wall-clock time) on a laptop with a 2.3 GHz Intel Pentium processor, out of which about 3 seconds were used for writing files for postprocessing.This implies that the execution time of the current model is about 100 times faster than real time.It is to be noted that the aforestated performance was recorded when the code was written of Fortran90.When the same code was rewritten using MATLAB and executed, the performance deteriorated significantly, and compute times of only about 5-8 times faster (as opposed to 100 times faster) than real time were achieved.

Thermal Effects.
In order to demonstrate the model's ability to simulate thermal effects, the same load cycle considered in the preceding study was now simulated but with the inclusion of the energy equation.Prior to performing the simulation, however, the thermal properties of the materials had to be estimated.Table 2 already shows the thermal properties used for the simulations.These were extracted from Srinivasan and Wang [3].These authors also provided curve-fits for the rate of change of the open circuit potential with temperature (as appearing in (10)) for both Li  C 6 and Li  Mn 2 O 4 electrodes.During a discharge cycle, reaction (R1) is endothermic, while reaction (R2) is exothermic.The net effect of the two reactions, however, is to generate heat irrespective of charge or discharge.This is clearly evident if the battery is operated under adiabatic conditions.Figure 7 shows the temperature evolution of the battery under such conditions.A 4.5 K rise in temperature is observed during the 15-minute hybrid drive cycle considered earlier.Although not shown here, even with Newton cooling at one of the ends, the temperature distribution was found to be very uniform because the battery is only 111 m thick and the resulting conduction resistance across it is very small.Therefore, even though the heat generation within the battery, as shown by (10), is not uniform, a very small spatial variation of temperature is observed within the battery because the heat spreads very quickly over a distance of 111 m.To truly understand thermal effects, multidimensional modeling is warranted.

Summary and Conclusions
A model for predicting the voltage characteristics of a lithium ion battery subjected to a prescribed load was developed.The model completely resolves both time and space across the battery and is, therefore, capable of capturing any spatial nonuniformity of relevant quantities such as concentrations, electric potentials, and temperature within the battery.Although the model is validated and demonstrated here for a Li  C 6 -Li  Mn 2 O 4 battery, in principle, it is general enough to be applicable to any lithium ion battery provided that the material properties and kinetics are appropriately characterized and known.The model was validated against experimental measurements for the entire charge and discharge cycle of a Li  C 6 -Li  Mn 2 O 4 battery and was able to accurately predict the voltage characteristics.It was then successfully demonstrated for a load input typical of that of an HEV automotive drive cycle.
The computational efficiency of the model was noteworthy.A 900-second drive cycle was computed in less than 12 seconds of CPU time on a laptop, indicating that the model can be integrated with real-time control algorithms.
Future work will involve incorporation of appropriate submodels to describe aging of the battery so that the life cycle of a battery can be predicted.Extension of the model to complex three-dimensional geometries is currently underway.

Figure 1 :
Figure 1: Schematic one-dimensional representation of a lithium ion battery with a single pair of electrodes under discharge conditions.

Figure 2 :
Figure 2: Algorithm employed for solution of the governing equations.

Figure 3 :
Figure 3: Convergence behavior at the tenth time step.

Figure 4 :
Figure 4: Computed cell voltage versus measured cell voltage for the entire charge-discharge cycle of a Li  C 6 -Li  Mn 2 O 4 battery with 1C charge/discharge rate: (a) charge, (b) discharge.

Figure 5 :
Figure 5: Spatiotemporal evolution of key quantities during the beginning, middle and end of the discharge process with 1C rate of discharge: (a) electrolyte phase concentration, (b) electrolyte phase potential, and (c) overpotential.

Figure 6 :
Figure 6: Computed voltage response of a 6 Ah battery subjected to an arbitrary load (current) over a 15-minute period.The prescribed load (current) is shown by the dotted red line, while the predicted voltage is shown by the solid green line.Positive current denotes discharge.

Figure 7 :
Figure 7: Spatiotemporal evolution of temperature within the battery under adiabatic conditions during a 15-minute hybrid drive cycle.

Table 1 :
Summary of boundary conditions required for the solution of the governing conservation equations.
Initial conditions for C e , C s and T Guess C e , C s ,  e ,  s , and T of whole computational domain Solve linearized equation for C e

Table 2 :
Values of various parameters for the Li  C 6 -Li  Mn 2 O 4 battery simulated in the present work.