Nonlinear State-Variable Method for Solving Physics-Based Li-Ion Cell Model with High-Frequency Inputs

A nonlinear state-variable method is presented and used to solve the pseudo-2D (P2D) Li-ion cell model under high-frequency input current and temperature signals. The physics-based governing equations are formulated into a nonlinear state variable method (NSVM)

Physics-based Li-ion cell models are useful tools for analyzing the battery cell performances, characterizing material properties, supporting cell design, and developing battery management systems (BMS). 1 The development of macroscopic, physics-based models of lithium ion cells and batteries began with a relatively simple solid phase lithium ion diffusion model, extended from Atlung et al.'s 2 work in 1979 and named by Santhanagopalan et al. 3 5 in their work, a dual-foil Li-ion cell was represented by a pseudo-two-dimensional domain; and therefore, this model is usually referred as the "pseudo-2D model (P2D)". 3 This model contains several coupled transport phenomena described by nonlinear partial differential equations (PDEs). Due to the numerical complexity of these models, various reduced-order models (ROMs) have been proposed as substitutes for either the single particle (SP) or the P2D full-order model (FOM) in practical online simulation applications. [6][7][8][9][10][11][12][13][14][15][16] Also, ROMs have been published to simulate capacity fading 17,18 and for use in parameter estimations. [19][20][21][22][23][24] The single particle model dates back to 1979, when Atlung et al. 2 presented a relatively simple cell model for the Li/TiS 2 couple. Their model includes the diffusion of lithium ions in different shaped particles of TiS 2 , which they called a solid solution. They used their model to predict several aspects of this cell including the specific energy. Their model was used later by Haran et al. 25 to model the impedance of a spherical metal hydride particle. After that, Subramanian et al. 26 extended Atlung et al.'s model to include a variable diffusion coefficient. Also, Subramanian et al. 27 developed an approximate solution for the concentration distribution in a spherical particle, which reduced the computation time needed to simulate cycling of a cell. Next, Ning and Popov 28 used Atlung et al.'s model to simulate cycling of a graphite, MCMB (Meso Carbon Micro Beads)/ LiCoO 2 cell with formation of a film on the graphite electrode on charge. This simple model provided the capability to estimate the capacity fade for tens of thousands of low rate cycles in less than an hour of computation time. In 2006, Santhanagopalan et al. 3 presented a review of several lithium ion cell models including the extension of Atlung et al.'s model and named this model the single particle (SP) model. Zhang and White 18 used the SP model to make capacity fade predictions for small cells cycled at very low rate. Safari et al. 29 used the SP model to study capacity fade as did Deshpande et al. 30 who extended the SP model to include crack growth in a graphite particle followed by film formation in the resulting cracks with loss of lithium ions. Rahimian et al. 31 used the SP model to study several possible mechanisms that could cause capacity fade of lithium ion cells. Wang et al., 32 Santhanagopalan et al. 16 and Subramanian et al. 33 solved the single particle diffusion equations by polynomial approximation (2 nd order or higher). Prasad and Rahn 11 applied the SP model for identifying the aging parameters in lithium ion batteries. Guo and White 34 presented an approximate solution to the spherical diffusion problem with a fixed flux which provides high accuracy with only a few terms in a series when compared to the analytic solution with at least 100 terms. This approximate solution provides significant savings in computation time.
In order to predict the performance of lithium ion cells operating at higher rates, more complicated models were developed. In 1982, West et al. 4 presented a Nernst-Planck model for porous insertion electrodes for a Li/TiS 2 cell. They pointed out the importance of coupling the transport of lithium ions between the insertion electrode and the electrolyte. In 1993, Mao and White 35 presented a Nernst-Planck model of the Li/TiS 2 cell which included transport of lithium ions through the separator and demonstrated quantitatively how the utilization of the TiS 2 electrode could be improved by decreasing the thickness of the separator. Also in 1993, Doyle et al. 36 published a concentrated solution theory model for a lithium/polymer/TiS 2 cell. Shortly after that in 1994 Fuller et al. 5 published a concentrated solution theory model for a dual lithium ion insertion cell ("dual foil"), which is known as "pseudo-2D model (P2D)". In 1996, Doyle et al. 37 published a comparison of their model predictions to experimental data from a Bellcore plastic lithium ion cell. Unfortunately, the computation burden associated with using the program dual foil is too great to use their program for extensive simulations needed for parameter estimation and cell life predictions. Consequently, many papers have been published to provide methodology to reduce this computation burden.
To simplify the P2D models, assumptions have been made to decouple the nonlinear PDEs to reduce the numerical complexity. Doyle and Newman 6 derived analytic expressions for simplified models for two limiting cases: a constant reaction rate in the porous electrode, and no concentration gradient. Kumar, 9 Dao et al., 13 and Tanim et al. 15 assumed a constant reaction rate in each electrode assumption in their models. These results are of limited utility. Smith et al., 12 and Plett et al. 14 used Taylor series to linearize nonlinear Bulter-Volmer equation, which only applies when the overpotential is close to zero volt. Also, Smith et al. 12 (p2568) and Plett 14 (p218) assumed that it is reasonable to decouple the electrolyte potential from the electrolyte concentration by neglecting the concentration term of the electrolyte charge conservation equation. Models with these two assumptions lose the nonlinear features between the electrolyte potential and concentration, and introduce errors in the simulation results.
Another approach that has been used to reduce the computation burden for solving the P2D model equations is to use different mathematical techniques. Cai and White 38-40 used the proper orthogonal decomposition method and the orthogonal collocation method to solve the model equations for lithium ion cells to reduce the computation needed for simulation. Subramanian and his coworkers 7,8,10,41,42 applied coordinate transformation, orthogonal collocation and model reformulation techniques to facilitate the efficient simulation of the P2D model equations. In later work, Subramanian et al. 43 extended the polynomial approximation approach to the diffusion equation in electrolyte phase. Their coordinate transformation was aimed to rescale each domain to [0,1], which will reduce the problem from three regions to a single region and decrease the required computation. Similarly, the model reformulation technique was a series of mathematical operations of the governing equations, which will reduce the total numbers of governing equations without losing any accuracy. Han et al. 44,45 used several first order processes to approximate the diffusion equation and a quadratic function to simulate the electrolyte concentration. Dao et al. 13 used the Galerkin method with cosine functions and constant reaction rates to solve the P2D model equations. Howey et al. 46 used a spectral orthogonal collocation method to solve the P2D model equations. These mathematical approaches to approximate the solutions of the partial differential equations will introduce errors since a limited number of trial functions were used, such as orthogonal collocation, polynomial approximation, etc. These methods can save time and memory in simulating the steady input cell operations (i.e. constant current/voltage modes), but lose accuracy if the input signal varies rapidly with time (i.e., high-frequency current pulses). Smith et al. 12 developed a ROM by a control-oriented method: by linearizing the material thermodynamics and electrochemical kinetics, the P2D model is first converted to a single-input transfer function in the Laplace-domain, and then a state variable model (SVM) is derived using approximated poles and gains. Smith's SVM approach provides excellent convergence at high current rates (up to 50 C) and has been adopted into the multi-scale multi-dimensional framework 47,48 to model the distributed thermal behavior over large-format cell geometry. Using a similar approach to SVM, Plett and his coworkers 49,50 derived a ROM from the Laplace transfer function by using a discretetime realization algorithm (DRA) and the resulting model was used to simulate the Urban Dynamometer Driving Schedule (UDDS). Plett's approach achieves a significant speedup for the high-frequency simulation as compared with the nonlinear FOM. Following Smith's SVM method, Chen et al. 51 reported a reduced order model based on the uni-form concentration distribution in both solid and electrolyte phases. However, the SVM methods also have limits: to take the Laplace transfer, P2D model must be linearized around a specific SOC and the temperature must be assumed constant, therefore, their approach will lose accuracy when implemented over wide operation windows.
In this work, we present a discrete-time algorithm to solve the nonlinear P2D model equations. The partial differential equations for mass transfer are numerically translated into an ordinary-differentialequation (ODE) system which is solved using a 1 st order exponential integrator approach 52,53 at discrete time points. The charge conservation equations are included as algebraic constraints, and two procedures are presented to solve these nonlinear equations. It is shown that linearization of the Butler-Volmer equation can provide a good approximation for potentials to reach a fast convergence of the nonlinear Butler-Volmer equation. With our approach, the high-frequency drive cycles for Li-ion cells can be simulated through a 0%∼100% SOC window with 20 • C temperature rise, with both high accuracy and computation time efficiency. Our method has also been used to estimate model parameters from synthetic voltage data from the FOM for the US06 drive cycle by using a nonlinear least-square regression technique. A detailed description of our method is presented in the following sections.

Mathematical Model
Description of the P2D model.-The P2D modeling domains for a Li-ion cell are illustrated in Figure 1, where the linear dimension x goes through the thicknesses of electrodes and separator, the radius dimension r is defined between the center and the surface of spherical particles, R ex is the contact resistance between the cell and the external terminals, and V ex is the voltage drop across the external terminals. The governing equations and boundary conditions for the P2D model are listed in Table I, and the symbols for the variables are listed in Table II. The constant parameter values listed in Table III are cited from reference; 38 in our model, certain physical properties are assumed to vary with temperature through the Arrhenius expression: where ψ denotes a physical parameter, ψ| T =T ref is the value for that parameter at reference temperature, E * is the activation energy, and T ref = 298 K is the reference temperature. The Arrhenius coefficients for the temperature-dependent parameters are listed in Table IV.
State-variable model.-The P2D model (See Table I) using our nonlinear state variable method can be transformed into a nonlinear state-variable model (NSVM). The transport phenomena in the solid  Table I. Governing equations and boundary conditions for the pseudo-2D model.

Model Equations
Solution phase diffusion

Butler-Volmer equation for electrochemical reactions
Solid phase diffusion Not applied in the separator region (l n < x < l n + l s ) and electrolyte phases can be expressed by the following state variable equation:ẋ = Ax + Bd [2] where x is the state variable vector, A and B are coefficient matrices, and d is the source vector. These are explained in detail below. The charge conservation and electrochemical kinetic equations can be written as a non-linear algebraic constraint: where u is the input signal vector and y is the output signal. Unlike the work by Smith et al., 12 Lee et al. 48 and Leek et al. 49 in which the state equation was simplified into a single-input model; in our algorithm, the source term d is a multi-variable vector and is determined by the nonlinear constraint Equation 3. The state variable x includes the concentration variables for both solid and electrolyte phases, and the source term d contains the surface electrochemical current densities at the nodes in the discretized x coordinate. The output y is defined as the external voltage drop based on the applied current density (i app ) and cell temperature (T), Table II. List of variables for pseudo-2D model.

Variable name Symbol
Li + concentration in electrolyte (mol/m 3 ) c l Electric potential of electrolyte (V) φ l Electric potential of solid phase (V) φ s Surface electrochemical current density (A/m 2 ) Both the applied current density (i app ) and the cell temperature (T) can be measured in real-time, so these two variables are taken as the input vector u: In the following sections, we will introduce the transformation of the P2D model into the NSVM expressed by Equations 2 and 3.
Electrolyte diffusion.-The finite element 54 method was used to yield an ODE system from the electrolyte mass transfer equation (See Table I): where M, K, and F are, respectively, the mass, stiffness, and force matrices, c l is the vector containing the discretized electrolyte concentration values distributed through different nodes in the electrodes and separator domains, j s is the vector containing the discretized surface electrochemical current density values distributed through different meshed elements in the two electrode domains. The dimensions of matrices and vectors in Equation Form the eigen-decomposition matrix for M −1 K: where D = diag (λ 1 , λ 2 , · · ·) is a diagonal matrix formed by eigenvalues (λ i ) and V is the matrix whose column vectors are the eigenvectors of M −1 K. Multiply Equation 7 by V −1 to obtain: , and d = j s , so that the electrolyte mass transport equation is in the state equation form, The concentration vector can be obtained from x: where x is included in the algebraic constraint (3). It can be shown from the finite element formula that in the matrix D = diag (λ 1 , λ 2 , · · ·), each eigenvalue is proportional to the electrolyte diffusivity D l,bulk and is therefore scale by a factor of exp (where E * is the activation energy for D l,bulk ) during a temperature change; and as a result, extra computation for the eigen decomposition 8 is not required at each time step.
Solid phase diffusion.-To transform the solid particle diffusion PDE into a state equation, we used an approximate transfer function approach presented earlier, as shown in Appendix A, which is also used by Smith et al., 12 Lee et al. 48 and Leek et al. 49 The resulting transfer function is: Whereθ * ,θ avg , andδ s are respectively the dimensionless Laplace transform for the concentration at the surface of the particles, the average concentration in the spherical particle, and the surface flux for solid phase diffusion, and β is the dimensionless Laplace variable. We used a finite Heaviside expansion G (β) to approximate G(β): where a i and b i are adjustable parameters obtained by minimizing the absolute value |G(β) − G (β)| through frequency response after substituting β = ω j (where ω is the dimensionless angular frequency) in Equation 12. Therefore after parameter optimization, Equation A17 can be expanded into the following series: [13] whereQ i is the eigenfunction defined as: So the inverse Laplace transform L −1 {Q(β)} = Q(τ) forQ i (β) is equivalent to the solution of following ODE: Scale the dimensionless time τ to time t through expression A5 and use Equation A8 to substitute the dimensionless flux δ s , and Equation 15 can be converted into the dimensional time domain: The inverse Laplace transform for θ avg can be obtained from Equation A16, and then the derivative of θ avg to time could be derived as: a 2 , · · · ,a N ,0), and B = 1 F Rscs,max [b 1 ,b 2 , · · · ,b N , − 3] τ , and the solid phase diffusion equation is formulated to be in the same form as the state Equation 2: The dimensionless surface concentration θ * is expressed as θ * = [1, 1, · · · , 1] x and is included in constraint Equation 3. As the surface electrochemical current density j s (t, x) is expressed as the discretized vector j s through the electrode domains, the reformulated solid phase diffusion equations are repeatedly implemented for each element of j s .
Electrical submodels.-Initial solution.-Our formulation of the electrical submodel (the coupled charge conservation and Butler-Volmer equations) begins with the linearization of electrochemical kinetics. The nonlinear Butler-Volmer equation is expressed as follow: where, j 0 is calculated by: [19] with c 1 and c s are the concentrations of Li species in the electrolyte and solid phases respectively, k r is the reaction rate constant, and c s,max is the maximum concentration of Li + in the solid phase. Taking 1 st order Taylor expansion for the Butler-Volmer Equation 18, an approximate expression for the surface electrochemical current density j s can be obtained: The superscript indicates the variable is obtained by linearized BV function. Substitute Equation 20 into the charge conservation Equation 3, and the electrical submodel can be simplified as follows: i app is negative for discharge cycle and positive for charge cycle. Equations 21 through 24 can be solved linearly through the meshed geometry but the solutions are inaccurate due to the linear approximation, therefore we use the symbols φ l and φ s to distinguish them from full-order solutions. Figure 2 shows the general procedure for solving the electrical submodel: the linear solution from Equations 21 through 24 is used as initial values for a refinement subroutine from which variables j s and V ex = φ s | x=ln+ls+lp − φ s | x=0 + i app R ex can be evaluated with much higher accuracy. In this work, we developed two approaches to refine the solutions for electrical submodel: Method 1: Using a Newton loop shown in Figure 3a to solve the Butler-Volmer equation iteratively, and this method is called as "implicit method" Method 2: Approximate the Butler-Volmer equations into quadratic equations, and solve the equations explicitly by following the steps shown in Figure 3b; therefore this method is called as "explicit method" Implicit method.-The implicit method uses the iterative loop as shown in Figure 3a. At the 1 st iteration, the electrolyte and solid phase potential values are set equal to the linear solutions φ l and φ s : The initial value for j s is calculated using Equation 20, and the initial reference potential values are defined as follows: [26] where φ ref l and φ ref s are respectively the absolute references for electrolyte and solid phase potentials. The residue for Butler-Volmer equation, Res BV , is expressed as: where η = φ s − φ l − U denotes the overpotentials for electrochemical reactions. If the electrode domains are discretized into N e nodes, Equation 27 will be implemented at each node, and therefore Res BV is expressed by a N e × 1 column vector Res BV [29] where Res n and Res p are scalars. Therefore, the full residue vector for the electrical submodel is expressed as follows: where Res is a (N e + 2) × 1 column vector.  [31] where the Jacobian matrix Jac has (N e + 2) × (N e + 2) dimension. The step change for the unknown variables, dy, is evaluated as follow: Journal of The Electrochemical Society, 164 (11) E3001-E3015 (2017) and the unknown variables are updated for the next iteration through the following equation: With the updated j s , φ ref l , and φ ref s values, the electric potentials φ l and φ s can be solved from the charge conservation equations with the imposed boundary conditions: ) unless CC License in place (see abstract Explicit method.-As shown above, the implicit method involves the derivation of analytical Jacobian matrix, which is usually difficult for in-house code scripting. Therefore, as shown in Figure 3b, we developed a simplified approach to solve the electrical submodels explicitly without iterative loops. In the explicit approach, the overpotential is split into two parts: η =η + η ref [36] The first term on the right-hand-side of Equation 36 is expressed as: [37] and as shown in Equation 37,η varies with both time and space. The second term on the right-hand-side of Equation 36 is expressed as: [38] and in each electrode region, η ref is a variable only varying with time. Therefore, the Butler-Volmer equation can be written as follow: [39] Next, substitute Equation 39 into the electrical limiting Equation 28, to obtain the following hyperbolic equations for each electrode region: where the unknown variable X is expressed as: [41] and coefficients A and B are, respectively, expressed as: From Equation 40, the unknown variable X can be solved explicitly: [43] and η ref can be derived from Equation 41: According to Equation 38, the reference potentials can be calculated from η ref values in the two electrode regions: The terminal voltage V ex can be evaluated from Equation 35, and j s is calculated by the Butler-Volmer equation. According to the flowchart shown in Figure 3b, j s and V ex can be calculated approximately by only going through several sequential steps and no iteration is needed. Formulation of electrical equations.-As shown in Figure 2, the electrical submodel takes i app as a direct input, while the state variable x and input temperature T are also involved to evaluate certain coefficients (i.e. j 0 , U , κ eff , and κ d ). The outputs of the electrical submodel include j s (expressed as vector j s over the discretized nodes) and V ex ; as defined previously, d = j s and y = V ex , so the electrical submodel is formulated same as the constraint Equation 3: Simulation on discrete time domain.-For discrete time simulation, it is usually assumed that all coefficients and the source terms in state Equation 2 are constant during a small time interval [t k−1 ,t k ]. As shown in sections Electrolyte diffusion and Solid phase diffusion, the coefficient matrix A in Equation 2 is diagonal: [46] where λ i (i = 1, 2, 3, · · · N s ) are the diagonal elements evaluated at time point t k−1 . Before moving to the next time point t k = t k−1 + t, the following matrix exponential operators are defined: and the state variables at the next time point, x k , can be calculated through the following equation: [49] where x k−1 , d k−1 , and B are respectively the variable and coefficient values at t k−1 , and t is the time interval between t k−1 and t k . Note that if λ i = 0, the corresponding diagonal element in Equation 48 is estimated by taking the limit for λ i → 0: At time point t k , the source term d k and output y k can be obtained by solving the algebraic constraint described in Electrical submodels section: where u k is the input vector at t k . This simulation procedure is illustrated by a flowchart in Figure

Results and Discussion
Particle diffusion parameters.-As stated in Solid phase diffusion section, several dimensionless parameters (a i , b i , where i = 1, 2, · · · N ) in Equation 12 should be optimized by frequency response so that the transfer function for the particle diffusion can be approximated by a Heaviside series. To determine these parameter values, the dimensionless Laplace variable β is expressed into imaginary frequency: where ω is the dimensionless angular frequency and j is the unit imaginary number. Substitute Equation 52 into Equations 11 and 12, scan ω from 10 −6 through 10 6 , evaluate the dimensionless transfer functions G(β) and G (β) at each frequency; and the profile of G (β) can be fit to G(β) using a least square approach. The optimization results in Figure 5 show that a good agreement between G (β) and G(β) can be achieved with only 4 terms (N = 4), and the optimized a i and b i values are presented in Table V. These parameters are implemented in Equation 16 to simulate the solid phase diffusion in both anode and cathode.   We use this approximate transfer function approach only for the solid phase diffusion; in this case, the transfer function G(β) contains only one dimensionless group β, so the values for a i and b i are unaffected by other parameters or temperature. If this approach is extended to the full-cell scale, the resulting transfer function will take the form like G (s, κ eff , U, T ), so the optimized ROM parameters would implicitly depend on temperature, SOC, and electrolyte concentration. As a result, the ROM developed by a full-cell transfer function must be re-evaluated for different operating windows, however our NSVM method has no such limitations.
Drive cycle simulation.-The NSVM developed in this work was coded with MATLAB into two versions: the implicit model solves the algebraic constraint 51 using the method shown in Implicit method section, while the explicit model uses method shown in Explicit method. The cell domain is meshed into 12 quadratic finite elements (5 elements for each electrode and 2 elements for separator), and both NSVM versions includes 75 states (c l (25 nodes from 12 quadric elements) + c s (5 states (N = 4 + θ avg ) × 10 electrode elements)). Using the linear assumptions given by reference 12 (i.e. a fully-linearized Butler-Volmer equation and constant electrolyte conductivity), a statevariable model (SVM) was also developed with MATLAB to compare to our NSVM. A rigorous pseudo-2D model was developed as baseline using livelink between MATLAB and COMSOL5.2a, in which COMSOL is called by MATLAB as a numerical solver. As a result, the SVM, the NSVM, and the baseline model are compared under the same MATLAB environment.
Input profiles.-The input current signal was synthesized using the speed profile of US06 drive cycles, 55,56 the maximum discharge current rate was scaled to 2.5C (see Figure 6a). A repetition of eight drive cycles were simulated to cover the full 0∼100% SOC window. A temperature profile rising from 25 • C to 45 • C was also artificially synthesized with random noise (see Figure 6b).
Simulation results.-With the synthesized input profiles, the two versions of NSVM and the baseline model were implemented for the drive cycle simulation with full SOC window, the sampling time interval t was set at 1 sec, and a 2.75 V cutoff voltage was applied as the stop condition. The simulated drive cycle hit the stop voltage limit at t = 4471 sec, so there were totally 4472 solution sets computed during the simulation. The implicit NSVM requires 8.95 sec to simulate this drive cycle (2.0 msec per solution) and the explicit NSVM requires 8.50 sec (1.9 msec per solution). There is very little difference between the two NSVM versions in terms of simulation speed, so the iterative loop for the implicit model does not necessarily add to the computational load. However, the COMSOL model requires 470 sec (107 msec per sample) to process this simulation on the same PC, so the NSVM achieves a 50:1 speedup. Compared with NSVM, the SVM does not show much benefit in terms of computation time efficiency (8 sec for one simulation or 1.7 msec per solution).
The simulated terminal voltage responses (V ex ) are presented in Figure 7a, and the absolute error for the two NSVM versions and SVM are plotted in Figure 7b. According to these results, accuracy for the two different NSVM versions basically resides at the same level (absolute error < 15 mV or relative error< 0.5%); however, the SVM shows higher error (>30 mV) as compared with the NSVM. The error plots also show that the NSVM has relatively lower accuracy at the early and final states of discharge, and we attribute this to the steeper dU dθ slopes. The time-integral approach provided by Equations 47 through 50 requires all coefficients and the source term be constant over the small time interval [t k−1 ,t k ], a greater | dU dθ | value will increase variation to related parameters and therefore weaken the ground for this assumption. To check the simulation results for solid phase diffusion, the particle surface concentration θ was averaged through x dimension in each electrode(the anode-average θ * is defined as 1 ln ln 0 θ * dx and the cathode-average θ * is defined as 1 lp ln+ls+lp ln+ls θ * dx). The electrodeaverage θ vs time profiles are plotted in Figure 8, and both NSVM versions provided excellent agreement with the baseline results; therefore, the approximation approach presented in Solid phase diffusion section provides high accuracy for dynamic simulation. Simulation results for electrolyte diffusion are shown in Figure 9, the boundary concentration vs time plots (Figure 9a and 9b) and the distributed concentration profiles (Figure 9c) all confirm the good accuracy of NSVM. Results for electrolyte potential φ l are available in Figure 10, the reference electrolyte potential φ ref l = φ l | x=0 is plotted against time in Figure 10a, and two distributed φ l profiles at t = 332 sec and t = 3896 sec (the 2.4C high rate pulses occur at these time points) are presented in Figures 10b and 10c. According to Figure 10a, the NSVM basically fit with the baseline model and the maximum deviation is about 8 mV, however, this error may cause observable offsets to the potential distribution under a high current pulse at low SOC ( Figure  10c where SOC < 10%). In Figure 11, distributed profiles for the surface electrochemical current density j s are plotted at t = 332 sec and t = 3896 sec, and the deviation between NSVM and baseline model is higher near the electrode/separator interface because the gradients of electrolyte potential are largest at these interior boundaries.   Table IV are taken as the truth values for these parameters. The synthetic data for US06 drive cycle were created by running the baseline model with the truth parameter values. We initially scaled these parameters by a factor of 10, and then implemented NSVM for US06 simulation with the scaled parameter values; a MATLAB subroutine for nonlinear least-square regression was used to fit the NSVM to the synthetic data by iteratively adjusting the parameters. Comparisons between the simulated voltage profiles before and after the optimization are presented in Figure 12, and the optimized curves have much better fit to the synthetic data. Table VI lists the ratios of estimated parameters to truth values, and these numbers should ideally converge to 1; therefore the nonlinear parameter estimation in this case deviates from expectation by 2%∼17%. To check the cause of inaccuracy, we implemented NSVM with the true parameter values, and the simulated voltage results have an average error of 2.5 mV from the synthetic data; however, when the values listed in Table VI were implemented in NSVM, the average error dropped to 1.2 mV; therefore, these parameters were over-optimized.

Conclusions
The NSVM can be used to implement the pseudo-2D model for a Li-ion cell on a discrete time domain with multiple inputs. The formulated equations include the nonlinear features for pseudo-2D model, and the simulation speed of NSMV is approximately 50 times faster than that of a COMSOL baseline model. In general, the NSVM has excellent accuracy in predicting the concentrations of solid and electrolyte phases, but shows some error for the electrical potentials. The reason is due to the accumulation terms ∂c ∂t in mass transfer equations, which make the concentration results less sensitive to instantaneous current changes than the solutions of charge conservation equations. For the voltage response at cell terminals, which is determined by coupled mass transfer and charge conservation physics, the NSVM has ± 0.5% agreement with the rigorous baseline model. Two types of NSVM, the implicit and explicit versions, are presented; the implicit NSVM solves Butler-Vomer equations by an iterative Newton method, while the explicit NSVM approximates the nonlinear electrochemical kinetics by a quadratic equation. Both NSVM versions achieve approximately the same levels of accuracy and time efficiency, but the explicit NSVM reduces the complexity of code development. The NSVM can be used to characterize the parameter values in the pseudo-2D model by nonlinear regression.

Appendix A: Solid Phase Diffusion
The governing equation and boundary conditions for the solid phase diffusion are given as follows: Considering the temperature-dependency of diffusivities, D s can be a variable changing with time since the input temperature varies with time, and Laplace transform does not apply to the above equations; therefore, we need to first make the model dimensionless.
The modified variables are defined as follows: where θ andr are, respectively, the dimensionless concentration and radius coordinate. Using the expressions shown in A3, the governing equation can be converted as follow: The dimensionless mass flux at the surface of particle (r = R s ), δ s , is defined as follow: and the dimensionless boundary/initial conditions can be derived as: The Laplace transform can be defined on the dimensionless time domain: where f (τ) is a given function of τ and β is the dimensionless Laplace variable. Take the dimensionless Laplace transform for Equations A7 and A9 to obtain: The solution for the dimensionless concentrationθ(r , β) is: and the dimensionless concentration at the particle surfacer = 1 can be solved as: The dimensionless particle-average concentrationθ avg is defined as follows: andθ avg can be solved by multiply Equation A13 by 3r 2 and integrate from 0 to 1: Subtract equation A16 from A14 to obtain the following equation: where variableθ * (β) −θ avg (β) stands for the transient response of solid diffusion.