Modeling of Li-Ion Cells for Fast Simulation of High C-Rate and Low Temperature Operations

Evaluation of Li-ion cell performance and life requires the ability to predict behavior at extreme conditions, such as low temperatures and high C-rates. Most electrochemical models assuming constant electrolyte diffusion properties fail to accurately predict the electrode and electrolyte potential at such conditions. This study presents a physics-based Extended Single Particle Model (ESPM) designed specifically for accurately predicting the behavior of a Li-ion cell at extreme conditions, incorporating concentrationdependent properties in the electrolyte diffusion dynamics. Since the proposed model aims at supporting long-term simulation, virtual design and optimization studies, minimization of the computational complexity is achieved through analytical Model Order Reduction (MOR) based on a Galerkin projection method. Results show that the implementation of the concentration-dependent diffusion properties leads to significant improvement of model accuracy at extreme conditions. The Reduced Order Model (ROM) can be simulated significantly faster than numerical methods with no loss of accuracy, supporting simulation of long-term usage cycles (10-year) and remaining useful life calculations. © The Author(s) 2016. Published by ECS. This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License (CC BY, http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse of the work in any medium, provided the original work is properly cited. [DOI: 10.1149/2.0761605jes] All rights reserved.

Lithium ion batteries are considered the state of the art for energy storage in electric and hybrid vehicles. When batteries operates at low temperature conditions, for instance during a cold start of an electric vehicle, the highly reduced ionic conductivity and diffusivity lead to the formation of large gradients in the electrolyte concentration within the cell domain. 1,2 Large concentration gradients can also be established at high C-rate conditions, which typically occur when the battery is subject to fast charge/discharge current loads in cold weather. Therefore, it is important to accurately predict the cell electrochemical behavior and terminal voltage at such conditions, as performing fast charging procedures at low temperatures could severely shorten the cycle life due to lithium deposition on the anode. 1,3,4 Electrochemical battery models based on first principles have shown the ability to accurately predict the concentration dynamics and terminal voltage. [5][6][7][8][9][10][11] In particular, the pseudo two-dimensional model (or the P2D model) based on Porous Electrode Theory 5 and Single Particle Model 12 (SPM) have been widely used for modeling the dynamic response of Lithium ion cells to variable input current profiles. Recently, Extended Single Particle Models (ESPM) [13][14][15][16][17] incorporating the electrolyte dynamics have been developed to capture the cell behavior under high C-rates conditions, albeit with simpler mathematical structure than P2D models.
While electrochemical models generally predict well the cell terminal voltage as function of current and temperature, they are characterized by the presence of coupled PDEs and nonlinear algebraic equations, increasing the mathematical complexity and leading to computational challenges when applied to fast simulation and to estimation or control design. For this reason, several approximation methods based on Model Order Reduction (MOR) techniques have been recently applied to electrochemical models using various techniques, including Padé approximation, 18 residue grouping, 19 polynomial approximation 20 and Proper Orthogonal Decomposition 21 (POD). In recent times, the Galerkin Projection method has also been applied to analytically approximate the diffusion PDEs describing the Lithium transport in the solid and electrolyte phases. [22][23][24][25][26][27][28] The main advantage of the Galerkin projection method is that it is able to predict the concentration profile in the whole cell domain while preserve the physical nature of model parameters, which allows to use the resulting Reduced Order Model (ROM) for state and parameter estimation. In addition, the accuracy is comparable to numerical MOR methods such as POD and Balanced Truncation, at a fraction of the computational cost. 29 When modeling the cell behavior at high C-rate and low temperature conditions, the concentration dependence of some electrochemical parameters, such as the electrolyte diffusion coefficient, should be considered to characterize the cell performance accurately. In fact, assuming constant diffusion coefficient under the above circumstances would lead to errors of 20% or more in predicting in the Lithium concentration at the boundaries and, ultimately, the cell terminal voltage. 30 The presence of concentration-dependent and temperature-dependent diffusion properties results into nonlinear diffusion PDEs, significantly increasing the numerical complexity in the resulting electrochemical models and preventing the use of MOR techniques for linear or quasi-linear systems, such as Balanced Truncation and Padé approximation.
In Ref. 15, the concentration-dependent diffusion coefficient is implemented and the nonlinear electrolyte diffusion PDE is reduced to a Differential Algebraic Equations (DAE) system and is solved by using the electrolyte average concentration. Ref. 23 incorporated the concentration-dependent diffusion coefficient in the pseudo 2D electrochemical-thermal models for performance simulation. A reduced-order nonlinear DAE system was obtained using transformation and orthogonal collocation reformulation techniques and was solved by adaptive solvers in time. Alternatively, the case of concentration-dependent diffusivity can be handled directly with discretization methods, albeit usually producing less efficient models not suitable for long-term life simulation, control and estimation applications. In Ref. 31, the orthogonal collocation on finite elements method and POD was applied to obtain a reduced-order P2D cell model. However, the resulting numerical models were generated for a particular model parameter set. As the parameters change due to battery aging, these numerical models might fail to capture the performance degradation. This limits the capability of these numerical models for long-term life simulation and estimation. This paper proposes a novel physics-based electrochemical model designed specifically for accurately predicting the behavior of a Li-ion cell at low temperature and high C-rate conditions. Under the framework of an ESPM, this work shows how to incorporate concentrationdependent and temperature-dependent properties in the PDEs describing the Lithium transport in the electrode and electrolyte region, resulting into an improved characterization of the cell behavior at extreme conditions. Since the proposed model is tailored to support virtual design and optimization studies of Li-ion battery packs for electric and hybrid vehicles, minimization of the computational complexity with no loss of accuracy or physical consistency are critical requirements. This is here achieved through a MOR method based on a Galerkin projection A667 of the nonlinear PDEs onto a functional subspace. The reduced order model can be simulated much faster than real-time to support simulation of long-term usage cycles (10-year) and remaining useful life calculations.
This paper is structured as follows. The next section introduces the nonlinear ESPM obtained by incorporating the concentrationdependent diffusion ROM in the electrolyte phase. The temperature dependence of the electrochemical parameters and variations of the open-circuit voltage (OCV) at different temperatures are also implemented into the nonlinear ESPM. The model is then validated against experimental data at different C-rate and temperature conditions. Next, to illustrate the impact of the concentration-dependent diffusion coefficient at low temperature and high C-rate conditions, the nonlinear ESPM is also compared with a linear ESPM that assumes constant diffusion properties in the electrolyte solution. Finally, the influence of the variation in OCV at low temperatures on model accuracy is also investigated in this work.

Model Development
The starting point of this work is the Extended Single Particle Model (ESPM), which approximates a porous electrode with a set of decoupled PDEs. [13][14][15][16][17] A schematic of the cell model is illustrated in Fig. 1.
The intercalation in the solid phase is modeled according to the Fick's law, leading to: with boundary conditions: where the subscript k = p indicates the positive electrode and k = n the negative electrode. The other symbols can be found in the List of symbols section. The intercalation current density J k is assumed to be uniform, and can be expressed as: where the subscript s indicates the separator.
To improve the accuracy in predicting the diffusion dynamics in the electrolyte solution, a concentration dependent diffusion coefficient is incorporated into the corresponding PDE: [4] with zero flux conditions at the boundaries: Since ε e,k , D e,k (c e ) and J have different values in the positive electrode, the separator and the negative electrode, additional boundary conditions should also be imposed to satisfy the continuity of the concentration and flux at the solid/liquid interface: The concentration in solid and liquid phase are then coupled by the kinetic overpotential η k , given by the Bulter-Volmer equation: The exchange current density i 0,k is defined as: where K is the kinetic rate constant, c s,max,k is the saturation concentration of the electrode and c s,sur f,k is the electrode surface concentration. The potential in the liquid phase φ e is related to the liquid concentration in the following way: The integration of Eq. 9 over x yields: The diffusional conductivity κ ef f D can be expressed as: where β is the activity coefficient and is assumed to be a constant in this work. The block diagram of the cell model is shown in Fig. 2, from which the cell output voltage can be obtained as: where the open-circuit potentials U k (k = p, n) are determined as a function of the surface concentrations in the solid phase. The parameters used in the model can be found in Table I. It should be noted that the model equations described above are in an isothermal condition. The model parameters, including temperature-dependent parameters such as D s, p , D s,n , K p , K n , D e , κ and κ ef f D 30,32,33 are initially calibrated at 25 • C to build a reference for benchmarking purposes. The Open Circuit Potential (OCP) curves relative to Li/Li + for the positive and negative electrodes are obtained experimentally from Refs. 34 and 35, and are shown in Fig. 3. In this study, the upper voltage limit of the cell is designed as 4.1 V, and the initial lithiation y in Li y NiMnCo (NMC) is chosen as 0.5 to match this upper voltage limit.

Concentration-Dependent Parameters
Key to correctly predict the cell behavior at high C-rate conditions is the inclusion of a concentration dependent diffusion coefficient in the electrolyte solution. To this extent, experimental data and empirical correlations have been previously reported in Ref. 30 for LiPF 6 electrolyte solutions at different temperature conditions. The main results of this study are shown in Fig. 4.
From prior simulations, the maximum values of the concentration gradients in the electrolyte solutions are typically limited between 300 mol/m 3 and 3500 mol/m 3 , 15,21,22,25,30 which allows for   approximating the functional dependence of the diffusion coefficient from the Lithium concentration with a linear expression: where h 1 and h 2 are the polynomial coefficients of the linear fit. The effective diffusion coefficient is then obtained by considering porosity and tortuosity: , ε e denotes the porosity and brugg is the Bruggman coefficient. Note that ε e,k , brugg k are spatially dependent and are defined in Table I.

Temperature-Dependent Parameters
It has been shown that the performance of batteries has a strong dependence on the temperature. In the solid phase, the temperature dependent parameters such as solid diffusion coefficients D s, p , D s,n and the kinetic rate constant K k can be modeled based on the Arrhenius equation: where denotes the temperature dependent electrochemical parameter, re f is the reference value of the parameter and E act is the activation energy. The values of re f and E act for D s, p , D s,n , K p and K n are documented in Table II. The temperature dependent parameters in the liquid phase are electrolyte diffusion coefficient D e , electrolyte ionic conductivity κ and electrolyte diffusional conductivity κ

Temperature
electrolyte concentration of 1000 mol/m 3 . The linear fitted electrolyte diffusion coefficient in Fig. 4 will be used in the nonlinear electrolyte diffusion model, and the resulting temperature dependent coefficients for the linear fit are listed in Table III.
The temperature dependence of electrolyte ionic conductivity κ and electrolyte diffusional conductivity κ ef f D are also extracted from Ref. 30 and can be expressed as: where [18] In this work, κ and κ ef f D are assumed concentration-independent, and only temperature-dependent.

Temperature Dependence of Open-Circuit Potential
The open-circuit potential of the NMC-based positive electrode and graphite-based negative electrode and cell capacity used in this study are initially calibrated at 25 • C. However, it has been shown that the cell OCV curves shift slightly when the cell is operated at different temperatures due to entropy effects. 36,37 Therefore, the variation of cell OCV with temperature should also be considered in the nonlinear ESPM to accurately characterize the cell behavior at different temperature conditions.
The variations of OCPs of NMC and graphite were investigated as part of Ref. 37. Since the cell capacity in this work is slightly different than the one reported in Ref. 37, the U/ T curve as a function of SOC in Fig. 5 was modified based on the ratio of negative and positive Where the OCP model is integrated into the nonlinear ESPM, the reference temperature was set to 25 • C, so that U at different temperature conditions can be directly calculated as a function of SOC. Finally, the "corrected" OCV for the full cell can be written as: c s,sur f, p,k (t)/c s,max,k , and U/ T combines the variations of OCPs for the anode and cathode using the interpolation function described in Ref. 37.

Model Order Reduction of the Nonlinear Diffusion PDEs
The Galerkin projection method is used to approximate the Lithium diffusion dynamics in the electrolyte solution. Application of this method with constant diffusion coefficients has already been found in Ref. 22. In Ref. 23, the concentration-dependent diffusion PDE is reduced to a system of DAEs. Boundary conditions are satisfied by adding linear and/or quadratic terms to a Fourier series expansion in basis functions. Recently, this method has also been applied along with coordinate transformation to deal with the non-zero, time-dependent boundary conditions in the solid phase. [25][26][27][28] In this paper, a general, "top-down" process for projection-based model order reduction is proposed to analytically approximate the nonlinear diffusion PDE into a low-order ODE system.
First, trial solutions for the electrolyte concentration in the positive electrode, the separator and the negative electrode are firstly assumed as:ĉ [20] where ψ i (x k ) are basis functions, p i (t), s i (t) and n i (t) are the timevariant coefficients of the basis functions for each region, N is the number of the basis functions, and x k is the dimensionless length for each region: The diffusion PDE and boundary conditions using dimensionless coordinates can be written as: The last two equations are a result of the continuity of concentration and temperature at the solid/liquid interfaces and consequently the effective diffusion coefficient D e,k in Eq. 14.
In this study, the shifted and normalized Legendre polynomials 24 are chosen as basis functions and an order of truncation N = 3 is chosen here for illustration purposes: Note that this choice of trial solutions should satisfy the boundary conditions in Eq. 23. Substitution of Eqs. 20 and 24 into Eq. 23 yields six algebraic equations with respect to the time-variant coefficients p i (t), s i (t) and n i (t): The residuals are obtained by substituting the trial solution Eq. 20 into Eq. 22: Substituting Eqs. 20 and 14 into Eq. 26 gives: ) unless CC License in place (see abstract Here k i is used to represent the coefficients p i (t), s i (t) and n i (t) for clarifications. Applying the Galerkin projection method to Eq. 27 gives: Note that the inner product ψ i (x k ), R c e,k (x k , t) , (i = 1, 2, 3) is calculated in the three domains, as the model parameters such as porosity and Bruggman coefficient are spatial dependent in the positive electrode, the separator and the negative electrode.
Then nine nonlinear ODEs are obtained by substituting Eq. 27 into Eq. 28 and calculating the integrals with respect to x k : 10,k I [29] where k = p, s and n, and Here the nonlinear terms k 2 i and k i k j (i, j = 1, 2, 3) are preserved from Eq. 27. All the coefficients C ψ i 1,k , . . . , C ψ i 10,k in Eq. 30 can be calculated analytically in the pre-processing phase. Also, the orthogonality properties of the Legendre polynomials can be applied to simplify terms in Eq. 30. Once the calculation is complete, each coefficient is an analytical function of model parameters in Table I. It is also worth noting that if the diffusion coefficient is assumed to be constant, the coefficients C ψ i 1,k , . . . , C ψ i 6,k will equal to zero, leading to a linear ODE system.
To obtain the electrolyte concentration, nine equations are needed to solve for the nine coefficients in Eq. 20. Since the six algebraic equations in Eq. 25 must be considered to satisfy all the boundary conditions, the three remaining equations can be obtained by choosing p 1 (t), s 1 (t) and n 1 (t) as three independent variables in Eq. 29. The system is now reformulated as: [32] where the expression for θ( p 1 (t), s 1 (t), n 1 (t)) can be found in Eq. 33. Once p 1 (t), s 1 (t), n 1 (t) are solved by Eq. 31, the other six coefficients can be calculated using Eq. 33. Finally, the electrolyte concentration in each region can be obtained using Eq. 20.
x [36] where N p , N s , N n and N represents the number of nodes in the positive electrode, the separator, negative electrode and the whole cell domain, respectively, and c e (i, t) is the concentration at the ith node (i = 1, 2, · · · N ). It should be noted that the continuity of the concentration at the solid/liquid interface are automatically satisfied by using FDM. Therefore, the total number of ODEs in Eq. 35 is N − 4 since four algebraic equations are generated by the boundary conditions. Numerical simulations were run in MATLAB on a PC with a 2.60 GHz i5-3320M CPU and 4 GB RAM. The presence of the high-order terms in Eq. 31 and Eq. 35 leads to a nonlinear model that exhibits numerical stiffness, hence the variable-step ODE15s solver was used to solve the ROM and FDM model. A maximum step length of 1 s was used for the solver. Figure 6 compares simulation results where the electrolyte concentration profiles at the end of 1 C discharge at 25 • C obtained with the ROM and FDM with different numbers of nodes are compared. It can be observed that the FDM model becomes insensitive to the number of nodes after about 500 elements. In addition, even if the nonlinear diffusion ROM only uses three states [ p 1 (t), s 1 (t), n 1 (t)] in this study, it agrees very well with the high-order FDM solution. More importantly, the CPU time for the ROM is almost 6 times faster than the FDM model, as shown in Table IV. It can be also noted that the 3rd-order ROM only slightly outperforms the 100-node FDM in terms of computational speed due to the presence of a much larger number of nonlinear terms in Eqs. 31 and 32.
The nonlinear diffusion ROM is then compared with the FDM model using 500 nodes at different C rates and different temperature conditions. As shown in Figs. 7 and 8, the diffusion ROM has a good agreement with the FDM model at different conditions, and only marginal differences in the concentration at the boundaries can be observed at high C-rate conditions.

Influence of the Concentration-Dependent Diffusion Coefficient on Model Accuracy at Low Temperature and High C-Rates
The influence of the concentration-dependent diffusion coefficient on the model accuracy is investigated in this section. For comparison purposes, a simplified electrolyte diffusion model is used as benchmark, assuming a constant diffusion coefficient. Applying the Galerkin projection method to this PDE results into a linear ROM. 25 Figure 9 compares the concentration distribution in the electrolyte obtained by the linear and nonlinear electrolyte diffusion ROMs. The concentration-dependent diffusion leads to significant differences in the prediction of the Lithium concentration within the cell domain. Larger concentration gradients in the electrolyte solution are also observed in the distribution map obtained by nonlinear ROM as expected. For example, as ions accumulate at the negative electrode during charge, the concentration in the negative electrode domain gradually increases, and the diffusion process will be decelerated as a result. As shown in Fig. 9, the implementation of the concentration-dependent coefficient leads to the prediction of larger concentration gradients at A comparison of the electrolyte overpotential predicted by the linear and nonlinear diffusion ROMs is shown in Fig. 10. The electrolyte overpotential is obtained by calculating the potential difference between the cell boundaries using Eqs. 9 and 10. Note that the potential at x = 0 is set at 0 V since the only the potential differences are relevant to the calculation of the cell terminal voltage. 14 As shown in Fig. 10, the electrolyte overpotential predicted by the nonlinear ESPM is higher due to larger concentration gradients at the cell boundaries. Then the electrolyte overpotential is used to calculate the cell terminal voltage following Eq. 12.
By incorporating the nonlinear and linear diffusion ROMs in the electrolyte phase, a nonlinear ESPM and a linear ESPM are developed separately for comparison purposes. The spherical diffusion PDEs in the positive and negative electrodes in Eq. 1 are also reduced via Galerkin projection, and implemented into both the nonlinear and linear ESPMs. Details on this step can be found in Ref. 25. A prismatic 5 Ah cell with 3.7 V nominal voltage was considered in this study. The battery was cycled under different conditions, including pulse charge/discharge current inputs and a US06 driving cycle. All the input current profiles are described in Fig. 5. Current cycles were applied with an Arbin BT-2000 battery tester and temperature control was achieved by a Thermotron Industries SE-300-2-2 environmental chamber. The chamber temperature was set to constant values of -10 • C, 0 • C and 25 • C during the cycling experiments. Figure 11 compares the voltage predicted by the nonlinear ESPM with the experimental data at different test conditions. As shown in Figs. 11a, 11b and 11c, the battery is subject to pulse charge/discharge current inputs at -10 • C, 0 • C and 25 • C respectively, and the predicted voltage agrees well with the measured voltage. To further validate the nonlinear ESPM under dynamic input conditions with high charging/discharging rate, an additional experimental dataset of US06 driving cycle are used to compare with the simulated voltage in Fig.  11d. The agreement of the predicted voltage with the measurements indicates that the model is still able to maintain a high level of accuracy even though the magnitude of the current input reaches to as high as 150 A (30 C-rate).
The prediction errors produced by the linear ESPM and nonlinear ESPM at different test conditions are also shown in Fig. 11. In addition, RMSEs in voltage are calculated and summarized in Table V. Both the error plots and RMSEs show that the voltage prediction errors are reduced at each temperature and current input condition by the nonlinear ESPM. It is also worth noting that the implementation of the concentration-dependent diffusion coefficient in the nonlinear ESPM results in significant improvement of model accuracy at low temperature conditions. As illustrated in Fig. 9 and Fig. 10, larger electrolyte overpotentials are induced by larger concentration gradients at low temperatures, and consequently contribute to the terminal voltage predictions.
The computation time required by the linear ESPM and nonlinear ESPM for each simulation are also listed in Table V. Although the nonlinear ESPM is slower than the linear ESPM due to the presence of the nonlinear ODEs, the simulation time for the nonlinear ESPM is still at least 330 times faster than real time. Figure 12 illustrates the influence of the implementation of the OCV variations on model accuracy. Another nonlinear ESPM is developed here which does not consider OCV variations at different temperatures. The voltages predicted by the two models are compared together with the experimental data and the corresponding errors are shown in Fig. 12. As the reference temperature is set to 25 • C, a comparison of the model error is only shown at -10 • C and 0 • C. When U/ T is involved, the predicted voltage agrees better with the measurements. The voltage RMSEs of the nonlinear ESPM without U/ T are 34.8 and 14.2 mV respectively at -10 and 0 • C, which are slightly higher than the corresponding RMSEs in Table V. Particularly, model errors are almost reduced to 0 at the relaxation part (from 1.0 to 1.7 hours, or after 6.4 hours in Fig. 12a; from 0.8 to 1.6 hours, or after 6.1 hours in Fig. 12b), which indicates that the improved nonlinear ESPM predicts more accurate OCV of the battery. However, the decrease in the RMSEs is not significant as expected since the magnitude of U/ T in Fig. 5 is relatively small. Nevertheless, it  indicates the integration of change of OCV at different temperature conditions does provide better model accuracy than otherwise.

Conclusions
This paper describes a nonlinear ESPM developed specifically for extreme operating conditions by incorporating a concentrationdependent diffusion dynamics in the electrolyte solution. The computational complexity of the nonlinear diffusion PDE is reduced through a Galerkin projection of the PDE to a low-dimensional functional subspace. The resulting electrolyte diffusion ROM is validated against the numerical solution via FDM at different operating conditions. Simulation results show that the accuracy of the nonlinear diffusion ROM is comparable to FDM at a fraction of the computational cost. Then the nonlinear ESPM is benchmarked against with a linear ESPM assuming constant electrolyte diffusivity, as well as experimental data at different temperature and high C-rate conditions. Results show that the nonlinear ESPM is more accurate at all temperature and current input conditions due to the improved consistency of the model with the physical behavior of the system. The nonlinear ESPM with improved accuracy can be used to support simulation of long-term usage cycles (10-year) and remaining usable life calculations at high rates and low temperatures.