Modeling of a Pouch Lithium Ion Battery Using a Distributed Parameter Equivalent Circuit for Internal Non-Uniformity Analysis

Dafen Chen 1,2, Jiuchun Jiang 1,2,*, Xue Li 1,2, Zhanguo Wang 1,2 and Weige Zhang 1,2 1 National Active Distribution Network Technology Research Center, Beijing Jiaotong University, No. 3 Shang Yuan Cun, Haidian District, Beijing 100044, China; 11117360@bjtu.edu.cn (D.C.); 13117372@bjtu.edu.cn (X.L.); zhgwang@bjtu.edu.cn (Z.W.); wgzhang@bjtu.edu.cn (W.Z.) 2 Collaborative Innovation Center of Electric Vehicles in Beijing, Beijing Jiaotong University, No. 3 Shang Yuan Cun, Haidian District, Beijing 100044, China * Correspondence: jcjiang@bjtu.edu.cn; Tel.: +86-10-5168-4056; Fax: +86-10-5168-3907


Introduction
The lithium ion battery is one of the most promising candidates for the energy storage system (ESS) in electrical vehicles.To reduce the cost, prolong the life and ensure the safety of the ESS, the optimal design of the battery cell and pack with robust battery management system (BMS) is essential.Two kinds of battery models are proposed to fulfill the optimal design and to be implanted in the BMS, respectively.The first kind of model is physics based.A model based on the theory of intercalation electrodes and concentrated solutions was proposed by Doyle et al. [1].Lou et.al.[2] built an extended single particle model for higher rates simulation.Kim et al. [3,4] developed a multi-scale multi-domain model framework, which resolved electrochemical, thermal and electrical coupled physics at varied length scales in the lithium ion battery.Three-dimensional electrochemical-thermal models were proposed and widely used by researchers [5][6][7][8] to analyze heat generation and temperature distribution in the lithium ion battery.A two-dimensional model was reported by Shin's group [9][10][11] for the thermal behavior and scale-up design of the lithium-polymer battery.These models are usually computationally complex and time consuming.In addition, a large number of parameters is needed for model inputs, and some of the parameters are extremely hard to obtain.Thus, such models are only suitable for cell designing and analyzing.For example, a three-dimensional electrochemical model was proposed by Kohneh et al. [12] and used to examine the effect of parameters, such as current collector thickness and tab location, on reducing non-uniform voltage and current distribution in the cell.Liu et al. [13] developed an integrated computational method that considered the mechanical, electrochemical and thermal behaviors of the lithium ion battery to study the nail penetration problem.Samba et al. [14] used the electrical-thermal model to design the battery thermal management system.Though simplified approaches, such as the simplified electrochemical multi-particle model and homogenous pseudo two-dimensional model, were developed by Mastali et al. [15] to decrease the computational time, the speed and simplicity of three-dimensional electrochemical-thermal models are still of concern.The second kind of models is the equivalent circuit model (ECM).Here [16][17][18][19][20], the battery is usually regarded as a mass point.Only a few parameters are in the model and can be derived from the external characteristics of the cell [21,22].Therefore, they are suitable to be implanted in the battery management system (BMS) for the state of charge or the state of health estimation [23][24][25][26].Lumped-parameter thermal models were added to ECM to predict the thermal characteristics of the cell by Lin et al. [27] and Forgez et al. [28], which made the model more comprehensive.
Meanwhile, in order to reach higher power and energy density, the size of the battery cell is growing.This leads to the imbalance of the potential, current density, temperature and state of charge (SOC) becoming significant, so the battery can no longer be treated as a mass point.Large SOC differences up to 5.3% were reported by Fleckenstein et al. [29] in a LiFePO 4 cell.Taking the computational cost into consideration, a new model to predict the internal imbalance of the battery is desired at this stage.Therefore, the distributed parameter equivalent circuit model (DPECM) is proposed in our former works [30].However, the former proposed model could only solve the electric related variables, and the model is only validated by the terminal voltage.As the temperature of lithium ion batteries is vital for their life span, safety and performance, thermal models need to be established and implemented into the former proposed model, to analyze the distribution and variation trend of internal current, voltage and temperature.
On the other hand, the validation of those models that could obtain the internal variables is somehow technically difficulty, because the distribution of the internal potential, current density or SOC cannot be measured directly.Wang's group [31,32] proposed an in situ measurement method for current distribution in a lithium ion battery.The structure of the cell was changed because the positive electrode was segmented along the length of the electrode sheet; while temperature is much easier to measure compared to other variables inside a battery.Forgez et al. [28] and Li et al. [33] successfully embedded thermocouples into cylindrical and pouch cells, respectively.They proved that the cell characteristics affected by the embedded sensors were negligible.To validate the proposed model in our research, a customized cell with nine T-type thermocouples embedded inside the battery was provided by the battery manufacture.The customized battery was tested for the parameter identification and validation of the temperature distribution inside the battery.
The remaining parts of this paper are arranged as follows.Section 2 describes the experimental setup, cell dimensions and temperature sensor locations.Section 3 introduces the modeling method and validation.The internal non-uniformity is discussed in Section 4. The conclusions are summarized in Section 5.

Experiments
The cell studied in this paper is a 35 Ah LiNi 1/3 Co 1/3 Mn 1/3 O 2 (NCM) pouch cell, which has a total of 36 laminated layers.Each laminated layer is composed of positive foil, positive material, separator, negative material, and negative foil, and the characteristics of each component are presented in Table 1.The top view from cell thick direction and the cell dimensions are shown in Figure 1.The length and width of the main body are 195 mm and 165 mm, respectively.Positive and negative tabs locate at two opposite sides of the cell, whose length and width are 40 mm and 80 mm, respectively.The total thickness of the cell is only 14 mm, which is much smaller than the length and width of battery, so the differences of local SOC, current density, and potential along thick direction are ignored in this paper.In order to acquire internal temperature reliably, nine T-type thermocouples are embedded inside the cell during the manufacturing process of the cell using the same method proposed by Li et al. [33].The thermocouples were sandwiched by two pieces of separators and then placed in the middle of the thickness before vacuuming and sealing of the cell.Nine thermocouples are uniformly distributed on the middle plane.

Ly=195mm
Lx=165mm Lty=80mm T( T(3,3) In order to obtain the model input parameters and verify the accuracy of the model, two sets of experiments are executed.One is for parameter estimation, and the other is for validation.

Parameters Identification of the Lumped First-Order Resistor-Capacitor (RC) Model
A typical lumped first-order resistor-capacitor (RC) model is shown in Figure 2. It consists of a voltage source, a resistor and a parallel branch of the resistor and capacitor.The voltage source represents the open circuit voltage (OCV) of the cell.Resistor R Ω represents the ohmic resistance in the cell.The parallel branch is utilized to model the chemical diffusion in the cell.In order to obtain the parameters of the first-order RC lumped model, pulse discharge is a widely-used method for lithium ion batteries.In this paper, the cell was put into an environment chamber during the test to keep the ambient temperature constant, and then, discharge pulses were applied to the cell by battery test equipment: Arbin BT-2000.Each discharge pulse depleted 5% of the total capacity and was followed by 2 hours of rest.To make the parameters suitable for different current rates, the discharge pulse at each SOC contained multiple current rates, including the current fragments of 0.5C, 1C and 2C, as shown in Figure 3a.The corresponding voltage response of the cell is shown in Figure 3b.The lumped first-order RC model parameters can be estimated according to the current and voltage curve using the method proposed by Huria et al. [34,35].The parameters were determined using the parameter estimation tool in Simulink Design Optimization.The tool iteratively simulated the discharge voltage profile while comparing the simulation results with experimental data.The nonlinear least squares algorithm was used to minimize the sum of square error.The fitting result of the discharge pulse at SOC = 0.55 is presented in Figure 4.The fitted results agree well with the experimental data.The estimated parameters of U ocv , R Ω , R ct and C ct are shown in Figure 5.

Entropy Change Measurement
The entropy change of the lithium ion battery can describe the generated reversible heat.It can be calculated using Equation ( 1): where ∆S is the entropy change of the battery; F is the Faraday constant; U is the equilibrium potential of the battery; T is the temperature.The entropy change could be calculated through the OCV variation with temperature.The method proposed by Forgez et al. [28] was used in this paper to obtain the coefficient of OCV variation with temperature.The cell was discharged to an objective SOC and rested to the steady state.Then, the temperature of the environment chamber was changed to different values, as shown in Figure 6a.In order to extract the temperature coefficient from these data, the voltage was fitted by the function V(t, T) = A + BT + Ct with A, B and C as constants and B corresponding to the temperature coefficient ∂U/∂T.The fitting results for ∂U/∂T estimation are shown in Figure 6b.The ∂U/∂T results are shown in Figure 7.

Validation Experiments
In order to analyze the voltage and temperature characteristics, constant current discharge at 1C, 2C and 3C rates was performed using Arbin BT-2000.The temperature was measured by an Agilent 34901A.The voltage and temperature during the discharge process with different rates are shown in Figure 8a. Figure 8b shows the temperature measured by the internal thermocouples.The thermocouples at T(1, 1) were broken, so the data are not shown here.At the end of discharge, temperature at the center of the cell (T(2, 2)) is the highest.The center temperature close to the positive tab (T(1, 2)) is higher than that close to the negative tab (T(3, 2)), which might be caused by the different properties of positive and negative tabs.Although the geometries of positive and negative tabs are the same, the positive aluminum tab has a smaller heat conductivity coefficient and larger electrical resistivity compared to the negative Copper tab.Therefore, more heat is accumulated around the positive tab, leading to a higher temperature around positive tab.The temperature values of T(2, 1) and T(2, 3) are almost the same, owning to the symmetric cell structure.

Modeling and Validation
The model is comprised of two parts.One is the electrical part, which is used to solve the electrical-related variables, such as the potential, current density and SOC.The other is the thermal part to solve the thermal-related variables, including heat generation, conduction and dissipation.

Electrical Model
The schematic of the electrical part in DPECM is shown in Figure 9.According to the layered structure and electrochemical principles of pouch batteries [36][37][38], the proposed model consists of five parts, including the positive and negative tabs, the positive and negative foils and the main material body.The current flow of positive and negative tabs only exists in the horizontal plate.The positive and negative foils not only have current flowing in the horizontal plate, but also have current pouring into or pumping out from the sub-models.The main material body, including the positive and negative materials, electrolyte and separator, etc., is divided into a group of sub-models.Each sub-model is represented by a first-order RC equivalent circuit.The first-order RC model is adopted in this paper, because Hu et al. [39] indicated that the first-order RC equivalent circuit model was the most suitable candidate for an NCM cell after comparing twelve battery models.One branch of the sub-models was magnified in the right side of Figure 9 to show the details of first-order RC model.It consists of a voltage source, a resistor and a parallel branch of the resistor and capacitor.The current of the sub-model located at the i-th row and j-th column in Figure 9 could be calculated as Equation (2): where the subscript (i, j) represents the node location; i is the current of the sub-model; e is the equilibrium potential of the sub-model; ϕ − is the node potential on the negative foil connected with the sub-model; ϕ + is the node potential on the positive foil connected with the sub-model; r Ω is the ohmic resistance of the sub-model; u p is the polarization voltage of the sub-model, which could be calculated as: where u p is the initial polarization voltage; r p is the polarization resistance of the sub-model; c p is the polarization capacitor of the sub-model.The resistance value r p in the sub-model is calculated using the equivalent resistance of the lumped first-order RC model multiplied by sub-model numbers.The capacitance value c p in the sub-model is calculated using the equivalent capacitance of the lumped first-order RC model divided by sub-model numbers.
In this paper, u p is iteratively updated.Therefore, Equation (3) needs to be discretized.The discretized result of Equation (3) using the Euler method is: where ∆t is the step time; u p,t is the polarization voltage at instant t; u p,t+∆t is the polarization voltage at instant t + ∆t; r p,t , c p,t and i t are the polarization resistance, polarization capacitance and current of the sub-model at instant t, respectively.Assume that there are N nodes in total on the positive foil and tab, as shown in Figure 9.According to Kirchhoff's current law, the sum of currents flowing through each node is zero.Therefore, a typical node equation is: where the subscript (i, j) represents the node location; ϕ + (i,j) is the potential of the node located at (i, j); r +(i−1,j)(i,j) , r +(i+1,j)(i,j) , r +(i,j−1)(i,j) and r +(i,j+1)(i,j) are resistances between two adjacent nodes.As the locations of all nodes are even in this paper, all of the resistances between two adjacent nodes are equal on the positive foil.Therefore, Equation ( 5) could be simplified as: where r + represents the resistance between two adjacent nodes.Define a connectivity function as: The relationship between the current of all sub-models and the potential of all nodes located on the positive foil and tab is: where the node potential equations of the positive foil and tab could be written as: where A + is an N × N matrix representing the node connectivities of the positive foil and tab; Φ + is an N-dimensional node potential vector of the positive foil and tab; I is an N-dimensional sub-model current vector.Those nodes on the positive tab do not have sub-models, so the currents of the sub-models in the corresponding places of I are set to zero.Similarly, the node potential equations of the negative foil and tab could be deduced as: where A − is an N × N matrix representing node connectivities of the negative foil and tab; Φ − is an N-dimensional node potential vector of the negative foil and tab; r − is the resistance between two adjacent nodes on the negative foil.Those nodes on the negative tab do not have sub-models, so the currents of the sub-models in the corresponding places of I are set to zero.
The current of all sub-models could be organized as: where E is the N-dimensional equilibrium potential vector of sub-models; V P is the polarization voltage vector of sub-models; R Ω is an N × N-dimensional diagonal matrix of sub-model ohmic resistances that: The resistance is set to be infinity for those nodes that do not have sub-models.Define: In order to solve the coupled sub-model currents and node potentials, combine Equations ( 9)- (11), and the following equation is obtained: It could be simplified as: In order to solve Equation ( 15) to obtain Φ + and Φ − , boundary conditions and initial values are needed.In this paper, the boundary conditions are: the node potentials at the edge of the negative tab are set to zero; the nodes at the edge of the positive tab share the external current evenly.The initial value of V p is set to 0. The initial value of each element in E is equal to the equilibrium potential when SOC = 1.The initial values of r Ω , r p and c p are also set to be the corresponding values when SOC = 1.The currents of sub-models I could be obtained by substituting the Φ + and Φ − into Equation (11).
The current flow in the positive foil and tab could be calculated as: where J + is the current flow in the positive foil and tab; γ + is the electrical conductivity of the positive foil and tab.Similarly, the current flow in the negative foil and tab could be calculated as: where J − is the current flow in the negative foil and tab; γ − is the electrical conductivity of the negative foil and tab.

Thermal Model
The schematic of the thermal part in DPECM is shown in Figure 10.Each sub-model consists of a parallel branch of the current source and capacitor, a resistor and a voltage reference [28].The current source represents the heat generation source.The capacitor represents the heat capacitance of the sub-model.The voltage reference provides the ambient temperature of the sub-model.The resistor represents the resistance of heat conduction and convention from the sub-model to the ambient environment.The heat conduction in the cell is modeled by the resistors between two adjacent nodes.The values of the resistors are calculated through the thermal properties of each layer using the method proposed by Chen et al. [40].
The heat generation of each sub-model is composed of irreversible and reversible heat as Equation ( 18): where qg,i is the heat generation rate of the i-th sub-model; T is the temperature of the sub-model.Assume that the number of thermal nodes in Figure 10 is M. As the positive and negative tabs do not have connected sub-models, where the generated heat of those nodes are set to zero, so the heat generation rate of all of the sub-models could be described as: The heat conducted from other nodes could be calculated as: where Qc is an M-dimensional vector representing the conductive heat; T is an M-dimensional vector representing temperatures of all nodes; A t is an M × M matrix representing the thermal connectivities of nodes in the thermal model, which is defined as: .., M and a = b) is the connectivity function of the nodes in the thermal model.
The heat generation rate of the i-th node in the positive foil and tab is: Therefore, the heat generation rate of all of the nodes on the positive foil and tab is: Similarly, the heat generation rate of all of the nodes on the negative foil and tab is: The temperatures of all of the nodes can be calculated as Equation ( 24): where T 0 is the initial temperature vector; ∆t is the step time; m is the equivalent mass of the battery block corresponding to the sub-model; C is the equivalent specific heat capacity of the battery.

Solution
The above proposed model is built and solved using MATLAB language (M-code).The solution of the potential and current density of each sub-model can be obtained by Equations (11) and (15).The advantage of this method is that the variables in the electrical part are coupled directly.Then, the V p of the next step is predicted.At the same time, the outputs of the electrical model are used to calculate the heat generation rate, so the thermal part of the model can be solved.The flow chart of the solving method is shown in Figure 11.

Results and Discussion
The simulated and experimental temperatures at the cell center (T(2, 2)) and positive side center (T(1, 2)) during 1C, 2C and 3C discharge are shown in Figure 12.The simulated temperature conforms with the experimental data very well, and the simulated temperature error is less than 1 °C.The accurate temperature results indicate that the proposed model has the capability of predicting internal situations in the cell.

Non-Uniformity Distribution at the Initial State
During discharge, according to the structure shown in Figure 9, the current from sub-models flows through the positive foil and then is collected at the positive tab end.After the current passes the external circuit and comes back to the negative tab end, it is allocated to sub-models through the negative foil.At the beginning of battery charge/discharge, SOC and temperature are uniform in all sub-models since the battery is rested to steady state, but the current distribution becomes non-uniform when the current flows to relatively narrow tabs.Figure 13a,b shows the current density in the positive and negative current collectors, respectively.The direction of arrows represents the current direction, and the length of arrows represents the magnitude of the current flow.Therefore, the current close to tabs is larger than that in the opposite side.Figure 13c,d shows the potential of the positive and negative collectors, respectively.Since the thickness of each foil is uniform, large potential drops are produced close to tabs because of larger current flow in those areas.Figure 13e shows the current density in sub-models.The potential differences between two end nodes (ϕ + − ϕ − ) of sub-models close to tabs are smaller than that of the central part.At the same time, the polarization voltages of all of the sub-models are zero, and the equivalent potentials of all of the sub-models are equal at the beginning of the discharge.Therefore, the current flows through the sub-models, decided by Equation ( 2), close to tabs are higher at the beginning of the discharge process.The temperature close to tabs, shown in Figure 13f, is higher due to higher current density in foils and sub-models.

Non-Uniformity Evolution During Discharge
The non-uniformity becomes more complex after the initial state.It changes with the current profile and ambient temperature.The distribution of the current density and temperature in the sub-models could change dramatically, which sometimes even turns into an inverse distribution compared with the initial state.
The distributions of current flowing through sub-models and temperature at the end of the 1C discharge are displayed in Figure 14. Figure 14a indicates that the current density close to tabs becomes smaller than those at the cell center, which is inverse compared to the situation at the initial state, as shown in Figure 13e.This is associated with the non-uniform current accumulation and potential change in the sub-models.The current of the sub-model is calculated by (e − (ϕ + − ϕ − ) − u p )/r Ω , so the current is decided by both the potential (e) and the voltage difference between positive and negative foils (ϕ + − ϕ − ).The voltage difference between positive and negative foils close to tabs is always smaller than that at the cell center during discharge because of the existence of foil resistance and the current flow direction.While the accumulation of larger current in the sub-model causes lower local SOC during the discharge process, the local SOC of the sub-models close to tabs becomes lower than that at the cell center.Therefore, the potential difference of sub-models close to tabs and the cell center increases gradually.The influence of the potential difference becomes stronger along with the accumulation of non-uniform current distribution and finally dominates the current of the sub-models, which inverses the current density distribution.Figure 14b indicates that the temperature at the cell center becomes the highest, which is mainly due to the poor heat dissipation at the center.To further illustrate the evolution of internal non-uniformity during discharge, three points (P(1,2), P(2,2) and P(3,2) in Figure 1) are selected to be analysis points.The current density and local SOC at those three points are calculated and presented in Figure 15. Figure 15a indicates that the difference of the local current density declines at the beginning of the discharge.The current at P(3,2) is the largest, and the current at P(2,2) is the smallest.The current becomes uniform at around SOC = 0.38.After the current crosses the point at SOC = 0.38, the current difference rises again, leading to the current at P(2,2) being the largest.The current difference decreases again after SOC around 0.2.
The current non-uniformity will cause a local SOC difference of each sub-model as shown in Figure 15b.The local SOC difference increases before the current cross point and decreases later, which can be explained by its definition.The SOC, defined as: is decided by the initial state and current flow through the sub-model.The local SOC difference reaches the maximum value at the current cross point and then is eliminated by the inversed current density distribution.The appearance of the current cross point is caused by the relatively flat characteristics of de/dSOC and dr/dSOC between 0.8 > SOC > 0.2, as shown in Figure 16.The effects of resistance in foils are counteracted by the parameter difference caused by the local SOC difference at the current cross point.After the current cross point in Figure 15a, the increase of the de/dSOC difference becomes larger and dominates the current distribution, which enlarges the difference of the current density.It is worth noting that the direction of currents in sub-models remains the same, while the shape of the current density distribution changes from the bowl-shaped distribution to the peak-shaped distribution.The rapid increase in cell impedance after SOC < 0.2 makes the difference in current density decrease again.

Limitation of Technique
As mentioned in Section 1, the distributions of internal potential, current density or SOC can hardly be measured directly without changing the geometry and characteristics of the commercial lithium ion battery.Though the proposed model has the capability of predicting internal temperature, current density and potential, only the distributions of internal temperatures can be measured even using a customized battery.Therefore, only the internal temperatures predicted by the proposed model are validated at the present stage.

Conclusions
In this paper, a DPECM model consisting of a directly-coupled electrical model and a thermal model for a lithium ion battery is proposed.The model is simpler than the 3D electrical-chemical thermal model, and the parameters are easier to obtain, which is similar to the lumped ECM model.However, the proposed model has the capability of simulating the internal distribution of the current density, potential and temperature inside a cell.The model is validated by comparing the simulated temperature with measured results at test points in a customized cell that is embedded with nine T-type thermocouples.The simulated temperature error is less than 1 °C, indicating that the model can be applied to predict the internal characteristics of the battery.Summarizing the simulation results, the following conclusions can be obtained: 1.The initial non-uniformity of current density inside a 35-Ah NCM pouch cell is caused by the electrode foil resistance and relatively narrow tab. 2. The current flow through the sub-models close to tabs is higher than that at the center at the beginning of discharge, because the potential drops in positive and negative foils.3. A current cross point exists during the constant discharge process.The local SOC difference increases with the reduction of the average SOC until the average SOC arrives at that point and then decreases until the average SOC arrives at zero. 4. The current flow through the sub-models close to tabs becomes lower than that at the center at the end of discharge due to the accumulation of the local SOC difference, rapid OCV drop and rapid resistance increase.

Figure 3 .
Figure 3. Voltage and current profiles for parameter estimation.
Temperature results at 3C discharge

Figure 8 .
Figure 8. Voltage and temperature experimental results.
when node a and node b are connected) 0 (when node a and node b are disconnected) a, b ∈ 1, 2, ..., N and a = b b ∈ 1, 2, ..., N and a = b) is the connectivity function of the positive foil and tab; ϕ + n (n ∈ 1, 2, • • • , N) is the potential of the n-th node on the positive foil and tab; i n (n ∈ 1, 2, • • • , N) is the current of the n-th sub-model.Define:

Figure 11 .
Figure 11.The flow chart of the solving method.

Figure 12 .
Figure 12.Temperature validation at different discharge rates.

Figure 13 .
Figure 13.Non-uniformity at the beginning of the 1C discharge.

Figure 14 .
Figure 14.Non-uniformity at the end of the 1C discharge.

Table 1 .
Properties of each battery component.