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

Meng Guo; Xinfang Jin; Ralph E. White

doi:10.1149/2.0021711jes

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).¹ 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² work in 1979 and named by Santhanagopalan et al.³ as the single particle (SP) model in 2006. Another most dominant Li-ion cell full-order model (FOM) was also extended from Atlung et al.'s⁴ work and developed in 1994 by Fuller et al.;⁵ 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)".³ 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–16 Also, ROMs have been published to simulate capacity fading^17,18 and for use in parameter estimations.^19–24

The single particle model dates back to 1979, when Atlung et al.² presented a relatively simple cell model for the Li/TiS₂ couple. Their model includes the diffusion of lithium ions in different shaped particles of TiS₂, 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.²⁵ to model the impedance of a spherical metal hydride particle. After that, Subramanian et al.²⁶ extended Atlung et al.'s model to include a variable diffusion coefficient. Also, Subramanian et al.²⁷ 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²⁸ used Atlung et al.'s model to simulate cycling of a graphite, MCMB (Meso Carbon Micro Beads)/ LiCoO₂ 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.³ 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¹⁸ used the SP model to make capacity fade predictions for small cells cycled at very low rate. Safari et al.²⁹ used the SP model to study capacity fade as did Deshpande et al.³⁰ 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.³¹ used the SP model to study several possible mechanisms that could cause capacity fade of lithium ion cells. Wang et al.,³² Santhanagopalan et al.¹⁶ and Subramanian et al.³³ solved the single particle diffusion equations by polynomial approximation (2^nd order or higher). Prasad and Rahn¹¹ applied the SP model for identifying the aging parameters in lithium ion batteries. Guo and White³⁴ 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.⁴ presented a Nernst-Planck model for porous insertion electrodes for a Li/TiS₂ 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³⁵ presented a Nernst-Planck model of the Li/TiS₂ cell which included transport of lithium ions through the separator and demonstrated quantitatively how the utilization of the TiS₂ electrode could be improved by decreasing the thickness of the separator. Also in 1993, Doyle et al.³⁶ published a concentrated solution theory model for a lithium/polymer/TiS₂ cell. Shortly after that in 1994 Fuller et al.⁵ 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.³⁷ 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⁶ derived analytic expressions for simplified models for two limiting cases: a constant reaction rate in the porous electrode, and no concentration gradient. Kumar,⁹ Dao et al.,¹³ and Tanim et al.¹⁵ assumed a constant reaction rate in each electrode assumption in their models. These results are of limited utility. Smith et al.,¹² and Plett et al.¹⁴ used Taylor series to linearize nonlinear Bulter-Volmer equation, which only applies when the overpotential is close to zero volt. Also, Smith et al.¹² (p2568) and Plett¹⁴ (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.⁴³ 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.¹³ used the Galerkin method with cosine functions and constant reaction rates to solve the P2D model equations. Howey et al.⁴⁶ 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.¹² 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 discrete-time 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.⁵¹ reported a reduced order model based on the uniform 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-differential-equation (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;³⁸ in our model, certain physical properties are assumed to vary with temperature through the Arrhenius expression:

where ψ denotes a physical parameter, 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.

**Figure 1.** Schematic of a Li-ion cell sandwich containing dual insertion porous electrodes separated by an ionically conducting membrane. The porous electrodes are formed by spherical insertion particles.

Table I. Governing equations and boundary conditions for the pseudo-2D model.

Model	Equations
Solution phase diffusion
Charge conservation in electrolyte
Charge conservation in solid phase

Butler-Volmer equation for electrochemical reactions
Solid phase diffusion

Table II. List of variables for pseudo-2D model.

Variable name	Symbol
Li⁺ concentration in electrolyte (mol/m³)	c_l
Electric potential of electrolyte (V)	ϕ_l
Electric potential of solid phase (V)	ϕ_s
Surface electrochemical current density (A/m²)	j_s
Exchange current density (A/m²)	j₀
Open circuit potential (V)	U
Electrochemical overpotential (V)	η
Li⁺ concentration in solid phase (mol/m³)	c_s

Table III. Constant parameter values for P2D model.³⁸

	Value
	Negative		Positive
Parameter	electrode	Separator	electrode
Maximum Li capacity c_{s, max} (mol/m³)	26390		22860
Radius of particle R_s (m)	12.5 × 10^{− 6}		8.0 × 10^{− 6}
Thickness of components l_nl_sl_p (m)	100 × 10^{− 6}	52 × 10^{− 6}	183 × 10^{− 6}
Porosity of electrode ε_l	0.357	1.0	0.444
Volume fraction of active material ε_s	0.471		0.297
Specific surface area a (1/m)	3ε_s/R_s		3ε_s/R_s
Solid electronic conductivity σ (S/m)	100		3.8
Universal gas constant R (J/mol/K)	8.314
Faraday constant F (C/mol)	96485
Reference temperature T^ref (K)	298
Electrolyte ion transfer number t⁺	0.363
1C rate current density i_1C (A/m²)	17.5
External connection resistance R_ex (Ωm²)	0.005
Average electrolyte concentration	2000
Activity coefficient of solution phase f_±	1

Table IV. Temperature-dependent properties for certain transport and kinetic parameters.

Parameter		Value at reference temperature T^ref	Activation energy E* [J/mol]	Units
Solid phase diffusivity D_s	Anode	3.9 × 10^{− 14}	2.0 × 10⁴	m²/s
	Cathode	1.0 × 10^{− 13}	2.0 × 10⁴
Reaction rate constant k_r	Anode	2.334 × 10^{− 11}	3.0 × 10⁴	m^2.5/mol^0.5/s
	Cathode	2.334 × 10^{− 11}	3.0 × 10⁴
Diffusivity in bulk electrolyte D_{l, bulk}		7.5 × 10^{− 11}	2.0 × 10⁴	m²/s
Electric conductivity in bulk electrolyte κ_bulk		See below^*	3.0 × 10⁴	S/m

^*κ_bulk = 4.153 × 10^{− 2} + 5.007 × 10^{− 4}c_l − 4.7212 × 10^{− 7}c²_l + 1.5094 × 10^{− 10}c³_l − 1.6018 × 10^{− 14}c⁴_l (where c_l in[mol/m³]) T^ref = 298K

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 and electrolyte phases can be expressed by the following state variable equation:

where is the state variable vector, and are coefficient matrices, and 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 is the input signal vector and y is the output signal. Unlike the work by Smith et al.,¹² Lee et al.⁴⁸ and Leek et al.⁴⁹ in which the state equation was simplified into a single-input model; in our algorithm, the source term is a multi-variable vector and is determined by the nonlinear constraint Equation 3. The state variable includes the concentration variables for both solid and electrolyte phases, and the source term 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),

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 :

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⁵⁴ method was used to yield an ODE system from the electrolyte mass transfer equation (See Table I):

where , , and are, respectively, the mass, stiffness, and force matrices, is the vector containing the discretized electrolyte concentration values distributed through different nodes in the electrodes and separator domains, 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 6 are determined specifically by the mesh pattern of the finite element method. Each element in and is a lumped variable depending only on time. Left-multiply Equation 6 by to obtain:

Form the eigen-decomposition matrix for :

where is a diagonal matrix formed by eigenvalues (λ'_i) and is the matrix whose column vectors are the eigenvectors of . Multiply Equation 7 by to obtain:

Let , , , and , so that the electrolyte mass transport equation is in the state equation form,

The concentration vector can be obtained from :

where is included in the algebraic constraint(3). It can be shown from the finite element formula that in the matrix , each eigenvalue is proportional to the electrolyte diffusivity D_{l, bulk} and is therefore scale by a factor of (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.,¹² Lee et al.⁴⁸ and Leek et al.⁴⁹

The resulting transfer function is:

Where , , and 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 (where is the dimensionless angular frequency) in Equation 12. Therefore after parameter optimization, Equation A17 can be expanded into the following series:

where is the eigenfunction defined as:

So the inverse Laplace transform for 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:

Let , , , and , 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 and is included in constraint Equation 3. As the surface electrochemical current density j_s(t, x) is expressed as the discretized vector through the electrode domains, the reformulated solid phase diffusion equations are repeatedly implemented for each element of .

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₀ is calculated by:

with c₁ 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 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"

**Figure 2.** General solution flowchart for electrical submodel which includes the charge conservation equations and the Butler-Volmer equation (See Fig. 3 for details of Refinement block).

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:

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 . The two reference potentials, ϕ^ref_l and ϕ^ref_s, are extra unknown variables which are determined by the electrical limiting equations defined as follows:

and the residues for Equation 28 are expressed as:

where Res_n and Res_p are scalars. Therefore, the full residue vector for the electrical submodel is expressed as follows:

where is a (N_e + 2) × 1 column vector. The equation system includes totally N_e + 2 unknowns: the discretized j_s values at N_e elements (j_s is expressed by a N_e × 1 column vector ) and the two reference potentials, ϕ^ref_l, and ϕ^ref_s. Next, we use Newton's method to decide the j_s, ϕ^ref_l, and ϕ^ref_s values that make the norm for residue vector . Derive the Jacobian matrix for the residue vector as:

where the Jacobian matrix has (N_e + 2) × (N_e + 2) dimension. The step change for the unknown variables, , is evaluated as follow:

and the unknown variables are updated for the next iteration through the following equation:

With the updated , ϕ^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:

With the updated solution for ϕ_l and ϕ_s, repeat steps described by Equation 26 through 34 until value is small enough (<10⁻⁸), then exit the loop and output the and ϕ^ref_s values from the last iteration. The terminal voltage V_ex is evaluated through the following equation:

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:

The first term on the right-hand-side of Equation 36 is expressed as:

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:

and in each electrode region, η^ref is a variable only varying with time. Therefore, the Butler-Volmer equation can be written as follow:

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:

and coefficients A and B are, respectively, expressed as:

From Equation 40, the unknown variable X can be solved explicitly:

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 and input temperature T are also involved to evaluate certain coefficients (i.e. j₀, U, κ_eff, and κ_d). The outputs of the electrical submodel include j_s (expressed as vector over the discretized nodes) and V_ex; as defined previously, 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 in Equation 2 is diagonal:

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, , can be calculated through the following equation:

where , , and 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 and output y_k can be obtained by solving the algebraic constraint described in Electrical submodels section:

where is the input vector at t_k.

This simulation procedure is illustrated by a flowchart in Figure 4 and can be described as follow:

(1)
At the initial time point (k = 0), the state vector equals to the preset initial values (in this work, the initial values are c_l = 2000 mol/m³, θ_n = 0.5635, and θ_p = 0.1706);
(2)
Evaluate the coefficients and operators in the model with and ;
(3)
Solve Equation 51 to obtain and y_k
(4)
Update the index k = k + 1;
(5)
Calculate from Equation 50 using and ;
(6)
Repeat steps 2) through 5).

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⁶, 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.

**Figure 5.** Comparisons between G(β) and G'(β) for dimensionless frequency ω sweeping from 10⁻⁶ through 10⁶ (plots only show values for ω > 10^{− 1} to have better resolution).

Table V. Optimized parameter values in approximated transfer function G'(β) with N = 4.

	i = 1	i = 2	i = 3	i = 4
a_i	35058.7	1382.966	141.595	22.32279
b_i	−268.261	−30.9242	−7.59606	−2.59525

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¹² (i.e. a fully-linearized Butler-Volmer equation and constant electrolyte conductivity), a state-variable 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).

**Figure 6.** Synthesized current and temperature profiles for practical drive cycles: (a) current profile scaled to C rate, (b) temperature profile.

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 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 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 and the cathode-average θ* is defined as ). The electrode-average θ 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.

**Figure 8.** Electrode-averaged θ vs time for (a) anode and (b) cathode.

**Figure 9.** Plots for electrolyte concentrations: (a) concentration vs time at boundary x = 0, (b) concentration vs time at boundary x = l_n + l_s + l_p, (c) concentration vs x at t = 655 s and t = 2845 s.

**Figure 10.** Plots for electrical potentials of electrolyte phase: (a) electrolyte potential vs time at boundary x = 0, (b) electrolyte potential vs x at t = 332 s, (c) electrolyte potential vs x at t = 3896 s.

**Figure 11.** Plots for surface electrochemical current density vs x: (a) at t = 332 s, (b) at t = 3896 s.

Parameter identification

As shown in Simulation results section, both versions of NSVM show excellent accuracy; also, the NSVM is much faster than the rigorous baseline model. Another advantage for NSVM is that the model formulation is unaffected by parameter values. Therefore, it can be used to estimate key parameter values in the pseudo-2D model. In this work, the solid phase diffusivities and reaction rate constants at reference temperature ( and ) are chosen as adjustable parameters, and values listed in 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.

**Figure 12.** The voltage responses simulated by NSVM compared with baseline results: (a) before parameter optimization, (b) after parameter optimization.

Table VI. Optimized key parameter values in NSVM through least-squares fitting.


	Anode	Cathode	Anode	Cathode
Implicit NSVM	0.857838	1.020308	1.022884	1.095047
Explicit NSVM	0.833355	1.102079	1.023376	1.103187

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 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.

List of Symbols

a	Active surface area, (m²/m³)
a_i	Approximated poles for the solid phase diffusion transfer function 12
	Diagonal stiffness matrix for the general state Equation 2
A	Coefficient involved to derive the explicit solution 42
	Forcing matrix for the general state Equation 2
B	Coefficient used to derive the explicit solution 42
b_i	Approximated gains for the solid phase diffusion transfer function 12
c₁	Lithium ion concentration in electrolyte, (mol/m³)
c_s	Lithium ion concentration in the solid phase, (mol/m³)
c_{s, ini}	Initial Lithium ion concentration in the solid phase, (mol/m³)
c_{s, max}	Maximum lithium ion concentration in the solid phase, (mol/m³)
	Vector form time-dependent electrolyte concentration, (mol/m³)
	Source vector in the general state Equation 2
	Diagonal stiffness matrix in formulated electrolyte diffusion Equation 9, (s⁻¹)
D_{l, bulk}	Diffusivity of Li⁺ in bulk electrolyte, (m²/s)
D_s	Lithium ion diffusivity in the solid phase, (m²/s)
E*	Activation Energy for transport and kinetic parameters, (J/mol)
	Forcing matrix in the finite element Equation 6 for electrolyte diffusion, (mol/C)
F	Faraday constant, 96485 (C/mol)
f_±	Mean activity coefficient of the electrolyte
G(β)	Full-order dimensionless transfer function
G'(β)	Finite Heaviside expansion of G(β)
	Vector form of surface electrochemical current density in the electrode regions, (A/m²)
j_s	Surface electrochemical current density in the electrode regions, (A/m²)
j₀	Exchange current density, (A/m²)
j	Unit imaginary number
k_r	Reaction rate constant, (m^2.5/mol^0.5/s)
i_app	Applied current density, (A/m²)
	Stiffness matrix in the finite element Equation 6 for electrolyte diffusion, (m/s)
l_n	Thickness of negative electrode, (m)
l_s	Thickness of seperator, (m)
l_p	Thickness of positive electrode, (m)
	Mass matrix in the finite element Equation 6 for electrolyte diffusion, (m)
N_e	Total number of meshed elements in the two electrode domains, N_e = 5 + 5 = 10
N	Total number of approximated Heaviside expansion terms for solid phase diffusion, N = 4
N_s	Total number of states, N_s = 75
	Laplace transform of eigenfunctions for solid phase diffusion
	Time derivative of inverse Laplace transform of eigenfunction (1/s)
Q_i	Inverse Laplace transform of eigenfunction for solid phase diffusion
R	Gas constant, 8.314 (J/mol/K)
R_s	Radius of the solid particles, (m)
Res_BV	Residue for Butler-Volmer equation, (A/m²)
	Residue vector for spatially discretized Butler-Volmer equation, (A/m²)
r	Radial coordinate of the particle in the electrodes, (m)
	Dimensionless radius coordinate
T	Temperature, (K)
t	Time, (s)
t⁺	Transference number of lithium ion species
	Input signal for NSVM
U	Open circuit potential, (V), listed in Table I
V_ex	External potential drop across the cell, (V)
	Matrix for eigenvectors
x	x coordinate, (m)
	State variable vector in Eq. 2
y	Output signal
	Jacobian matrix
	Step change for the unknown variables
R_ex	External resistance, (Ωm²)
X	Unknown variable, defined by Eq. 41

Greek

ψ	Temperature-dependent transport/kinetic parameters (Equation 1, Table IV)
θ	Dimensionless concentration of solid phase
θ_ini	Initial dimensionless concentration of solid phase
θ*	Dimensionless concentration at the surface of solid particles
	Laplace transform of dimensionless concentration of solid phase
	Dimensionless Laplace transform for the surface concentration
θ_avg	Dimensionless particle-average concentration of solid phase
	Dimensionless Laplace transform for the average concentration
	Dimensionless Laplace transform for the boundary mass flux for solid phase diffusion
β	dimensionless Laplace variable
τ	Dimensionless time
η	Electrochemical overpotential, (V)
	Electrochemical overpotential defined in Eq. 37 (V)
η^ref	Reference electrochemical overpotential defined in Eq. 38, (V)
ϕ_s	Electric Potential in the solid phase, (V)
ϕ_l	Electric Potential in the electrolyte phase, (V)
ϕ^ref_s	Absolute reference for solid phase potential, (V)
ϕ^ref_l	Absolute reference for electrolyte potential, (V)
ϕ'_s	Linear solution of electric Potential in the solid phase, (V)
ϕ'_l	Linear solution of electric Potential in the electrolyte phase, (V)
κ_bulk	Electric conductivity of bulk electrolyte, (S/m)
κ_eff	Effective conductivity of the electrolyte phase, (S/m)
κ_d	Conductivity resulted from diffusion in the electrolyte phase, (S/m)
σ_eff	Conductivity of the solid phase, (S/m)
	Matrix exponential operator defined in Eq. 47
	Matrix exponential operator defined in Eq. 48
	Dimensionless angular frequency
ε_l	Porosity of the void phase in the separator
ε_s	Porosity of the solid phase in the electrode
λ_i	ith diagonal element in Eq. 46
λ'_i	ith eigenvalues of matrix D for electrolyte diffusion equation in Eq. 8

Subscripts

i	Term i
l	Liquid phase (Electrolyte phase)
s	Solid phase
n	Negative electrode
p	Positive electrode
k	Time step k
eff	Effective
bulk	Bulk property

Superscripts

ref	Reference state

: 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 θ and are, respectively, the dimensionless concentration and radius coordinate. Using the expressions shown in A3, the governing equation can be converted as follow:

Next, the dimensionless time, τ, is defined by the integral of with time:

and the left-hand-side of equation A4 can be derived as follow:

Substitute Equation A6 into A4 and the time-dependent variable term can be canceled form both sides to yield:

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: