Simplified Electrochemical and Thermal Model of LiFePO₄-Graphite Li-Ion Batteries for Fast Charge Applications

Literature on electrochemical and thermal modeling of battery systems is quite extensive. In the porous electrode theory, the electrode is treated as a superposition of two continua, namely the electrolytic phase and the solid matrix. The solid matrix is modeled as microscopic spherical particles, where the Li ions diffuse and react at the surface of the spheres. A classical 1D representation of a Li-ion battery is shown in Fig. 1. The complete electrochemical system is composed of three porous media, namely the negative electrode, the separator and the positive electrode. The porosity of the three regions is filled by the electrolyte liquid phase. Considering the LiFePO₄-graphite system, the electrochemical storage reactions in charge can be represented by Eqs. 1–2.

During discharge operating conditions, the reverse reactions occur in the cell. The nominal physical and chemical phenomena occurring in Li-ion systems can be expressed by the mass conservation of Li⁺ species (Eqs. 3–4), charge conservation (Eqs. 5–6) and electrochemical kinetics (Eqs. 7–8) that are presented in Table I. All these equations are the framework of the P2D electrochemical model. Non isothermal electrochemical models have also been addressed by Gu and Wang⁵ in a 1D model at the electrode level whereas other researchers have introduced a simple lumped-parameter thermal model for the full cell.^6,7 As presented in Table II, the energy balance (Eq. 14) integrates both generated (Eq. 15) and exchanged with environment (Eq. 16) thermal fluxes. Generated thermal flux takes into account both the irreversible and reversible contributions. Coupling between the electrochemical and thermal models is performed thanks to a classical Arrhenius law expressed by Eq. 17 for all mass transport and kinetic parameters. Readers can refer to the aforementioned literature for further detailed description of the mathematical developments. Based on the P2D mathematical structure, the theoretical set-up of the simplified model is detailed below.

Table I. 1D electrochemical model equations.⁷

Physical and chemical mechanisms	Eq.	Boundary conditions
Solid phase: conservation of Li+ species
Electrolyte phase: conservation of Li+ species
Solid phase: charge conservation
Electrolyte phase: charge conservation
Electrochemical kinetics
Electrode overpotential
Electrolyte ionic diffusivity
Electrolyte ionic conductivity
Electrolyte ionic diffusional conductivity
Solid phase electronic conductivity
Specific interfacial surface area

Table II. Heat transfer and energy balance equations.

Heat transfer and energy balance	Eq.
Energy balance
Thermal flux generated during operation
Transferred thermal flux to environment
Coupling between 1D electrochemical and lumped thermal models
Arrhenius law applied to mass transport and kinetic parameters Ψ

Zoom In Zoom Out Reset image size

Within the framework of electrochemical models, one has to consider the properties of electrode materials in detail. Indeed, even with the P2D model, non satisfactory results can be obtained if specific features of materials are not considered. Among all the electrode materials used in Li-ion batteries, LiFePO₄, which is a phase transformation material as initially reported by Padhi et al. in 1997,⁸ is one of the most promising chemistry due to its intrinsic safety, low cost and electrochemical performances. Nevertheless, LiFePO₄ material presents specific features like partial solid solution regions or differences in capacity restitution according to the charging or discharging path which makes the modeling of this material difficult. During the past decade, electrochemical modeling works on LiFePO₄ have been addressed to understand the complex mechanism of Li ions insertion/extraction that proceeds through a two-phase transition between a Li-poor Li_xFePO₄ phase (α-phase) and a Li-rich Li_yFePO₄ phase (β-phase). Srinivasan and Newman⁹ reported the first shrinking-core model, but in order to obtain a good agreement with experimental data on discharge curves at different current rates, they had to consider a particle size distribution (PSD) effect with spherical particles of two different sizes. The original shrinking-core model was then extended to the planar geometry by Wang et al.¹⁰ who obtained good results on discharge operating conditions up to 20C. The authors extended the original shrinking-core model by introducing the Li diffusion in both phases and the interface mobility as a model parameter. This model was then used to design new galvanostatic and potentiostatic intermittent titration techniques (GITT and PITT) to determine the diffusion coefficients of the respective phases in the composite (α+β).¹¹ Nevertheless, shrinking-core models can be difficult to program and issues related to the management of multiple interfaces can appear for complex charge/discharge profiles. Alternatively, Thorat et al.¹² proposed the concept of phase-change diffusivity to account for LiFePO₄ behavior in charge and discharge, introducing a concentration-dependent diffusion coefficient of Li in the solid phases. Moreover, the authors considered the resistive-reactant feature of insulating LiFePO₄, introducing the PSD effect with four different particle sizes. The resistive-reactant concept was further developed by Safari and Delacourt who studied, with a simplified electrochemical model, the (dis)charge path dependence, the asymmetry and the electronic resistance distribution of this material.¹³ Good agreement was obtained with experimental data until 1C charge and discharge rates. Later, Delacourt and Safari¹⁴ proposed a simplified mathematical model to explain the insertion/extraction mechanism of Li ions for LiFePO₄ material until 2C. In this simplified approach based on the SP model without considering PSD effect, the authors determined relationships between particle radius and current density to fit experimental charge and discharge curves. In this work, the radius dependencies upon current were tentatively explained with the mosaic concept.

In order to design a simplified electrochemical model, the AM developed by Di Domenico et al.⁴ was adopted and extended. The AM consists in neglecting the solid concentration distribution along the electrode and considering the material diffusion, for each electrode, inside a representative solid particle at which the Li concentration is equal to the mean value of the concentration over the electrode thickness, . Accordingly, the faradaic current density per unit volume, j_f, is equal to the mean value of the current density over the electrode thickness, . Therefore, the AM model considers a single particle for each electrode as the SP model but, instead of algebraic-differential equations in the SP model, the AM deals with Partial Differential Equations (PDE), well suited for the design of BMS functions. However, as previously mentioned, no distribution of the Li concentration in the electrolyte phase was considered in Ref. 4. In order to extend the AM mathematical structure, the mass conservation equations in the electrodes (Eq. 3) and in the electrolyte (Eq. 4) are used with the AM assumption, enabling the variation of Li concentration in the solid phase, , and in the liquid phase, c_e(z), to be computed at any time.

The current density per unit volume in each electrode satisfies the following spatial integrals:

where δ_n and δ_p are the thicknesses of the negative and positive electrodes, respectively. With this assumption, the charge conservation equations (Eqs. 5 and 6) can be solved analytically as a function of the spatial coordinate z, giving:

for the negative electrode region,

for the separator region, and

for the positive electrode region.

As described by Di Domenico et al.,⁴ the voltage of a battery at time t, V(t), is expressed as the sum of different terms, namely thermodynamic potentials, U, and overpotentials, η, appearing as current is passed through the system. According to Eq. 8:

In Ref. 4, the electrolyte concentration is uniform in the battery and, therefore, there is no diffusion overpotential in the electrolyte. Moreover, the thermodynamic potentials, U_p and U_n, of the positive and negative electrodes have been estimated from the Li concentration, c^s_s, at the surface of the single particle representing each electrode. In our work, the distribution of the electrolyte concentration is considered in order to deal with high load solicitations. Mass transport in both solid and liquid phases is taken into account by replacing the overpotential boundary values η(L) and η(0) by the average values of the overpotentials for each electrode. These overpotentials are due to the kinetic and diffusion phenomena occurring, respectively, at the solid-electrolyte interfaces and within the solid and electrolyte phases. In addition, U_p and U_n are estimated from the Li concentration, c^b_s, in the bulk (r = 0) of the single particle. The battery voltage can then be expressed as:

where θ^b_p and θ^b_n denote the normalized Li concentrations in the bulk of the positive and negative electrodes:

In these equations, c_s,n,max and c_s,p,max are the maximum Li concentration in the negative and positive electrodes. By the way, the normalized bulk concentrations θ^b_p and θ^b_n correspond to the coefficients x and y used in Eqs. 1 and 2.

To express the electrode overpotential, the approach used by Bergveld et al. was adopted in this work.¹⁵ The electrode overpotential η can be expressed as the sum of the kinetic and diffusion overpotentials at both electrodes:

This mathematical description of the charge transfer overpotentials facilitates the quantification, in first approximation, of the different contributions to the potential distribution in the battery. In the following, the kinetic and mass transport overpotentials are successively defined.

Firstly, to express the kinetic overpotential at each electrode as a function of the current intensity, the oxidation and reduction charge transfer coefficients α_ox and α_red are assumed to be equal to α = 0.5 for both electrodes. Then, considering Eq. 7 for the positive and negative electrodes with the assumption of the average model, and replacing the electrode overpotentials by the kinetic overpotentials, one gets:¹⁵

The positive and negative kinetic overpotentials can be expressed as:

where

In these equations, the exchange current densities are sometimes considered to be given by the following expressions:^1,2,4

However, for sake of simplicity and as performed in several reported battery models, the exchange current densities have been set as constants.^7,16

Secondly, the overpotentials due to mass transport (diffusion and migration) in the electrolyte and inside the solid active material have to be defined. In the electrolyte phase, the mass transport overpotential is obtained as the difference between the potentials of the electrolyte phase taken at the extremities of the electrodes. Eq. 23 expressed at z = L results in

Under the assumption of constant electrolyte concentration along the z axis, the last part of Eq. 38 can be seen as a representation of electro-migration inside the electrolyte. This term is generally known as ohmic resistance and can be determined thanks to electrochemical impedance spectroscopy for instance.¹⁷ This will be discussed in detail in the model calibration section. Neglecting the electric resistance of connectors and the conductivity of the solid phases in the electrodes, the ohmic resistance, R_ohm, can thus be expressed as a function of the geometric area of the electrodes, A, the thicknesses of the electrodes and separator, δ_n, δ_p, δ_sep, and the electrolyte effective conductivity, κ^eff, in the three regions:

The effective conductivities mentioned in Eq. 38 can be defined as a function of the Li concentration, c_e, and volume fractions, ɛ_e and ɛ_f, of the electrolyte and filler, and from the Bruggman exponent, Brugg, of each region, as follows:

In the solid phases, the mass transport overpotential is approximated as a function of the bulk, c^b_s, and surface, c^s_s, Li concentrations of the spherical particles. Neglecting migration in the solid phases, the diffusion overpotential in the negative and positive solid phases is expressed as in Bergveld et al.'s work:¹⁵

where θ^s_p and θ^s_n denote the normalized Li concentrations at the surface of the positive and negative electrodes:

The state of charge (SOC) of the full cell, SOC_bat, can be expressed as a function of the normalized bulk Li concentration in the positive electrode or in the negative electrode, for example:

with respect to the negative electrode, where θ^b_{n, 0%} and θ^b_{n, 100%} are the normalized Li concentrations in the bulk of the negative electrode, respectively for a fully discharged and fully charged cell.

Rearranging Eq. 26 with Eqs. 28–29, 32–33, 43–44 gives

or, with Eqs. 34–35, 45:

Interestingly, one can notice that the cell voltage is directly correlated to the specific design parameters of the system such as the electrode porosities, or to the physical parameters such as the electrolyte conductivity or the active material concentrations. This expression, which distinguishes the main electrochemical phenomena, could be useful in discussing and quantifying the impact of the variations of design parameters during aging.

Li-ion battery electrodes are often phase transition materials leading to the presence of multiphase regions as lithiation / delithiation proceeds, which greatly complicates the modeling work. Mathematical models for equilibrium potentials of Li-ion battery electrodes can be found in literature. For the main electrode materials like LiCoO₂, LiMn₂O₄, LiFePO₄, analytic expressions¹⁸ and theoretical thermodynamic laws^19,20 have been proposed. Hereafter, detailed information concerning charge and discharge equilibrium potentials of LiFePO₄ and graphite materials are given.

During lithiation / delithiation, lattice misfits between different phases are responsible of accommodation energy that can lead to hysteresis phenomenon. It was demonstrated that LiFePO₄ electrodes²¹ and graphite electrodes^22,23 show pronounced Open Circuit Voltage (OCV) hysteresis phenomena. Detailed thermodynamic explanations of this phenomenon can be found in Ref. 10 for the LiFePO₄ electrode. Dreyer et al.²¹ reported a permanent hysteresis of approximately 20 mV for LiFePO₄. In order to properly model hysteresis phenomena on both electrodes, specific experimental tests were performed on half-cells in this work. The hysteresis phenomena are represented in Fig. 2 for LiFePO₄ and graphite materials.

Zoom In Zoom Out Reset image size

In the present work, dynamic hysteresis is considered. As the diffusion phenomena in the solid particles have already been considered, the equations reported by Roscher and Sauer²⁴ are adapted to our model as follows:

Eq. 52 gives the evolution of the hysteresis factor Γ as a function of the current and of the hysteresis parameter χ, which is determined experimentally and is used to quantify the dynamics of the hysteresis.²⁴ The state value of Γ is limited to the range between 0 and 1, thus Eq. 52 is used within this restricted range. As can be seen in Eq. 49, when the values of Γ are 1 and 0, OCV_bat is respectively equal to OCV_ch and OCV_dch. From the thermodynamic equilibrium potentials of each electrode in charge and discharge, and based on the experimentally determined stoichiometries θ defined in Eq. 46, the open circuit voltage of the full cell, OCV_bat, can be reconstructed (Eqs. 46, 50–51) for both charge and discharge, as represented in Fig. 3.

Zoom In Zoom Out Reset image size

The thermal behavior of Li-ion cells has been widely investigated in order to optimize thermal management and to address potential security issues. In this work, a simple lumped parameter energy balance was considered as described in Table II. Both irreversible and reversible (entropic) contributions of thermal power generated are considered, as shown in Eq. 15. Thermal properties of electrode materials are also well reported in literature.^25,26 The entropic contributions presented in Eq. 15 have been investigated in Ref. 26 for LiFePO₄ and Ref. 27 for graphite materials. As reported by Reynier et al.,²⁷ one can observe an hysteresis phenomenon on the entropic heat for graphite electrodes. Entropic contributions can be important especially for low current rates, as shown in Ref. 28 and in the model validation part below. As temperature control is a key point in BMS design to ensure security and long life of the system, an estimation of the internal temperature of the cylindrical cell has to be introduced in the present model. Many authors have investigated the spatial thermal behavior of cylindrical Li-ion battery cells.^28,29 Rigorous investigation of the thermal behavior should be carried out through multi-dimensional modeling. However, to limit computing time, the simple thermal model developed by Forgez et al.³⁰ to estimate the internal temperature of the cylindrical system was used in this work. The following equation links the internal temperature, T_int, the skin temperature, T_skin, and the ambient temperature, T_amb:

The thermal resistance, R_th,in, modeling the heat transfer inside the cell, can be defined as a function of the thermal conductivity, λ_cell, and radius, r_cell, of the cell:

On the other hand, the thermal resistance, R_th,out, modeling the heat transfer between the cell and the environment can be expressed as a function of the Newton convective coefficient, h_conv, and the lateral surface of the cylindrical cell A_cell:

From Eqs. 53–55, the internal temperature can be expressed as:

In first approximation, Eq. 56 provides quantitative information about the impact of the cooling system (air cooling, liquid cooling) on the thermal gradients inside the cell. The values of the parameters of the simplified electrochemical and thermal model are presented in the following section.

Experimental work was performed on commercial LiFePO₄-graphite cells ANR26650 M1 (2.3 Ah) from A123 Systems. C/10, C/4, C/2, 1C, 2C, 4C, and 8C charges and discharges were performed in a climatic chamber at 0°C, 10°C, 23°C, and 33°C using a multipotentiostat (VMP, Biologic, Claix, France) outfitted with 20 A boosters for high C-rate measurements. For discharge experiments, the initial state of the cells was obtained by charging at a constant current (CC) of 1C until 3.6 V, then a constant voltage (CV) of 3.6 V was maintained until the current was lower than C/20, finishing with a rest period of 4 h (1C-CC-3.6V-CV until I < C/20, rest 4h). Discharges at different rates were then carried out from the same initial state. For charge experiments, the initial state of the cells was obtained by discharging at a constant current of 1C until a cell potential of 2 V, then a constant voltage of 2 V was maintained until the current was lower than C/20, finishing with a rest period of 4 h (1D-CC-2V-CV until I < C/20, rest 4h). Charges at different rates were then carried out from the same initial state.

The experimental characterization of the inherent properties of the materials under study is an integral part of the modeling effort.³¹ Indeed, parameterization of physics-based models requires multi-level and multi-physics experimental devices, even if some properties of the materials can be found in literature. Electrochemical characterization tests were conducted at both cell and electrode levels. In order to perform tests at the electrode level, some cells were operated in an argon-filled glove box after complete discharge at 2 V. Scanning electron microscopy was performed on both LiFePO₄ and graphite electrodes to determine some of the design parameters. Pieces of the recovered electrodes were then separately reassembled in coin cells with a Li-foil counter electrode to determine the open circuit voltage and the electrochemical window for each electrode. For the LiFePO₄-graphite system under study, both literature and measured data were collected and reported in Table III.

Without considering aging, the charge balance on the cyclable Li can be checked thanks to the following equations:

Using the data reported in Table III, Eqs. 57 and 58 give a value of 2.32 Ah for the charge capacity of the electrodes, Q_p and Q_n, which is in good agreement with the nominal capacity of 2.3 Ah given by the cell manufacturer.

Special electrochemical features of LiFePO₄-graphite systems such as (dis)charge path-dependence have been investigated and modeled, incorporating particle size distribution and concentration-dependent solid state diffusion coefficient.^9,12,13 Moreover, a thermodynamically-consistent correction of the activity could be introduced to refine the present model on the basis of the experimental open-circuit voltages presented in Fig. 2.²⁰ However, the diffusion coefficients in Table III were set as constants to keep the model simple enough for both the positive and negative electrodes.

To precisely model the diffusion overpotential in the liquid phase (Eq. 38), the electrolyte conductivity has to be determined accurately. In this work, due to a lack of data on the electrolyte composition of the LiFePO₄-graphite cell investigated, a simple chemical analysis was carried out to determine the composition of the solvent mixture. This analysis indicates ethylmethyl carbonate (EMC) and dimethyl carbonate (DMC) as the main components (ratio 1:1 (w:w)) of the mixture. Unfortunately, the conductivity dependence on Li concentration of this specific composition was not found in literature. However, this electrolyte could be compared to other electrolyte compositions reported in Table IV.

The electrolyte conductivity varies as a function of the solvent mixture composition, as shown in Fig. 4. The conductivity curves generally present a maximum around a Li concentration of 1 mol·dm⁻³. The conductivity models in the literature report values of the solution conductivity κ. The effective conductivity of the electrolyte in the porous media, κ^eff, can be determined from the Bruggman correction (Eqs. 40–42). The Bruggman exponent can be used to quantify the impact of the porous media on the electrolyte conductivity.³³ Indeed, as demonstrated by Doyle for the EC/DMC electrolyte mixture, the variation of the Bruggman parameter from 1.5 to 4.5 can strongly change the effective conductivity.³³ The higher the Bruggman coefficient, the lower the effective conductivity. Recently, Thorat et al. quantified the tortuosity in porous Li-ion electrodes and reported different values for the Bruggman exponents ranging from 1 to 3.3.³⁴

Zoom In Zoom Out Reset image size

In order to select the most appropriate electrolyte conductivity model from those reported in Table IV, electrochemical impedance spectroscopy (EIS) results performed at 25°C and a SOC of 50% were used. The high-frequency intercept with the real axis in the Nyquist representation of the impedance diagram was used to determine the approximated value of the ohmic resistance. From Eq. 39, and using the design parameters in Table III, the static resistive contributions of the electrolyte in the negative and positive electrodes and in the separator were calculated as a function of the Li concentration. Assuming a Li concentration of 1.2 mol.dm⁻³ and classical Bruggman exponents of 1.5 due to a lack of available data, a good agreement between the experimental EIS value and the predicted ohmic resistance with the electrolyte model "c" is obtained, as can be seen in Fig. 5.

Zoom In Zoom Out Reset image size

For the thermal model calibration, literature data on thermal parameters were collected and experiments were performed in an adiabatic calorimeter to estimate missing values (see Table V). Temperature dependency is rather complex for the electrolyte. In the present work, temperature dependencies were modeled by a simple Arrhenius law. This assumption can strongly limit the temperature range availability of the model especially under 0°C.

In this part, as a unique particle size is considered, the capability of the modeling concept proposed by Delacourt et al.¹⁴ for the understanding of lithium insertion and deinsertion in Li_yFePO₄ with a single particle approach is investigated. In Ref. 14, the authors fitted experimental charge and discharge curves of LiFePO₄ electrodes using a model with a current-dependent particle radius. In this study, the correlation between the solid state diffusion coefficient, the particle radius and the surface resistance was demonstrated. Mentioning a mosaic model to explain the uncommon dependency of the particle radius as a function of the current regime, the asymmetry between charge and discharge was shown for LiFePO₄. As the mathematical structure of the present model is similar to that of the SP model used in Ref. 14, this concept was tested within the framework of our simplified electrochemical and thermal model. This approach was adopted to represent the capacity restitution in charge and discharge. In this preliminary modeling work, the diffusion coefficient in the solid phase was taken as temperature-dependent only and did not vary with concentration. Charge and discharge experiments at C/4, C/2, 1C, 2C, 4C and 8C were performed at 10°C, 23°C and 33°C. Figs. 6 and 7 present the simulation results compared to the experimental data at 23°C.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As can be observed in Figs. 6 and 7, deviations between experimental data and model predictions for both cell voltage and skin temperature are under 2% for charge and discharge rates up to 2C over the full SOC range. Larger deviations can be seen for 4C and 8C rates. At a 8C rate the test in charge ended prematurely since the upper cutoff voltage of 3.6 V was reached. To obtain these results, the radius of the spherical particles at the positive and negative electrodes were chosen as adjustment parameters as done in Ref. 14. Indeed, the behavior of the negative electrode is extremely complex as insertion and deinsertion of Li ions occur.^20,23 The evolution of the fitted particle radii upon current is shown in Fig. 8. As can be observed in Fig. 8, the higher the current intensity in charge or discharge, the smaller the particle radius. This result could be explained by the fact that smaller particles are filled faster than larger particles. The deinsertion/insertion processes in smaller particles could control the overpotentials of the electrodes especially if it is considered that the extreme SOCs are reached quicker for smaller particles. The mosaic model concept may also explain this phenomenon. In this concept, as described in detail in Ref. 14, boundaries between Li-rich and Li-poor phases could be formed preferentially in place of the growth of existing domains at large current density. The nucleation of multiple phase boundaries could delimit smaller diffusion domains and thus increase the apparent electroactive surface area, hence leading to an apparent reduced particle radius.

Zoom In Zoom Out Reset image size

The present simulation demonstrates that the implementation of a current-dependent radius of the single particles in the simplified electrochemical and thermal model gives good results with respect to capacity restitution from C/4 until 8C rates. These results are in good accordance with those reported in Ref. 14 for LiFePO₄. Interestingly, an asymmetry between charge and discharge is observed for both positive and negative materials. Therefore, the model can be considered as validated within the calibration range of C/4 - 8C rates.