Bridging physics-based and equivalent circuit models for lithium-ion batteries

In this article, a novel implementation of a widely used pseudo-two-dimensional (P2D) model for lithium-ion battery simulation is presented with a transmission line circuit structure. This implementation represents an interplay between physical and equivalent circuit models. The discharge processes of an NMC-graphite lithium-ion battery under different currents are simulated, and it is seen the results from the circuit model agree well with the results obtained from a physical simulation carried out in COMSOL Multiphysics, including both terminal voltage and concentration distributions. Finally we demonstrated how the circuit model can contribute to the understanding of the cell electrochemistry, exemplified by an analysis of the overpotential contributions by various processes.


Introduction
For lithium-ion batteries, mathematical models not only constitute tools to estimate the performance of different battery components, as well as the cell or the battery pack, but also provide tools to strengthen the understanding of many physical properties, which determine the electrochemical response during the battery operation. An adequate model can be used to both interpret experiment results [1] and offer estimations of quantities that cannot be easily accessed through measurement, for example local overpotentials, capacity losses, morphology changes [2], internal temperature fluctuations and the growth of the solid electrolyte interphase (SEI) layer [3,4,5]. Modelling can aid battery diagnostics, and thereby help to prevent premature ageing of the cells [6,7,8].
From atomistic scale to system scale, there exist battery models with different complexities depending on the application. In battery management systems, empirical circuit models are commonly used which can fairly accurately describe the dynamic behavior of batteries and can be computed in real time with on-board microprocessors. However, the quality of a circuit model is directly dependent on parameterization data which requires massive amount of experiments at different state of charges (SOC) and temperatures [9]. Moreover, the effectiveness of the model decreases with the battery being aged since all the parameters are affected by the state of health (SOH). It is also a challenge to transfer these models to different generations of cells or types of batteries, which make them reliable for only a limited segment of battery devices. In many cases, Kalman filters are used together with the empirical circuit models to improve the accuracy [10,11]. While empirical circuit models can illustrate the dynamic behavior of batteries with a transparent structure, such a structure gives very little or no information about the physical parameters actually governing the behavior of the system.
A physics-based approach can instead be employed using the first principles-based lithium-ion battery model that was developed by Newman, Doyle and Fuller [12,13] and has been implemented into a number of commercial softwares, e.g. COMSOL Multiphysics. Newman's model is a Pseudo-two-Dimensional (P2D) model consisting of a set of partial differential equations (PDEs). The PDEs are typically solved with the finite element method [14] or the finite volume method [15] which leads to a high computational load. As a consequence of the long computational time, there has been a great effort on simplifying the P2D model into single particle models [16,17] or reduced-order P2D models with polynomial approximations [18]. These simplifications increase the calculation efficiency but also limits the conditions where the models are valid [19]. Furthermore, the physics-based models with complex mathematical equations remains challenging to use and interpret for researchers without a solid mathematics or physics background.
To reduce the gap between the empirical equivalent circuit models and the physics-based models, different physics-based circuit models have recently been implemented in several different ways. One simple approach is to use the transmission line structure with an elementary resistor network but without any link to the mass transport process in the cell [20,21]. One further step, employed by Sato et al., was to connect the elements in the circuit model with physical principles [22], but then the current distribution within the electrode is ignored in their model. A recent interesting approach used by Li et al. [23] is to construct a physics-based equivalent circuit model with passive electrical components, which is fast and accurate, however, the number of circuit elements, voltage sources and transformers make the model somewhat less easy to comprehend and implement.
So far, the physics-based circuit models reported in literature are either simplified [20,21,22] and thus cannot fully capture the battery cell behaviors described in Newman's P2D model, or with an advanced structure [23], somewhat too much exceeding the simplicity of the usual transmission line models. These facts are crucial flaws, and thus forms an important research gap, i.e. to have an easy understandable, but still accurate, circuit-based model of a Li-Ion cell.
To bridge this gap, in this work we present a novel implementation for Newman's P2D model with a transmission line circuit structure. The circuit structure not only solves the current distribution with the mesh current method, but also offers a clear visual illustration of the P2D model. In this implementation, the physics-based electrochemical equations are kept without simplifications or approximations, and thus guarantees the validity of the model. Apart from the demonstration of the novel implementation, one additional purpose is to quantify its accuracy towards a 'full physical model', in our case in COMSOL Multiphysics, in terms of step time and number of meshes.

Model implementation
The discharge process in a lithium ion battery cell is described in Fig. 1, where lithium ions move from the negative electrode to the positive electrode inside the battery and electrons move from the negative electrode to the positive electrode through the outer circuit In the charge process, the opposite flows will occur. The basic steps shown in the figure include electrons moving in the solid phase, charge transfer and mass transport of lithium ions in both solid particles and in the liquid electrolyte. These processes can be modelled by a set of equations, which are listed in Table 1 [13]. In this work, these processes are represented with the electronic resistance in the solid , charge transfer resistance , electrode potential U and electrolyte resistance , which are connected in a transmission line structure. The transmission line structure was originally proposed by de Levie [24] assuming that the porous electrode consists of cylindrical particles. The transmission line structure can, however, be generally applied for a porous electrode if only the current distribution in the through-plane direction is considered, regardless of the particle shapes. The electronic conduction follows Ohm's law in (2) and is modelled with a resistor in the circuit. The value of only depends on the electronic conductivity σs and the volume fraction of the solid matrix, so will remain the same during the simulation unless the model is coupled with, for example, a temperature or aging phenomenon. In commercial lithium-ion batteries, the electronic conductivity is improved by adding different types of carbon additives in the electrode, and can therefore in many cases be ignored. The charge transfer process is described by the Butler-Volmer equation in (4) and is represented by in the circuit. Under a very small current, i.e. close to equilibrium, the Butler-Volmer equation can be linearized from a Taylor series and thus where is the charge transfer current density per surface area, 0 is the exchange current density, is the Faraday constant, is the gas constant, is temperature and η is the overpotential caused by the redox reaction. The linearization yields the charge transfer resistance for the active surface area ′ as Table 1: Equations used in Newman's model.
The charge transfer resistance for the electrode area is scaled with the specific surface area and the length of the mesh element ℎ where can be estimated with the volume fraction of the solid and the particle radius rs The linearization above is only used to initialize at the equilibrium state. During the simulation, will be updated according to (4) with the corresponding local current density and lithium ion concentrations.
After the charge transfer process, a concentration gradient will be built up both in the particles and in the electrolyte. The concentration gradient together with the potential gradient are the driving forces for the mass transport process. The mass transport in the particle is expressed by Fick's law in (8) and the surface concentration , will determine the electrode potential in the circuit. Finally, represents the resistance of the mass transport in the electrolyte, which is described with the concentrated electrolyte theory in (3) and (7). In the equilibrium state, there is no concentration gradient and ∇ = 0 meaning that can be initialized with only the electrolyte conductivity . Later in the simulation, will be updated according to (3) based on the electrolyte current and concentration. All the transport parameter values are corrected with (9) thereby taking the porosity of the electrode into consideration In this work, a decoupled quasi-dynamic simulation approach is implemented, where the current distribution is solved in an algebraic way and is considered to be static within the time step. With the formed static boundary conditions, the concentrations are solved dynamically, as shown in Fig. 2. The current distribution is solved with the mesh current method by the circuit structure, and the concentration distribution is solved with the finite difference method. The two methods are explained in the sections below.

Solving the current distribution
In a porous electrode, the current is distributed unevenly within the electrode, and the current density is normally higher close to the separator and lower close to the current collector. In this work, the current distribution is described with the transmission line circuit structure and solved with the mesh current method, instead of the (5) and (6). After initializing the component values in the circuit through (1)-(4), the resistance triangular matrix and voltage vector for the transmission line structure in Fig. 1 can be generated as where is the input current density. The current distribution can be calculated as where the elements in the current vector refer to the mesh currents 1 , 2 , ..., 1 −1 shown in Fig. 3. Thereby , and can be calculated with Kirchhoff's circuit laws. The resulting and will be used as input and boundary conditions in (7) and (8)

Solving the concentration distribution
Within one time step, the current solved from the mesh current method is assumed to be constant, and the PDEs (7) and (8) can be solved numerically to obtain the concentration distribution. Both (7) and (8) are parabolic PDEs in one space dimension with Neumann boundary conditions. We have chosen the finite difference method for the spatial discretization because of the simple geometric feature for one dimensional problems. Other numerical methods, such as the finite element method, the finite volume method, and the difference potential method can however also be used [25].
There are two distinct difficulties in solving these two PDEs. In (7), the material porosity is discontinuous at the two material interfaces. For an accurate spatial discretization, the finite difference stencils should not cross the material discontinuities. In our method, we discretize (7) in each subdomain (negative electrode, separator and positive electrode) separately by finite difference operators with a summation-by-parts (SBP) property [26]. At the material interfaces, we impose physical interface conditions such that the concentration and its flux are continuous. These interface conditions as well as the Neumann boundary conditions are imposed numerically by the simultaneous-approximation-term (SAT) method [27]. As a result, the semi-discretized equations are energy stable. The resulting system of ordinary differential equations is solved by the MATLAB built-in function ode15s.
Equation (8) is reduced from the three dimensional heat equation in spherical coordinates when the solution is independent of the polar angle and azimuthal angle. The numerical difficulty here is the singularity at = 0. In our method, we first multiply (8) by on both sides and then approximate the spatial derivatives in (18) by the SBP finite difference operators [28]. We again use the SAT method to impose the Neumann boundary conditions. In this case, the energy stable semi-discretized equations are a system of differential algebraic equations, and are also in this case solved by using the MATLAB built-in function ode15s.
After the current and concentration distributions are established, the resistance matrix will in the next step be updated according to (1)-(4) with the updated current and concentration values.

Comparison to a physics-based model
With the implementation described above, an NMC-graphite lithium-ion batterya widely used commercial battery chemistry especially for electric vehiclesis simulated at different discharge currents. The parameters used to describe the cell and its components are listed in Table 2 [29]. A simulation with the same cell parameters is performed with COMSOL Multiphysics as the comparison reference.
The oxidation reaction on the negative electrode side is and on the positive electrode side the reduction reaction is The simulated battery voltage profiles during discharge are presented in Fig. 5a, showing an excellent agreement with the results from COMSOL Multiphysics. The differences are shown in Fig.  5b in percentage. The difference is somewhat larger at higher currents, but still indicates a high accuracy for the transmission line-based model. Besides the terminal voltage, the lithium ion concentration distributions are also critical in the simulation as the concentration is one key factor if aging phenomenon or mechanical stress are to be coupled in the model [30,31]. The concentrations in the electrolyte is shown in Fig. 6a. At t = 0 s, the initial concentration in the electrolyte is 1 mol/dm 3 . At t > 0 s, the battery starts to discharge and lithium ions travel from the negative electrode to the positive electrode and thus a concentration gradient starts to build up. This concentration gradient will quickly reach an equilibrium state in the beginning and then stay constant during the rest of the discharge time. Fig. 6b shows how the concentration at the particle surface changes with time. The initial lithium ion concentration is 10% in the negative electrode and 90% in the positive electrode which corresponds to a high SOC. During discharge, the negative electrode is being de-lithiated and the surface concentration decreases. As the current tends take the path of least resistance, the charge transfer current density is higher at the separator side (46.6 µm) than the current collector side (0 µm). Therefore the surface concentration decreases faster at the separator side. The opposite process happens in the positive electrode, i.e. the particles are lithiated during discharge and again with a higher current density at the separator side. For the results shown in Fig. 5a, 20 mesh elements are used in the negative electrode, separator and positive electrode, respectively. In general, a finer mesh can give numerically more accurate results but on the other hand increases the computation time. In Fig. 7a the results with different meshing are presented using a discharge current of 2 mA/cm 2 . For the case of 20 meshing elements, a comparison of 2 nd order and 4 th order space discretization is also shown. As can be seen, the 4 th order space discretization provides a higher accuracy especially when the voltage profile starts to change rapidly at a low state of charge. However a minimum number of 8 mesh elements (9 grid points) is required in each domain to implement the 4 th order space discretization. Another factor that affects the accuracy and simulation time is the choice of the time step, especially when a quasi-dynamic approach is used, as in this work. For the results in Fig. 5a, the time step was 5 seconds and Fig. 7b shows the results for different choices of time steps. This indicates that under a constant current, the time step does not affect the result significantly.

Contributions to the overpotential
This current implementation with a transmission line structure serves as an interplay between physics-based models and equivalent circuit models. On one hand, it is built upon the first-principle electrochemical equations with physically based material parameters. On the other hand, it contains an illustrative yet simplistic structure that can be easily interpreted and straightforwardly solved. We exemplify here how to this transmission line circuit model can aid the understanding of the contributions of different processes to the overpotential.
In the previous calculations, the resistance elements , and represent the processes of electronic conduction, redox reaction and the mass transport in the electrolyte. The diffusion process in the solid has not been directly introduced, since the electrode potential is determined by the particle surface concentration, and which is reflected in the circuit. Here, however, we introduce the term , to represent the diffusion process in the solid, as shown in Fig. 8. With the presence of the resistance , , the voltage source is determined from the average concentration in the active material particles instead of from the surface concentration. A corresponding battery model is built using the parameters in Table 2 and discharged from 4.15 V with 2 mA/cm 2 . When the battery is discharged to 3.6 V, the values of all the resistance elements are presented in Fig. 9 (note that the y-axis is in logarithmic scale). One observation is that the resistance caused by the mass transport in the electrolyte is a few orders of magnitude lower than the resistance caused by the diffusion in the solid, as the diffusion coefficient in the liquid is much higher than the diffusion coefficient in the solid ( ≫ ). However the overpotential caused by mass transport in the electrolyte is significant, as demonstrated in Fig. 10.  One way to interpret this result is to compare the diffusion length. With 20 mesh elements in the positive electrode, the diffusion length is around 2 µm for both the liquid and solid phases, and the diffusion coefficient affects the resistance values ( and , ). However, the total diffusion length for the electrolyte phase is the electrode thickness, which is around 20 times longer than the diffusion length in the solid. This, in turn, leads to the overpotential being significant. Similarly, another interpretation is to focus on the electrical circuit structure. The solid diffusion resistance elements , are connected in parallel, and only the branch current goes through each resistor. Contrarily, are connected in series, and the overpotential on each resistance element thereby ultimately adds up. As a consequence, the total overpotential becomes comparable to the overpotential caused by , . This example shows that although the diffusion in the liquid is way faster than in the solid, the resulting overpotential caused by the electrolyte resistance cannot be ignored. Especially in the case of thick electrodes, the mass transport process in the electrolyte could be a limiting factor.

Conclusion
A new implementation of Newman's P2D model in MATLAB is presented in this work, constructed as an interplay between the physical and equivalent circuit models. A classic transmission line structure is used to replace the equations that describe the current distribution within the electrode. The concentration distributions are solved with the finite difference method and a decouple quasidynamic approach is used to combine the two solutions. The results from the circuit model agree very well with the result simulated with a commercial software COMSOL Multiphysics based on the finite element methods. This implementation closes the gap between physical and equivalent circuit models and is a useful tool to understand the processes inside the battery, as exemplified by the analysis of different contributions to the overpotential.