Nonlinear State-Variable Method (NSVM) for Li-Ion Batteries: Finite-Element Method and Control Mode

The ﬁnite element method (FEM) was used in our nonlinear state-variable method (NSVM) presented recently ( J. Electrochem. Soc. , 164 , E3001 (2017)). The details of the application of the FEM to solve the lithium ion pseudo-2D (P2D) model equations using the NSVM are presented here for several control modes (constant current, voltage, power, or load). Validation of the method was performed by comparison to rigorous full-order models and experimental data. The FEM based NSVM shows excellent performance, and the estimated cell parameters are determined with a high conﬁdence level. © The Author(s) 2017

The pseudo-2D (P2D) model for Li-ion cells 1 has been widely applied. This model includes two spatial dimensions and several partial differential equations (PDE) which are coupled by nonlinear electrochemical kinetics. 2 The basic solution approach for the P2D model is to discretize the spatially-distributed PDE into a differential-algebraic equation (DAE) system. The finite difference method (FDM), finite volume method (FVM), and finite element method (FEM) have been used to discretize the P2D model. [3][4][5][6][7][8][9][10][11] For physics-based simulation of Li-ion batteries, the FEM is an excellent solution approach especially for complex multi-scale cell domains. [12][13][14] Although there are commercial finite element software packages (e.g., COMSOL Multiphysics) with the capability of solving the model equations, in-house code developed with scripting languages (i.e., MATLAB, FORTRAN, and Python) are still important due to their low cost and flexibility.
In our previous work, 15 we presented a nonlinear state variable modeling (NSVM) algorithm to reduce the computational cost for solving the P2D model over a high-frequency time domain. According to that NSVM work, the state variables in the diffusion equations are solved using the following approach: x k = x k−1 + Bd k−1 [1] where x k and x k−1 are vectors for the concentration states at current and previous time instances, d k−1 is the source term vector at the previous time instance, and , , and B are matrix operators. The source terms are solved by a nonlinear algebraic constraint: where u k and y k are respectively the inputs and outputs of the model. The above mentioned approach enables the P2D model to be solved under a high-frequency current signal with significant time efficiency (50+ times faster than a COMSOL baseline model). However, Reference 15 does not show the full details of the NSVM due to the limits on article length; in Equation 1, operators , , and B come from a specific FEM formulation which was not presented in that article; in Equation 2, the detailed solution procedure is also based on the FEM formulation. Also, the published NSVM algorithm can only simulate current control and the results have not been compared with real data. Furthermore, some derivation steps involved in the FEM solution procedure can be applied potentially to develop a reduced-order-modeling (ROM) approach. Consequently, the following is presented in this article: 1) our procedure to implement the finite-element method for the P2D model equations; 2) extension of the NSVM to apply for voltage, power, and load control modes; 3) validation of the NSVM by comparison to experimental data.

Mathematical Model
As illustrated by Figure 1, the P2D model was developed based on a dual-foil representation for a Li-ion cell, in which the x dimension runs through the thicknesses of the domain layers (anode, separator, and cathode) and the r dimension represents the radius of the spherical active material particles. Several transport equations are defined, respectively, over the x and r dimensions and are coupled by the surface electrochemical reaction kinetics. These physics-based equations for the P2D model are listed in Table I. The solution approach for the solid particle diffusion equation (r dimension) was discussed previously. 15 Consequently, the following sections will focus on the equations in the x dimension.
Formulation of the governing equations.-Three partial differential equations (PDEs) in the x dimension in the electrode domains (diffusion and charge conservation in the electrolyte and charge conservation in the solid phase) are solved using the finite element method, and all the PDEs can be written in the following standard format: [3] where t is time, w(x, t) is a field variable, ρ(x) is the mass function, μ(x, t) is the transport coefficient, γ(x) is the forcing function, and j s (x, t) is the current density for surface electrochemical reactions. The governing equation for the electrolyte charge conservation is given by: where φ L is the electrical potential of the electrolyte, κ eff L is the effective electrical conductivity of the electrolyte, c L is the concentration of the electrolyte, t + is the electrolyte transference number, a is the specific surface area of the solid phase, T is temperature, R is the universal gas constant, and F is the Faraday constant. To make Equation 4 consistent with the standard form, the diffusional current density where φ ref L = φ L | x=0 is the reference electrolyte potential selected for the boundary x = 0. Therefore, Equation 4 can be simplified to contain only one field variable φ * L : and it can be derived from Equation 5 that φ * L = 0 at boundary x = 0: where Equation 7 serves as the imposed boundary condition to determinate the absolute electrolyte potential values. The charge conservation in the solid phase is described by the following equation: where φ S is the potential of the solid phase and σ eff S is the effective electrical conductivity of the solid phase. As the solid phase potentials are distributed through two separate regions (anode and cathode), two imposed boundary conditions are set for both the anode at x = 0 and the cathode at x = l n + l s + l p : [9] Table I. Equations for pseudo-2D (P2D) model of Li-ion cell.

Physics
Equations (Defined only for 0 ≤ x ≤ l n and l n + l s ≤ x ≤ l n + l s + l p ) c S,max = 0 (Defined only for 0 ≤ x ≤ l n and l n + l s ≤ x ≤ l n + l s + l p ) Butler-Volmer equation RT η (Defined only for 0 ≤ x ≤ l n and l n + l s ≤ x ≤ l n + l s + l p ), R film = 0 for cathode Journal of The Electrochemical Society, 164 (11) E3200-E3214 (2017) where φ ref S is the reference electrical potential of the cathode. Introduce the modified solid phase potential as follow: For cathode domain l n + l s ≤ x ≤ l n + l s + l p [10] The governing equation for φ * S is same as Equation 8, and the imposed boundary conditions for φ * S are made homogeneous: Accordingly, the overpotential of the electrodes can be expressed in terms of the modified potentials: [13] where U * denotes the modified open circuit potential expressed by [14] Note that the SEI resistance R film is neglected for the cathode and the definition of θ * is given in  [15] The boundary current density for the solid phase is defined as follows: where i B is the current density passing through the interfaces between the electrodes and the current collectors. Substituting boundary conditions 16 into the right-hand-side of Equation 15, the first limiting electrical equation is found to be: The diffusion of Li ions in the electrolyte phase is described by the following equation: where c L is the concentration of electrolyte, ε L is the volume fraction of electrolyte in the porous electrode, and D eff L is the effective diffusivity of the electrolyte in the porous electrode. As shown above, Equations 6, 11, and 19 are all consistent with the standard form given in Equation 3, and the interpretation for these equation variables are listed in Table II.
Mesh pattern for modeling domains.-The pseudo-2D model includes three computational domains: anode, separator, and cathode, and these domains are discretized into segmental elements ( 1 , 2 · · · shown in Figure 1). Let m 1 , m 2 , and m 3 , respectively, denote the numbers of elements in the anode, separator, and cathode, so that Figure  1 presents a mesh pattern with m 1 = 5, m 2 = 2, and m 3 = 5. Three node points are defined, respectively, at the two ends and the center of each element. Figure 2 illustrates the indexing of nodes and elements, where the three nodes of the i th element ( i ) are, respectively, labeled as x 2i−1 ,x 2i , and x 2i+1 , and node x 2i+1 is shared by two adjacent elements i and i+1 as the interior boundary. There are in total 2 × (m 1 + m 2 + m 3 ) + 1 nodes in the x dimension, where each domain contains 2m k=1,2,3 + 1 nodes. Several integer parameters are specifically defined to facilitate our derivations: where m denotes the total number of elements through the entire geometry, m denotes the number of elements included in the two electrode domains where electrochemical reactions occur, n 1 is the node index for anode/separator interface, n 2 is the node index for separator/cathode interface, n is the number of nodes included in the two electrode domains, and n is the total number of nodes through the entire geometry.
Basis functions.-Two types of basis functions are developed over the discretized geometry as shown in Figure 3. A series of piecewise quadratic polynomials p 1 (x), p 2 (x), · · · , p n (x), whose expressions are given as below, are defined through node points: In other elements The profiles for Equations 21a through 21d are respectively illustrated by Figures 3a, 3b, 3c, and 3d; it can be shown that the value of p i (x) at node point x j (i, j = 1, 2, · · · , n) makes a Kronecker delta function: A series of boxcar functions q 1 (x), q 2 (x), · · · , q m (x) are defined over the elements, whose expressions are given by: and Equation 23 is illustrated by Figure 3e. Table III lists the variable vectors in the discretized P2D model. According to Equations (R-1) through (R-3), vectors c L (t), L (t), and S (t) are made up by the nodal values of the field variables c L (x, t), φ * L (x, t), and φ * S (x, t); the volume-average for these field variables in each element,ĉ L,i (t),φ * L,i (t), andφ * S,i (t), are defined by Equations (R-4) through (R-6), and vectorsĉ L (t),ˆ L (t), andˆ S (t) are formed by these element average values as shown in Equations (R-7) through(R-9). The electrochemical current density j s (x, t) is assumed to be uniformly distributed in each element, and variables j s,i (t), which form vector j s (t) as shown in Equation (R-10), stands for the homogenized element value for j s (x, t). The vectors for basis functions are given by Equations (R-11) though (R-14), where vectors p asc (x) and q asc (x) are respectively formed by the node and element basis functions throughout the entire geometry(anode, separator, and cathode), and vectors p ac (x) and q ac (x) are respectively formed by the node and element basis functions in the two electrode domains. The effective electrolyte conductivity κ eff L (x, t) is evaluated

Discretized variables.-
where I stands for identity matrix and 0 stands for all -zeroes matrix (R-21) usingĉ L,i (t) in each element according to Equation (R-15), where the expression for the bulk electrolyte conductivity κ L,bulk (c L , T ) is given in Appendix A; and κ eff L (x, t) can be calculated by following expression: where κ eff L,i is the effective electrolyte conductivity for element i as defined by (R-15). In each element, the solid phase diffusion can be regarded as a submodel in which j s,i (t) serves as the single input, and θ * i denotes the surface state-of-charge value corresponding to j s,i (t); the solution procedure for the solid phase diffusion is shown in Appendix B. Kinetic variables related to the electrochemical reactions are also homogenized in each element, in Equations (R-16) through (R-18), variables U * i , η i , and j ex,i respectively denote the modified open circuit potential, the overpotential, and the exchange current density values in each element; and in Equations (R-19) and (R-20), variables U * i and η i form vectors U and η. With Equations (R-21) and (R-22), the overpotential vector η can be expressed as follow: The Butler-Volmer equation in element i is given by: where i = 1, 2, · · · , m 1 , m 1 + m 2 + 1, m 1 + m 2 + 2, · · · , m. Take 1 st order Taylor expansion for the right-hand-side of Equation 26: and the following approximation can be derived from Equations 26 and 27: where the expression for the kinetic resistance Z k,i is given by (R-23). Equation 28 can be expanded in matrix-vector format through a diagonal matrix Z k defined in Equation (R-24): ) unless CC License in place (see abstract   Equation 3, the field variable w(x, t) and source term j s (x, t) are, respectively, approximated by

Approximation for the dependent variables.-In
where expressions for the vectors p(x) and w(t) are given in Table  IV. The element-average values for w(x, t) are expressed byŵ(t) according to Table IV, and there exist a linear transform from w(t) tô w(t) given by:ŵ (t) = Hw (t) [32] in which the constant operator H is defined in Appendix C.
Weak form of the transport equations.-In this work, p(x) is used as test function to derive the weak form; left multiply Equation 3 by the test function p(x) and integrate through the domains to obtain, where the definition for operator dx is given in Table IV. On the right hand side of Equation 33, the first term can be expanded using integration by parts: where n B denotes the direction pointing out of the boundary, and Bi is the index set for the non-insulation boundary nodes (see Table IV).
Therefore Equation 33 can be written as follows: Substitute w(x, t) and j s (x, t) with the approximate expressions shown in Equations 30 and 31, and Equation 35 can be written as follows: where M is the mass matrix defined by K is the stiffness matrix defined by F is the forcing matrix defined by and L B is the boundary loading vector defined by Starting from the weak form Equation 36, we have presented in Ref. 15 the solution procedure for the electrolyte diffusion (in which L B = 0 for the insulated boundaries) over a discrete time domain. In the following sections of this work, we will focus on solving the coupled charge conservation and electrochemical kinetics under different control modes.  Table IV; therefore, the weak form Equation 36 can be written as follows: The reference potential Equations 7 and 9 are expressed as: where the definition for the reference boundary operator Br is given in Table IV. The limiting electrical Equations 17 and 18 are, respectively, expressed as follows: and the operators C 1 and C 2 are given by: and the expression for w(t) can be written as follows: According to Equations 32 and 48, the element-average potentials can be expressed as follows: where the transport resistance operator is given by: Equation 49 can be written, respectively, for the electrolyte and solid phases: Substitute expressions 51 into Equation 25, so that the overpotential can be written as a linear expression for j [52] Model linearization and approximate solution.-To solve the nonlinear Butler-Volmer equations with the Newton method, good initial guesses are important to reduce the numerical difficulties. Initial values can be determined from the linearized kinetics. Substitute the approximate expression 29 into Equation 52 to yield: where the operator Z is given by: Substitute expression 54 into the limiting electrical Equations 43 and 44 to yield Left multiply Equation 57 by [0, 1] to yield the expression for φ ref where the coefficients u B and z B are given by: According to Equations 54 and 57, values forj s (t),φ ref L (t), and φ ref S (t) can all be determined by setting i B (t).
Electrical control equations.-For the current control mode, i B (t) can be evaluated using the following equation: where I app (t) is the total current applied to the cell which is defined as positive in charge and negative in discharge, and A C is the total area of the cathode/current collector interface. For voltage control mode, the control equation is given as follows: where V app (t) is the applied voltage and the terminal current I ex and terminal voltage V ex are defined as follows: Substituting φ ref S (t) with expression 58 and according to Equations 62 and 63, the voltage control equation can be expressed as follow: [64] and i B (t) can be solved from Equation 64 as follow: For the power control mode, the control equation is given as follows: where P app (t) is the applied power. According to Equations 58, 62, and 63, the control Equation 66 can be expressed as follows: and i B can be determined from Equation 67 through the quadratic solution: Equation 68 contains a quadratic discriminant D m which is determined as follows: and the solution for i B according to Equation 68 is valid only for D m ≥ 0. Under high power conditions, the value for D m might be negative, which means the linearized Butler-Volmer kinetics cannot guarantee a valid solution; in this case try the following equation: For the load control mode, the correlation between terminal voltage V ex and the load resistance R load (t) is as follows: Note that the i B (t) values evaluated above are not accurate, they are only used to estimate the approximate solutionsj s (t),φ ref L (t), and φ ref S (t), which are used as initial values for the nonlinear step.
Nonlinear refinement of solutions.-According to Equation 26, the residue for Butler-Volmer equation in element i is expressed as follows: [73] where R BV,i is the Butler-Volmer equation residue and i = 1, 2, · · · , m 1 , m 1 + m 2 + 1, m 1 + m 2 + 2, · · · , m. The residue vector for Butler-Volmer equations R BV are given by: According to Equations 43 and 44, the residues for the electrical limiting equations, R LE , are given by: For the current control mode, the residue for the control equation, R CE , is given by: for voltage control mode, R CE given by: and for power control mode, R CE is given by: Define the unknown vector y and the model residue vector Res as follows: The Jacobian matrix Jac is expressed as the derivative of Res with respect to y: and the step change for the unknown vector y in each iteration can be calculated by: The solution for the nonlinear equations can be estimated with high accuracy using the Newton loop shown in Figure 4, where y (k) denotes the value of y at the k th iteration, and y is updated with Equation 81 until the absolute value |Res| is smaller than the tolerance δ (i.e., δ = 10 −8 ).

Results and Discussion
In the previous work, 15 we have examined the accuracy of our NSVM algorithm for current control mode through model-to-model comparisons. In this work, we focus on model vs experimental validation and simulation of different control modes.
Model validation.-The computational algorithm described in Mathematical Model section was applied in modeling the Samsung commercial cells (the nominal cell capacity is 2.6 Ah), where the cell design parameters and basic physical properties are given in Appendix A. Three cells, labeled as LIB4, LIB6, and LIB7, respectively, were tested at room temperature (T = 25 o C) using different protocols. The experimental details for these cells are listed in Table V. MATLAB was used to develop the code for our cell model.
In simulation of the mixed cell control modes, the voltage control algorithm was applied for the constant voltage steps and the current control algorithm was applied for other steps. The cell temperature change was neglected in these simulations, and the SEI film resistance was only applied for the anode. Using the least-squares method,   several cell parameters were estimated to fit the model to the test data. Among the estimated cell parameters, the solid phase diffusivities (D s ) and exchange current densities ( j ex ) for the anode and cathode materials were assumed to be the same for all three cells, while the anode SEI film resistances varied between the cells. The optimized model results with comparisons to test data are presented in Figure  5 and the estimated parameter values are listed in Table VI. These results show that our NSVM algorithm works well with current and voltage control modes and the error for the model is below 2% for most parts of the simulated voltage profiles.
Power and load control simulation.-Using the optimized parameter sets from the previous section, our NSVM algorithm was implemented to simulate a power control step. The high frequency power input was scaled from the drive cycle current profile given in Ref. 15, where the applied power P app is calculated as follow: [82] where I 1C = 2.6A is the 1 C rate current for the cell. The input profiles are shown in Figure 6, the power signal plot includes 8 periodic cycles and the peak power pulse of each cycle is around 17 W. The sampling frequency for the input power is 1 Hz or 1 data per second. The power control simulation was implemented for cell LIB4, and a rigorous model developed by COMSOL 5.2 was used as the baseline, in which the maximum time step is set to be 1 sec in accordance with the input sampling rate. The simulated voltage and current profiles are presented in Figure 7, and our algorithm shows excellent accuracy as compared with the baseline mode. However, it takes only 3.06 sec for our algorithm to finish the simulation as compared to the 762 sec for the baseline model. The simulated current profile in Figure 7a shows   that the magnitude for peak current in each cycle increases toward the end of discharge; the reason is that the power control mode requires I ex V ex = P app , as the cell voltage V ex drops with the depth of discharge, more current is needed to meet the power request. Other important results are presented in Figure 8 for analysis. The profile for the quadratic determinant D m , which is defined in Equation 69, is plotted in Figure 8a and the results show that D m takes negative values at high power pulses near the end of discharge; this suggests that a fully linearized model fails to provide valid solution at these points. In Figure 8b, the current values estimated by Equation 68 or Equation 70 are compared to those calculated through the nonlinear refinement loop (see Figure 4), the plots show that there are signif-icant deviations between the two simulated current curves near the end of discharge, and these results suggest that the refinement loop is necessary to achieve good accuracy. To simulate the load control mode, a dynamic profile for the resistance load R load (t) was synthesized and presented in Figure 9. The simulated current and voltage profiles for the load control discharge are presented in Figures 10a  and 10b. Similar to the power control results, the NSVM algorithm has excellent agreement with the baseline model, and the simulation takes only 3.24 sec for NSVM algorithm as contrast to the 510 sec for the baseline model. The error plots shown in Figure 11 confirm the accuracy of NSVM approach in the power and load control cases.  Numerical discussion.-As shown in Power and load control simulation section, the NSVM runs much faster than COMSOL and provides excellent accuracy. The time step approach presented by Equation 1 is the main reason for this significant improvement in simulation time efficiency. However, the FEM approach described in this article also contributes to this speedup. In COMSOL, the electrolyte and solid phase potentials (φ L (x, t) and φ S (x, t)) are solved iteratively with the electrochemical current density ( j s (x, t)) at each node; in NSVM, only j s (x, t) is solved in the Newton loop and potential values are derived as linear expressions of j s (x, t). In addition, COMSOL defines the discretized current density j s (x, t) values at each node while NSVM defines j s (x, t) in each element, this means vector j s (t) has smaller size in NSVM than in COMSOL. As j s (t) is coupled with the solid phase diffusion, NSVM also includes fewer solid phase concentration variables than COMSOL. The FEM approach in NSVM greatly simplifies the solution procedure for the constraint Equation 2.

Conclusions
In this work, we made in-depth discussion on the application of finite element method to the P2D physic-based Li-ion cell model. Galerkin's weighted residual approach with quadratic basis functions was employed for our finite element analysis. The discretized field variables are defined into two categories: node variables and element variables, and the conversion between these two types of variables can be performed through a simple linear transform. The Butler-Volmer equations and the solid phase diffusion equations are solved in each element rather than at each node, which significantly lowers the computation load in simulation. We also developed simulation approaches for different battery control modes (current, voltage, load, and power control modes), and identified the limitations for a fully linearized model in the power control mode. For a linearized model, the Butler-Volmer equations are approximated using a 1 st order Taylor series expansion, and the interface kinetic resistance is independent from the current density; when high power is applied, the linear polarization causes electrode potentials to drop sharply with increased current, and the cell might not meet the power request. Therefore, nonlinear compensation through are iterative Newton loop is important for the accuracy of model, because for nonlinear Butler-Volmer equations, the kinetic polarization drops exponentially with current density which prevents the cell voltage from being too low under high power. The model was also validated by test data with mixed current-voltage control modes and key parameters for the cell materials were estimated using least squares approach.

Appendix A
The basic cell design parameters and physical properties are listed in Table A1. The electrical conductivity of bulk electrolyte is a function of electrolyte concentration and temperature: The effective diffusivity and electrical conductivity for electrolyte is defined as follows: ) unless CC License in place (see abstract The open circuit potentials for anode and cathode materials (U ) are functions of the surface state-of-charge of solid phase (θ * ), and the U vs θ * profiles are presented in Figure A1. Figure A1. Open circuit potential profiles for active materials: (a) anode material; (b) cathode material.

Appendix B
The surface state-of-charge for solid phase, θ * i (t), is calculated as follows: where θ avg,i (t) is the particle-average state-of-charge for the solid phase, and Q i,k (t) (k = 1, 2, 3, 4) are eigenfunctions. The governing equation for θ avg,i (t) is as follows: where c max is the maximum concentration of the solid phase, and R s is the radius of the solid phase particles. The governing equation for Q i (x, t) (i = 1, 2, 3, 4) is as follows: where D s is the solid phase diffusivity, and values for coefficients a k and b k (k = 1, 2, 3, 4) are listed in Table B1. The details for the derivation of these equations are provided in Ref. 15. As shown above, the dependent variables for the solid phase diffusion are associated with the surface electrochemical current density j s (x, t) in each element, and therefore are defined as element variables.

Appendix C
According to the definition in Equation 23, the following rule can be derived for element basis function q i (x): therefore the product for q(x)q τ (x) yields a diagonal matrix as follow: The integral of q i (x) over the cell geometry equals to the length of element i : The integral of q i (x)w(x, t) over the cell geometry equals to the integral w(x, t) of over element i : therefore the integral of vector q(x)w(x, t) is expressed as follow: [C6] and the element-average vectorŵ(t) can be calculated as follow: therefore the linear transform from w(t) toŵ(t) can be expressed as follow: where the operator H is given by: