Core-shell enhanced single particle model for lithium iron phosphate batteries: model formulation and analysis of numerical solutions

Electrification and the transition to clean and green energy and circular economy have found their solution in lithium-ion battery (LIB) technology. Lithium iron phosphate batteries (or LFP), among the first LIBs to be commercialized [1], are today standard in China, used mostly for electric scooters and small electric vehicles (EVs). LFP batteries use lithium iron phosphate ($\mathrm{LiFePO}_{4}$) as the cathode material alongside a graphite carbon electrode as the anode [1]. LFP batteries do not decompose at higher temperatures, thus providing thermal and chemical stability, which results in an intrinsically safer cathode material than other commercially available chemistries such as NMC and LCO batteries [2]. The superior safety and charge/discharge rate characteristics of $\mathrm{LiFePO}_{4}$ have made it an attractive positive electrode material for LIB. In particular, the highly reversible transition between the olivine lithiated ($\mathrm{LiFePO}_{4}$) and delithiated ($\mathrm{FePO_{4}}$) phases allows thousands of charge/discharge cycles with high rate capability at an ideal voltage plateau around 3.45 V vs. lithium metal.

Specifically, during charge (and discharge) three phases are observed on the OCV vs Amph-sec characteristic: 1) Li-rich phase (only $\mathrm{LiFePO}_{4}$), 2) two-phase transition where the $\mathrm{LiFePO}_{4}$/$\mathrm{FePO_{4}}$ phases coexist, and 3) Li-poor phase (only $\mathrm{FePO}_{4}$), [3, 4, 5].

Developing a model incorporating the description of these two phases is crucial to understand lithium intercalation and deintercalation in LFP batteries and gain information on the electrochemical states from which one can build $SOC$ estimators.

According to [6], conflicting models exist to describe the phase transition mechanisms, lithium insertion mechanisms, and the existence of a two-phase region. In [7], a many-particle model is proposed where lithium is exchanged between individual particles, and sequential lithiation and delithiation are demonstrated. Transport and kinetic equations are ignored, making the model unsuitable when a careful description of the electrochemical states is needed or at high C-rate. The authors of [8] propose a domino-cascade model to describe lithium intercalation and deintercalation in the positive electrode lattice. In this model, the phase transition is modeled with a front moving, without any energy barrier, inside the crystalline structure. Mosaic models have also been investigated. In this context, smaller nucleation sites within a larger single particle, each undergoing a phase change during charge and discharge, are considered [9].

To model the $\mathrm{LiFePO_{4}/FePO_{4}}$ phase transitions, in [10] a core-shell approach, similar to the one shown in [11] for the discharge of a metal hydride electrode, is proposed. At the cost of adding some complexity to the electrochemical model – i.e., a mass balance equation – this technique allows for a detailed description of the electrode lithiation and delithiation dynamics. Diffusion is assumed to be isotropic and phase transitions are described upon a shell and core phase interacting via a moving boundary. To describe the path dependence during the battery operation, the core-shell model prescribes different phases to the shell and core depending on whether the battery is charging or discharging. In [12] and [13], it is demonstrated the successful implementation of this core-shell model for the single particle model (SPM) and pseudo-two-dimensional (P2D), respectively.

In this work, a core-shell enhanced single particle model (ESPM) is formulated to model two-phase transition dynamics during charge and discharge in $\mathrm{LiFePO_{4}}$/graphite batteries. By means of the core-shell modeling principle, both the one-phase at the beginning (Li-rich) and end (Li-poor) of charge (or discharge) and the two-phase region are modeled. A careful discretization of the electrochemical partial differential equations (PDEs) for an effective and robust numerical solution is proposed and used to convert the model into a system of ordinary differential equations (ODEs). A sensitivity analysis with respect to input current sampling time, positive and negative particles radial discretization points, and solver tolerances is performed. Given the uncertainty of some model parameters, the sensitivity analysis is performed for different realizations of such parameters, and the best combination of solver tolerances, discretization points, and input current sampling time is determined in a newly introduced probabilistic framework. This allows to obtain a fast solution of the governing equations while guaranteeing accuracy and numerical stability. Moreover, an optimization problem is formulated to identify the unknown model parameters under different C-rate. Ultimately, the proposed approach, while accounting for both solid phase and electrolyte dynamics, allows to reduce the computational burden with respect to a P2D model [13].

The remainder of the paper is organized as follows. Section 2 describes the physical principles related to intercalation and deintercalation in $\mathrm{LiFePO_{4}}$/graphite batteries. In Section 3, the core-shell ESPM governing equations are formulated. In Section 4, core-shell model equations are discretized and converted into a system of ODEs. Moreover, a convenient state-space representation is provided that, when compared to the traditional ESPM state space formulation, contains a new state for the moving core-shell boundary. In Section 5, results of the sensitivity analysis of numerical solutions are shown (details can be found in Appendixes A and B). Section 6 describes the parameter identification procedure. Finally, in Section 7 the model performance over charge and discharge experimental data is shown.

In lithium-ion batteries, lithium moves from the positive to the negative electrode during charging ($I<0$) and from the negative to the positive electrode during discharging ($I>0$). During charging, lithium leaves the positive electrode (deintercalation) and enters the negative one (intercalation). Conversely, during discharging, the negative electrode experiences deintercalation as lithium intercalates the positive electrode. For a LiFePO₄-graphite system, the half-cell reactions at the negative and positive electrode, are, respectively:

where $\mathrm{a}$ and $\mathrm{b}$ are coefficients defined between 0 and 1. According to [3], during the $\mathrm{LiFePO_{4}/FePO_{4}}$ phase transition most of the positive electrode reaction is dominated by the coexistence of two phases:

• •

  $\alpha$-phase: Li-poor phase $\mathrm{Li_{\alpha}FePO_{4}}$, with $\alpha=1-\mathrm{b}$ and $\alpha\simeq 0$;

• •

  $\beta$-phase: Li-rich phase $\mathrm{Li_{\beta}FePO_{4}}$, with $\beta=1-\mathrm{b}$ and $\beta\simeq 1$.

$\alpha$-phase and $\beta$-phase are defined in terms of the positive particle normalized lithium concentration (i.e., the ratio between the actual and maximum lithium concentration):

where $c_{s,p}^{\alpha}$ and $c_{s,p}^{\beta}$ are the lithium concentrations for $\alpha$-phase and $\beta$-phase (assumed uniform in the electrode), and $c_{s,p}^{max}$ is the maximum concentration. A pictorial representation of the positive electrode with the two phases (Phase #1 and Phase #2) is shown in Figure 1(a). The meaning of Phase #1 and Phase #2 changes depending on whether the battery is being charged or discharged. Furthermore, two open circuit potentials (OCPs) for the positive electrode are defined for the charge and discharge operations. In this work, the OCPs are provided by the industrial partner of the project (see, Figure 2). The negative electrode OCP is obtained from the Galvanostatic Intermittent Titration Technique (GITT) test [14]. Positive electrode OCPs – $U_{p}^{ch}$ and $U_{p}^{dis}$ in Equation (85) – are obtained fitting charge and discharge experimental OCPs measured at low C-rate (C/50). Experiments for positive and negative electrode’s OCPs are performed on the half-cell with lithium metal electrode as a reference.

During discharge (see, Figure 3(a)), lithium intercalates into the positive electrode and the $\alpha$-phase is formed. This corresponds to a rapid decrease of the OCP. As discharge continues, lithium concentration increases in the positive electrode and, once the normalized concentration $\theta_{p}^{\alpha}$ is reached, the formation of the $\beta$-phase starts. Here it is when the $\alpha$-phase starts transitioning into the $\beta$-phase and the coexistence of two phases leads to a constant OCP (Figure 3(b)). The transition ends at the normalized concentration $\theta_{p}^{\beta}$. After this point, the positive particle is all in its $\beta$-phase and the OCP decreases until the end of discharge is reached (Figure 3(c)). As shown in Figure 4, during charge, the process is reversed and the positive particle transitions from $\beta$-phase (Figure 4(c)) to $\alpha$-phase (Figure 4(a)). As shown in [7], at the start of the charging process the surface of the positive electrode is in $\beta$-phase, with a corresponding different OCP. Ultimately, the presence of two OCPs for the positive electrode is related to the voltage path dependence behavior.

The system of coupled PDEs in Table 1 is discretized using the finite difference method (FDM), for mass transport in the solid phase, and finite volume method (FVM), for mass transport in the electrolyte phase. This allows to convert the governing equations of the core-shell ESPM model into a system of ODEs and algebraic equations. The system of ODEs can be rewritten into a convenient state-space representation and solved relying on the Matlab numerical solver ode15s [17]. This solver is effective for nonlinear/stiff problems and suitable for the solution of electrochemical transport equations.

In this work, we focus on the numerical solution of the positive electrode solid phase mass transport dynamics, where the core-shell process takes place. In the next sections, Equations (71) and (73) in Table 1 are discretized relying on the FDM [16]. The discretization of negative electrode solid phase and electrolyte mass transport dynamics is provided in Table 4. For additional details, readers are referred to [18].

Previous works from the authors [15, 18, 19] show that a number of discretization points $N_{x,n}=N_{x,s}=N_{x,n}=10$ and $N_{r,n}=N_{r,p}=10$ generates a stable and accurate solution of the ESPM governing equations. In the core-shell ESPM, though, the presence of the moving boundary $r_{p}$ as an additional state is a potential source of numerical instabilities and lack of convergence of the solver. A sensitivity analysis is therefore conducted in this work to find suitable solver settings for a robust solution of the ODEs system.

The sensitivity analysis is carried out with respect to the following parameters: the radial coordinate discretization points ($N_{r,n}$ and $N_{r,p}$), the input current sampling time ($dt$), and the solver absolute and relative tolerances ($\mathrm{abstol}$ and $\mathrm{reltol}$). Sampling time and solver tolerances control the solver convergence, whereas the radial coordinate discretization points control the accuracy of the solution of the solid phase diffusion (i.e., the number of ODEs). In this work, $N_{r,n}=N_{r,p}=N_{r}$ and results of the sensitivity analysis are shown as a function of $N_{r}$. On the other hand, for the solution of the electrolyte phase concentration dynamics along the Cartesian coordinate, we consider $N_{x,p}=N_{x,s}=N_{x,n}=10$, as these parameters do not affect the solver stability.

Different combinations of solver settings $\mathcal{C}=\{N_{r},dt,\mathrm{reltol},\mathrm{abstol}\}$ were tested and the optimal combination of values, chosen based on a probabilistic framework explained in details in Appendixes A and B, is given by:

The optimal solver settings in Equation (43) are used for all identifications and simulations shown in the next sections.

In this work, the values of region thicknesses ($L_{p}$, $L_{n}$, $L_{s}$), maximum solid phase concentrations ($c_{s,p}^{max}$, $c_{s,n}^{max}$), porosity of electrodes and separator ($\varepsilon_{p}$, $\varepsilon_{n}$, $\varepsilon_{s}$), active volume fractions ($\nu_{p}$, $\nu_{n}$, $\nu_{s}$), and transference number ($t_{+}$) were provided by the industrial partner of the project. Readers can refer to [21] to find similar values for maximum solid phase concentrations, transference number, porosities of separator and negative electrode, and active volume fractions of separator and negative electrode. Porosity and active volume fraction of the positive electrode, and region thicknesses can be separately identified, for example, following the methodology in [22].

Unknown model parameters are identified using voltage vs capacity data for a $\mathrm{LiFePO}_{4}$/graphite pouch cell charged and discharged at C/12, C/10, C/6, and C/3 constant current at 25^∘C. Charge and discharge capacities at the different C-rate are computed from experimental data and listed in Table 5. A second cell, with the same technical specifications and tested under the same conditions, is used to verify the goodness of the identified parameters. Using the particle swarm optimization (PSO) algorithm¹¹1Matlab particleswarm function: https://www.mathworks.com/help/gads/particleswarm.html, the parameter vector $\Theta$, described in the next section, is identified.

The identification is performed following the line of [15] and [19], i.e., minimizing the following cost function in constant current charge and discharge conditions:

where $k\in\mathcal{K}=\{\text{charge},\text{discharge}\}$, $\mathcal{C}$ is the solver setting (determined according to Section 5), $j$ is the time index, $N$ is the number of samples, $SOC_{p}^{k}$ and $SOC_{n}^{k}$ are the simulated state of charge at the positive and negative electrodes, $V^{k}$ is the simulated voltage profile, $V_{exp}^{k}$ and $SOC_{exp}^{k}$ are the experimental cell voltage and state of charge from Coulomb counting, respectively. The weights $w_{1}$, $w_{2}$, and $w_{3}$ are user-defined dimensionless parameters here equal to one. The parameter vector $\Theta$ is the same for charge and discharge.

The identification problem is formulated as the following optimization problem:

$\mathbf{subject\ to}$

The minimization of the cost function is subject to the system dynamics $\mathrm{(a)}$, and inequalities $\mathrm{(b)}$ and $\mathrm{(c)}$ to ensure the two-phase region (defined between $\theta_{p}^{\alpha}$ and $\theta_{p}^{\beta}$) to be contained inside the positive electrode stoichiometric window $\theta_{p,0\%}$-$\theta_{p,100\%}$. Constraints $\mathrm{(d)}$ enforce the moving boundary $r_{p}$ to be greater than or equal to zero for $j\in[1,N-1]$. In this work, constant current charge and discharge profiles are considered and the cell always reaches the one-phase at the end of the charge or discharge process. Therefore, the moving boundary $r_{p}$ reaches zero for $j=N$.

The four stoichiometric coefficients defining the stoichiometric windows in the positive and negative electrode, respectively, are identified once at C/12 while imposing charge conservation. Charge and discharge capacities are computed according to the following formula:

where $i\in\hat{\mathcal{M}}=\{p,n\}$ indicates the positive and negative electrode. To enforce charge conservation, $Q_{i}^{k}$ must satisfy the following constraint:

where parameters $\overline{Q}$ and $\underline{Q}$ are suitable bounds defined as $(1+0.01)Q_{C/12}^{k}$ and $(1-0.01)Q_{C/12}^{k}$, with $Q_{C/12}^{k}$ the C/12 charge and discharge capacity computed from experimental data (Table 5). These bounds are used to ensure the numerical feasibility of the identification problem.

The model parameters identified for each C-rate of operation are summarized in Table 6, whereas the model parameters that have been provided by the industrial partner of the project (region thicknesses, maximum solid phase concentrations, porosities, active volume fractions, and transference number) are kept fixed across the various identifications.

Identification results for C/12 charge and discharge are shown in Figure 6. The solution of the identification problem leads to values for $J_{V}$ of $0.0021$ and $0.0031$ for charge and discharge, respectively. Figure 6 shows the behavior of the positive electrode solid phase concentration $c_{s,p}$ normalized with respect to $c_{s,p}^{max}$, and moving boundary $r_{p}$ normalized with respect to the positive electrode radius $R_{p}$, where $r_{p}/R_{p}=0$ corresponds to one-phase regions. The solid phase concentration starts in the one-phase $\beta$ (or $\alpha$) region as the battery starts being charged (or discharged). As the battery enters the two-phase region, $r_{p}/R_{p}$ becomes greater than zero and the core-shrinking process takes place, i.e., the core shrinks until the one-phase region at the end of charge or discharge is reached. During charge (discharge), the core is at uniform concentration $c_{s,p}^{\beta}$ ($c_{s,p}^{\alpha}$) and transitions to $c_{s,p}^{\alpha}$ ($c_{s,p}^{\beta}$), as the core shrinks.

In Figure 7, identification results for C/12, C/10, C/6, and C/3 are compared and corresponding values of the cost function are summarized in Table 7. For positive and negative particles, the model replicates the $SOC$ from Coulomb counting, with $J_{SOC_{n}}$ and $J_{SOC_{p}}$ below 0.005 (or, equivalently, $0.25\$[Ah]). It is noted that for any given C-rate, $J_{SOC_{p}}$ is always higher than $J_{SOC_{n}}$. This trend is attributed to the presence of the core-shell dynamics, which makes the numerical solution of the positive particle solid phase concentration and the minimization of the cost term $J_{SOC_{p}}$ challenging. Performances of the model with respect to the cell voltage profile are satisfactory for C/12, C/10, and C/6, with the cost function $J_{V}$ lower than or equal to 0.0054 (or $15\$[mV]). As shown by the increasing value of $J$ (last row in Table 7), the discrepancy between model outputs and experimental data increases while increasing the C-rate. As shown in Figure 8, the staging from the graphite OCP becomes more pronounced at C/3, which determines the discrepancy between experimental and simulated output voltage profiles. This prevents the proper convergence of PSO and leads to $J_{V}$ equal to 0.0215 (or $60\$[mV]). This behavior is due to limitations of the proposed core-shell ESPM model and will be further analyzed in future works.

The model performance is verified against experimental data acquired for a second LFP pouch cell tested at C/12, C/10, C/6, and C/3 constant current. Results are shown in Figure 9 and corresponding cost function values are summarized in Table 8. The performance of the model is satisfactory and trends are consistent with the identification results, i.e., increasing the C-rate the discrepancy between simulated and experimental data increases. In the worst case, $J_{V}$ reaches 0.0294, i.e., $80\$[mV].

This work presented the development of a core-sell enhanced single particle model for LFP batteries. Starting from the physical understanding of intercalation and deintercalation phenomena, governing equations are formulated and discretized using FDM and FVM within a core-shell modeling framework which assumes a growing shell, or lithium rich phase, and a shrinking core, or lithium poor phase, during discharge. The model parameters are identified via a dedicated optimization routine over different C-rate of charge and discharge. Finally, a state-space representation of the core-shell ESPM dynamics is provided, where it is shown that the positive electrode state vector depends on the moving boundary. A sensitivity analysis of the solid phase concentration discretization points, input current sampling time, and solver tolerances is performed to find a solver setting ensuring stable and robust numerical solutions.

The model developed in this paper provides a tool for the accurate description of LFP phase transition which can be used for state of charge estimation algorithm design. Future works will be devoted to model hysteresis and to the validation of the model under real-world driving conditions.

Relative ($\mathrm{reltol}$) and absolute ($\mathrm{abstol}$) tolerances represent relative and absolute error thresholds used to define the following convergence criterion for numerical solvers [17]:

with $e$ the estimated solver error and $y$ a generic variable. For practical purposes, the absolute tolerance should be “low enough” to be active only when the system is converging.

For the battery electrochemical model developed in this paper, the numerical convergence of the following variables is analyzed: solid phase concentrations ($\mathbf{c}_{s,p}$ and $\mathbf{c}_{s,n}$), electrolyte phase concentration ($\mathbf{c}$), and moving boundary ($r_{p}$). For a smooth convergence of the solver the following condition should be satisfied for all the four state variables:

In this appendix, the analysis is carried out considering the following scaling of the absolute tolerance

with the relative tolerances in the range $[1\times 10^{-9},1\times 10^{-3}]$. A scaling factor of 0.001 of the $\mathrm{abstol}$ is the default approach used in Matlab [24].

Starting from the positive particle solid phase concentration ($\mathbf{c}_{s,p}$), the product $|y|\times\mathrm{reltol}$ in Equation (52) is computed with respect to the maximum and minimum admissible concentrations:

where $\overline{\theta}_{p,0\%}$ and $\underline{\theta}_{p,100\%}$ are the upper and lower bounds of the corresponding stoichiometric coefficients, as shown in Table 6. The same procedure is followed for the negative particle solid phase concentration ($\mathbf{c}_{s,n}$), with the product $|y|\times\mathrm{reltol}$ defined as:

where $\overline{\theta}_{n,100\%}$ and $\underline{\theta}_{n,0\%}$ are the upper and lower bounds of the corresponding stoichiometric coefficients, as shown in Table 6. Results for Equations (54) and (55) as a function of $\mathrm{reltol}$ are shown in Figure 10. A relative tolerance higher than $1\times 10^{-5}$ leads to values of (54) and (55) higher than $1~{}[\mathrm{mol/m^{3}}]$ (or $6.94~{}[\mathrm{g/m^{3}}]$, given the molar mass of lithium equal to $6.94~{}[\mathrm{g/mol}]$). Hence, to ensure high accuracy solutions, the following constraint is enforced:

which is used in Tables 9 and 10, to define the relative tolerance upper bound.

Figure 11 shows the product $|y|\times\mathrm{reltol}$ for the electrolyte phase concentration $\mathbf{c}$ (a) and the moving boundary $r_{p}$ (b). For the electrolyte phase concentration, $|c_{0}|\times\mathrm{reltol}$, with $c_{0}=1200\ [\mathrm{mol/m^{3}}]$ the initial electrolyte concentration from [21], is computed and compared to the absolute tolerance. For the moving boundary, $|\overline{R}_{p}|\times\mathrm{reltol}$ and $|\underline{R}_{p}|\times\mathrm{reltol}$ are calculated, with $\overline{R}_{p}$ and $\underline{R}_{p}$ the positive particle radius upper and lower bounds defined in Table 6²²2$\overline{R}_{p}$ and $\underline{R}_{p}$ are scaled from meter to picometer..

Ultimately, constraint (52) is satisfied for solid phase concentrations, electrolyte phase concentration, and moving boundary.

The sensitivity analysis is divided in two steps. Step 1 is used to select values for sampling time ($dt$) and absolute tolerance ($\mathrm{abstol}$). Step 2 is used to select the pair $\{N_{r},\mathrm{reltol}\}$ ensuring a robust solution of the model equations. For each solver setting $\mathcal{C}$, the ODEs solution accuracy is defined with respect to the objective function in Equation (44), for the C/12 constant current discharge profile.

For each setting $\mathcal{C}$ in Table 9, the simulated voltage and $SOC$ for positive and negative particles are shown in Figure 12 and compared to C/12 experimental voltage and $SOC$ from Coulomb counting (black-dashed lines). For some combinations $\mathcal{C}$, the solution is unstable and a nonlinear-diverging $SOC$ profile is obtained (Figure 12(b)). This behavior, visible in the positive particle only, proves the challenges introduced by the core-shell dynamics and the need to perform the sensitivity analysis.