{'fraction_non_alphanumeric': 0.043786562666193934, 'fraction_numerical': 0.01517904626839213, 'mean_word_length': 4.653846153846154, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Core-shell enhanced single particle model for lithium iron phosphate batteries: model formulation and analysis of numerical solutions

In this paper, a core-shell enhanced single particle model for iron-phosphate battery cells is formulated, implemented, and verified. Starting from the description of the positive and negative electrodes charge and mass transport dynamics, the positive electrode intercalation and deintercalation phenomena and associated phase transitions are described with the core-shell modeling paradigm. Assuming two phases are formed in the positive electrode, one rich and one poor in lithium, a core-shrinking problem is formulated and the phase transition is modeled through a shell phase that covers the core one. A careful discretization of the coupled partial differential equations is proposed and used to convert the model into a system of ordinary differential equations. To ensure robust and accurate numerical solutions of the governing equations, a sensitivity analysis of numerical solutions is performed and the best setting, in terms of solver tolerances, solid phase concentration discretization points, and input current sampling time, is determined in a newly developed probabilistic framework. Finally, unknown model parameters are identified at different C-rate scenarios and the model is verified against experimental data.


Introduction
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 (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 LiFePO 4 have made it an attractive positive electrode material for LIB. In particular, the highly reversible transition between the olivine lithiated (LiFePO 4 ) and delithiated (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.
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 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 twophase transition dynamics during charge and discharge in LiFePO 4 /graphite batteries. By means of the core-shell modeling principle, both the one-phase at the beginning (Li-rich) and end (Lipoor) 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 LiFePO 4 /graphite batteries. In Section 3, the coreshell 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.

Physical principles
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 4 -graphite system, the half-cell reactions at the negative and positive electrode, are, respectively: where a and b are coefficients defined between 0 and 1. According to [3], during the LiFePO 4 /FePO 4 phase transition most of the positive electrode reaction is dominated by the coexistence of two phases: • α-phase: Li-poor phase Li α FePO 4 , with α = 1 − b and α 0; • β-phase: Li-rich phase Li β FePO 4 , with β = 1 − b and β 1.
α-phase and β-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 and c β s,p are the lithium concentrations for α-phase and β-phase (assumed uniform in the electrode), and c max s,p 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 ch p and U dis p 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 α-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 θ α p is reached, the formation of the β-phase starts. Here it is when the α-phase starts transitioning into the β-phase and the coexistence of two phases leads to a constant OCP (Figure 3(b)). The transition ends at the normalized concentration θ β p . After this point, the positive particle is all in its β-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 β-phase (Figure 4(c)) to α-phase (Figure 4(a)). As shown in [7], at the start of the charging process the surface of the positive electrode is in β-phase, with a corresponding different OCP. Ultimately, the presence of two OCPs for the positive electrode is related to the voltage path dependence behavior.
3 Cell governing equations 3.1 Model development As shown in Figure 1(b), the negative electrode is approximated as a single spherical particle and following [15], electrochemical dynamics are described by means of mass transport in the electrolyte and solid phase (Equations (64) and (66)), and charge transport in the electrolyte phase (Equation (65)). Mass and charge transport equations are used to model the electrolyte dynamics in the separator (Equations (67) and (68)).
As shown in Figure 1(b), the positive electrode is modeled as one single particle with two phases: Phase #1 (the core) and Phase #2 (the shell). Conventional mass and charge transport equations are applied to the electrolyte phase (Equations (69) and (70)). To describe the presence of two phases, the solid phase mass transport equation is enhanced to account for both the one-phase and two-phase regions. As described in Section 2, the positive electrode behavior is a function of charge and discharge conditions, that are characterized by two different OCPs: U ch p and U dis p . During discharge, the U dis p OCP is used and the positive particle is α-phase. The OCP decreases until the normalized concentration reaches θ α p and, after this point, the formation of the β-phase starts (Figure 3(b) and (e)). Two-phase coexist inside the particle whose behaviour is captured by the core-shell paradigm, where the core and shell of the particle are at α-phase and β-phase, respectively. The core is at a constant and uniform concentration c α s,p = θ α p · c max s,p and subject to a shrinking process, as discharge takes place, which replaces the α-phase with β-phase. This phenomenon is modeled by means of the following mass balance equation: Equation (4) describes the motion of the interface (or boundary) r p between the α-phase and βphase while assuming dr p /dt to be function of the concentration gradient ∂c s,p /∂r only. In Equation (4), c α s,p and c β s,p are constants defined as a function of θ α p and θ β p , respectively (see Equation (3)), c s,p is the solid phase concentration, D s,p is the constant solid phase diffusion coefficient, and r is the radial coordinate (i.e., the coordinate along the radius of the particle). The term sign(I) accounts for the fact that, during discharge, α-phase transitions to β-phase and that the opposite happens during charge. In Figure 3(e), the moving boundary r p is the distance between the center of the particle and the interface between α-phase and β-phase. The model is completed by the following boundary condition: c s,p r=rp = c β s,p enforcing the concentration at the interface r p to be always equal to c β s,p . Assuming the transition from Figure 3(d) to (e) to happen at the time instantt, the following initial conditions are introduced: Equation (6) is the initial condition for the moving boundary r p when, at time instantt, the system evolves from one-phase to two-phases. Att, the boundary r p is initialized at R p − , with a small enough coefficient (< 0.001R p ). Equation (7) defines the initial condition for core and shell solid phase concentrations at the transition time instantt. The core is at constant and uniform concentration c α s,p and, as the core shrinks, c α s,p is replaced by the β-phase concentration c β s,p . In the shell region (i.e., for r ≥ r p ), the solid phase concentration starts at c β s,p and rises according to the governing equation (71). Once the core is completely depleted, the particle is fully in β-phase and lithium intercalates until the end of the discharge process, as shown in Figure 3(f).
During charge, the opposite process, described in Figure 4, occurs. The core is in β-phase and shrinks until the whole particle is in α-phase. General equations for charge and discharge are provided in Equations (73), (75), and (76).
Governing equations are summarized in Table 1. Additional equations for the core-shell ESPMtransport parameters, active area, porosity, cell voltage, overpotential, and SOC -are summarized in Tables 2 and 3.
The transition from the one-phase to two-phase is performed when θ bulk p = θ * p , i.e., when the following mass balance is satisfied: with c s,p the one-phase concentration at the transition time instantt. In Equation (9), the normalized core concentration θ * p is equal to θ α p or θ β p depending on whether the battery is being discharged (Figure 3) or charged ( Figure 4). Similarly, in Equation (9), c * s,p is equal to c α s,p and c β s,p for discharge and charge, respectively. The transition from two-phase (Figures 3(e) and 4(e)) to one-phase (Figures 3(f) and 4(d)) is already taken into account by the core-shell formulation, i.e., the transition happens when the core is completely depleted and r p reaches zero.

Numerical solution
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].

Coordinate system transformation: normalization
Starting from the positive electrode core-shell model in Section 3, the transformation proposed in [10] is used to move from the radial coordinate system (r) to the normalized coordinate χ: where r represents a radial position in the particle and r p is the moving boundary. This transformation remaps the discretization of the shell region from [r p , R p ] to [0, 1], making the domain stationary while the boundary is moving. In Equation (6), the small enough avoids the rise of singularities during transitions from one-phase to two-phase. As proposed in [13], the left hand side of the positive particle diffusion Equation (71) is transformed from r to χ domain as follows: Introducing the following relationships: Equation (11) can be rewritten as: The change of coordinates applied to the right hand side of Equation (71) is obtained using the following transformations: Given Equation (14), the following reformulation is obtained: 2D s,p r ∂c s,p ∂r = ∂c s,p ∂χ Given Equations (13), (15), and (16), the solid phase mass transport in the positive particle (Equation (71)) is rewritten as: Starting from Equation (14), boundary conditions and mass balance in Equation (73) are rewritten in terms of the χ coordinate: where χ = 0 coincides with the moving boundary r p , χ = 1 corresponds to the surface or R p , and g(I) is the concentration at the moving boundary (as defined in Equation (75) of Table 2).

Discretization
The core-shell framework -Equations (17), (18), (19), and (20) -is discretized into N r,p points ( Figure 5). The right-sided and central finite difference schemes are used for the first and second derivative approximations: where u is a generic variable and l defines the index of the discretization point χ l : In Equation (23), r l , the radial position along the shell region, takes the following form:

Discretization of mass balance
Equation (20) is rewritten as: Approximating the term ∂cs,p ∂χ 0 as: and given that c s,p 0 = g(I) (from Equation (19)), the discretized version of the mass balance is obtained:

Discretization of boundary conditions
At χ = 0, the moving boundary condition in Equation (19) takes the following form: and the fixed boundary condition (18) is rewritten as: Approximating the left hand side of Equation (29) as: the following is obtained:

Discretization of the solid phase mass transport
Equation (17) is discretized relying on both the right-sided and central finite difference schemes. For l ∈ [1, N r,p − 2], the following expression is obtained: with M 3 and M 4 defined as: At l = N r,p − 1, the discretized solid phase mass transport equation takes the following form: From Equation (31), Equation (34) is rewritten as: It is worth mentioning that, according to Equation (28), at l = 0 the time derivative ∂cs,p 0 ∂t = 0.

State-space representation
Equations (27), (28), (32), and (35) constitute a system of coupled ODEs which can be rewritten in state-space form. This is a practical reformulation for the implementation of the core-shell governing equations and their numerical solution.
Given the vector of discretized solid phase concentration states: and the moving boundary r p , the state vector x is defined as: Introducing the variables: the positive particle state-space representation is given by: In Equation (39), matrices A 1 and A 2 take the following expressions: Similarly, vectors B and G are defined as: Compared to a conventional ESPM state-space representation, [18], the state vector x in Equation (37) has the moving boundary r p as the new state from which the positive particle solid phase concentration depends on. Overall, the positive particle is described by N r,p states, among which N r,p − 1 are associated with the positive particle solid phase concentration and 1 with the moving boundary. This is in contrast with the conventional ESPM formulation, where the discretization leads to N r,p states describing the positive particle solid phase concentration dynamics. State-space representations for negative electrode solid phase and electrolyte mass transport dynamics, derived as shown in [18], are summarized in Table 4. The electrolyte concentration (c) is discretized along the Cartesian coordinate x with N x,n , N x,s , and N x,p discretization points for negative particle, separator, and positive particle, respectively. The negative particle solid phase concentration (c s,n ) is discretized along the radial coordinate r, with N r,n discretization points defined between the radial coordinates 0 and R n ( Figure 5).

Sensitivity analysis of numerical solutions
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 (abstol and 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 C = {N r , dt, reltol, 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.

Parameters identification
In this work, the values of region thicknesses (L p , L n , L s ), maximum solid phase concentrations (c max s,p , c max s,n ), porosity of electrodes and separator (ε p , ε n , ε s ), active volume fractions (ν p , ν n , ν 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 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 1 , the parameter vector Θ, 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 ∈ K = {charge, discharge}, C is the solver setting (determined according to Section 5), j is the time index, N is the number of samples, SOC k p and SOC k n are the simulated state of charge at the positive and negative electrodes, V k is the simulated voltage profile, V k exp and SOC k exp 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 Θ is the same for charge and discharge.
The identification problem is formulated as the following optimization problem: subject to (a) Governing equations (Table 4) The minimization of the cost function is subject to the system dynamics (a), and inequalities (b) and (c) to ensure the two-phase region (defined between θ α p and θ β p ) to be contained inside the positive electrode stoichiometric window θ p,0% -θ p,100% . Constraints (d) enforce the moving boundary r p to be greater than or equal to zero for j ∈ [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 ∈M = {p, n} indicates the positive and negative electrode. To enforce charge conservation, Q k i must satisfy the following constraint: where parameters Q and Q are suitable bounds defined as (1 + 0.01)Q k C/12 and (1 − 0.01)Q k C/12 , with Q k C/12 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.
6.1 Identification at C-rate C/12, C/10, C/6, and C/3 Charge and discharge experimental data at C/12 are used to identify the parameter vector: where R n , R p , and A cell are geometrical parameters, D s,n and D s,p are solid phase diffusion coefficients controlling the mass transport in the negative and positive electrodes, θ n,100% , θ n,0% , θ p,100% , and θ p,0% define the stoichiometric window of the cell, and R l is the lumped contact resistance. For the implementation of the core-shell paradigm, the identification of θ α p and θ β p is key to characterize the transition from one-phase (α or β) to two-phases (α and β).
At C/10, the following parameter vector is identified: Stoichiometric and geometric parameters are set as constants from the Θ C/12 identification results. This is reasonable since the proposed model does not consider electrode swelling phenomena [1] and the geometry is not changing while increasing the C-rate. The stoichiometric window -defined by θ i,0% and θ i,100% -describes the maximum and minimum concentration of lithium-ion in the positive and negative electrodes, and it does not change while increasing the C-rate because of mass conservation. On the other hand, the discharged (or charged) capacity at C/10 is lower than the discharged capacity at C/12. In this situation, the stoichiometric window is not changing, but the cell is using a narrow range inside the θ i,0% and θ i,100% window. Therefore, constraint (e) in the identification problem (45) is valid only for the C/12 scenario, where the full stoichiometric window is used, and deactivated for the other C-rates used in the identification. The overpotential contribution is small at C/12 and C/10, making the parameter k n and k p not identifiable, for the values of k n = k p = 10 −11 [m 2.5 /(mol 0.5 s)] are used. These values are initial guesses that, as done for C/6 and C/3, should be refined through identification. At C/6 and C/3, the following parameter vector is identified: In addition to the diffusion coefficients and lumped resistance -also identified at C/10 -since the overpotential contribution increases at higher C-rate, the reaction rate constants k n and k p are identified. PSO is a nonlinear optimization algorithm that searches the optimal parameter vector inside a space defined by the parameters' upper and lower bounds and it requires an initial guess for the parameter vector to be identified. In this paper, the identification is performed sequentially, starting with C/12 and C/10, C/6 and lastly C/3 to follow. At C/12, reasonable bounds and initial guesses are obtained from the available literature on LFP batteries [13], [20], [21] and information provided by our industrial sponsor (Table 6, second, third, and fourth columns). For C/10, C/6, and C/3, the PSO initial guesses are defined by the solution of the identification problem at C/12, C/10, and C/6, respectively. This methodology allows to identify a sequence of consistent parameters, avoiding falling in local minima far away from the C/12 solution. Provided physical consistency of the solutions, C/6 and C/3 scenarios are characterized by faster diffusion and increased overpotentials (Equations (87) and (88)) compared to C/12 and C/10. To incorporate this physical knowledge in the identification problem, PSO lower bounds for diffusion coefficients and reaction rate constants are modified by using the identified parameters from the previous identification. Moreover, contact losses R l (of the order of few mΩ for LFP batteries [23]) are expected to increase as C-rate increases.

Results
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 max s,p , 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 β (or α) 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 (c α s,p ) and transitions to c α s,p (c β s,p ), 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 SOCn and J SOCp below 0.005 (or, equivalently, 0.25 [Ah]). It is noted that for any given C-rate, J SOCp is always higher than J SOCn . 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 SOCp 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].

Conclusions
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.

Appendix A: relative and absolute tolerance
Relative (reltol) and absolute (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 (c s,p and c s,n ), electrolyte phase concentration (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 abstol = reltol × 0.001 (53) with the relative tolerances in the range [1 × 10 −9 , 1 × 10 −3 ]. A scaling factor of 0.001 of the abstol is the default approach used in Matlab [24]. Starting from the positive particle solid phase concentration (c s,p ), the product |y| × reltol in Equation (52) is computed with respect to the maximum and minimum admissible concentrations: Maximum concentration: |θ p,0% · c max s,p | × reltol Minimum concentration: |θ p,100% · c max s,p | × reltol where θ p,0% and θ 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 (c s,n ), with the product |y| × reltol defined as: Maximum concentration: |θ n,100% · c max s,n | × reltol Minimum concentration: |θ n,0% · c max s,n | × reltol where θ n,100% and θ 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 reltol are shown in Figure 10. A relative tolerance higher than 1 × 10 −5 leads to values of (54) and (55) higher than 1 [mol/m 3 ] (or 6.94 [g/m 3 ], given the molar mass of lithium equal to 6.94 [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| × reltol for the electrolyte phase concentration c (a) and the moving boundary r p (b). For the electrolyte phase concentration, |c 0 | × reltol, with c 0 = 1200 [mol/m 3 ] the initial electrolyte concentration from [21], is computed and compared to the absolute tolerance. For the moving boundary, |R p | × reltol and |R p | × reltol are calculated, with R p and R p the positive particle radius upper and lower bounds defined in Table 6 2 .
Ultimately, constraint (52) is satisfied for solid phase concentrations, electrolyte phase concentration, and moving boundary.

Appendix B: sensitivity analysis
The sensitivity analysis is divided in two steps.
Step 1 is used to select values for sampling time (dt) and absolute tolerance (abstol).
Step 2 is used to select the pair {N r , reltol} ensuring a robust solution of the model equations. For each solver setting C, the ODEs solution accuracy is defined with respect to the objective function in Equation (44), for the C/12 constant current discharge profile.
Step 1: Selection of abstol and dt. The sensitivity is performed with respect to all combinations of solver settings in Table 9, considering the initial guess Θ 1 for the parameter vector ( Table  6, column four).
For each setting 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 C, the solution is unstable and a nonlineardiverging 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. Figure 13 shows the cost function J(Θ 1 , C) for abstol = reltol × 0.001, considering all the combinations in Table 9. Each contour plot is associated with a different current profile sampling time dt. The yellow region highlights where the solver is not converging. To this region, we enforce a fictitious high cost equal to J(Θ 1 , C) = 100. Lack of convergence is shown for high N r (which leads to high order ODEs systems) and low relative tolerances (requiring higher solution accuracy before the stopping criterion (51) is met). As shown in Figure 13, a sampling time increase does not change the convergence region of the solver. However, a reduction of one order of magnitude of the sampling time (e.g., from 10 to 1 [s]) leads to a one order of magnitude increment of the simulation time (from 3 to 30 [s]). This is shown in Figure 14, where, for each sampling time dt and various absolute tolerance scaling factors, the average simulation time (t sim ) is computed over the N r -reltol region in which the solver is converging (i.e., blue regions in Figure 13).
For a given sampling time and absolute tolerance scaling factor, the optimal {N r , reltol } 1 setting is the one minimizing the cost function J(Θ 1 , C): In Figure 15, the optimal setting for dt = 50 [s] and abstol = reltol × 0.001 is shown. In this scenario, the cost function is minimized for: Step 2: Selection of reltol and N r . The same analysis proposed in Step 1 is repeated for 600 realizations Θ ι of the model parameter vector (with ι = 1, ..., 600). Each realization Θ ι takes the following expression: For each Θ ι , parameters are randomly picked inside bounds defined in Table 6 and all combinations of solver settings in Table 10 are tested. As a result of Step 1, absolute tolerance scaling and sampling time are fixed (Equation (59)). Similarly to Equation (57), for each realization Θ ι , the optimal pair {N r , reltol } ι is the one minimizing the cost function J(Θ ι , C): The optimum of problem (61) is a function of ι and, for each vector Θ ι , a different optimal pair {N r , reltol } ι is found. Our goal is to define a solver settingĈ suitable for most scenarios, i.e., parameter vectors Θ ι . To this aim, we rewrite the problem in a probabilistic framework and, upon a frequentist approach, we define the "Probability of {N r , reltol }" as the probability that a certain configuration {N r , reltol} is optimal.
The probability is defined as P({N r , reltol }) and takes the following expression: where N {N r ,reltol } is the number of times, given the set of tested vectors Θ ι , a solver setting {N r , reltol } is optimal and N is the total number of realizations Θ ι considered in the analysis, i.e., 600.
The probability distribution is shown in Figure 17(a), together with the average simulation time t sim in Figure 17(b). Most of the optimal solutions are in the high relative tolerance region, between 1 × 10 −6 and 1 × 10 −5 . This is reasonable since, for these tolerances, the problem becomes easier to solve and convergence of the solver is more likely. The dark blue region inside the red polygon in Figure 17(a) highlights the solver settings for which there are no optimal solutions due to lack of convergence of the solver. This happens for high values of N r and low relative tolerances. From Figure 17(b), it is seen that increasing the number of discretization points N r leads to an increase of the computational time as a higher order ODEs system is to be solved.
Summing the probabilities in Figure 17(a) along N r , while grouping by the relative tolerance, leads to the cumulative probabilities in Figure 18, from which it is found that 66.4% of the optimal solutions are at reltol = 1 × 10 −5 , which is chosen as the best relative tolerance setting. According to Figure 19, N r = 70 is associated with the best trade-off between probability and computational time -normalized with respect to the maximum simulation time at N r = 100 -and is chosen as the best setting for the radial coordinate discretization points.
Ultimately, Steps 1 and 2 define the best solver settingĈ for the ODEs system solution: Appendix C: identifiability of model parameters To assess the sensitivity of the model output with respect to small perturbations of the identified parameters, a sensitivity analysis is carried out considering the C/12 constant current scenario. Identified parameters in Table 6 (column five) are used as nominal condition and perturbed to compute the sensitivity matrix. The sensitivity of each parameter, computed as the Euclidean norm along the columns of the sensitivity matrix, is shown in Figure 20(a). Following the rationale of [15], starting from the sensitivity matrix, a correlation analysis is performed. Figure 20(b) allows to understand if perturbations of different parameters result in the same output response. If this happens, correlation coefficients assume high values (> 0.8 or < −0.8), meaning that parameters can not be identified uniquely. Generally, this correlation analysis is used to reduce the vector of identified parameters. However, the cell under investigation is still in the design phase and most of its parameters are unknown, hence, results in Figure 20 should be interpreted just as indicators of the goodness of the identification. Moreover, the correlation analysis is performed a posteriori and is a function of the nominal parameter vector and perturbation strategy. The analysis of other methodologies for practical and structural identifiability is left to future works [25]. Reference stoichiometry ratio at 0% SOC (i ∈M), Concentration at the boundary ic k Core initial condition sign Sign function α Positive particle α-phase β Positive particle β-phase n Negative particle p Positive particle s Separator avg Average along x bulk Particle bulk ch, dis Charge and discharge ef f Effective Both positive and negative electrodes are modeled with a single spherical particle (b). The positive particle is characterized by two phases: Phase #1 (the core) and Phase #2 (the shell). Cartesian (x) and radial (r) coordinates are shown along with the thicknesses L n , L s , and L p for negative particle, separator, and positive particle, respectively.   Figures (a), (b), and (c) show the behavior of the positive electrode OCP during discharge. During the one-phase regions (a) and (c), the potential decreases. As shown in (b), the coexistence of two phases leads to a flat OCP. Figures (d), (e), and (f) show the single particle used to describe the one-phase regions ((d) and (f)) and, in particular, the transition from the α-phase to β-phase (e). The transition from one-phase (d) to two-phase (e) at timet is shown.   One-phase regions are associated with r p /R p = 0 and a monotonic decrease (charge) or increase (discharge) of the normalized solid phase concentration c s,p /c max s,p . r p /R p > 0 is associated with the coexistence of two-phase and, therefore, to the core-shrinking process. In this region, during charge, the concentration moves from θ β p to θ α p (vice versa, during discharge, θ α p transitions to θ β p ). For both charge and discharge, the corresponding phases described in Figures 4 and 3 Figure 17: Probability of {N r , reltol } to be optimal (a) and average simulation time t sim (b). The dark blue region inside the red polygon in (a) highlights the settings that are never optimal due to lack of convergence of the solver.    = 0 Positive electrode (p)

Current convention
Concentration at the moving boundary Core-shell initial condition Diffusivity and conductivity Active area Porosity