A physics-based fractional-order equivalent circuit model for time and frequency-domain applications in lithium-ion batteries

Journal of Energy Storage 64 (2023) 107150 Available online 22 March 2023 2352-152X/© 2023 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/bync-nd/4.0/). Contents lists available at ScienceDirect Journal of Energy Storage journal homepage: www.elsevier.com/locate/est Research papers A physics-based fractional-order equivalent circuit model for time and frequency-domain applications in lithium-ion batteries Pablo Rodríguez-Iturriagaa,∗, David Anseánb, Salvador Rodríguez-Bolívara, Manuela González b, Juan Carlos Vierab, Juan Antonio López-Villanuevaa aDepartment of Electronics and Computer Technology, Faculty of Sciences, University of Granada, Granada, 18071, Andalusia, Spain bDepartment of Electrical Engineering, Polytechnic School of Engineering, University of Oviedo, Gijon, 33204, Asturias, Spain ARTICLE INFO Keywords: Lithium-ion battery Equivalent circuit model Fractional-order model Physics-based model EIS ABSTRACT Equivalent circuit models (ECMs) remain the most popular choice for online applications in lithium-ion batteries because of their simpler parameterization and lower computational requirements in comparison to electrochemical models. Nevertheless, standard ECMs lack physical insight and fail to accurately reproduce cell behavior under a wide range of operating conditions. For this reason, the development of physics-informed ECMs becomes essential so as to provide a better description of the physical processes while maintaining a reduced computational complexity. In this article, we propose a novel physics-based ECM derived directly from an electrochemical model, so that there is a clear correlation between circuit states and internal battery states, as well as circuit and physical parameters. The proposed model yields an RMS error below 1.46 mV for cell voltage, 0.28% for the surface concentration in the active material particles, 0.6% for the electrode-averaged electrolyte concentration and 0.32 mV for the charge-transfer overpotentials. Another key feature of this model is the relationship between circuit parameters and those identified in frequency-domain tests, which allows us to characterize and validate the model experimentally. We understand that the presented model constitutes an alternative to standard ECMs as well as electrochemical models as it combines advantageous characteristics from both of them. 1. Introduction As concerns over energy supply and environmental issues grow larger worldwide, electrochemical energy storage has become a subject of intense research. In particular, rechargeable lithium-ion batteries have materialized as the leading storage solution for a number of applications [1–3], due to their high energy density, high specific energy, and low self-discharge [4]. Therefore, the development of comprehensive battery models is critical for their online monitoring by a Battery Management System (BMS) [5], thus ensuring their safe operation as well as a prolonged useful lifespan by accurately estimating internal battery states [6]. Three broad groups are commonly considered in the field of battery modeling: physics-based models, equivalent circuit models (ECMs) and data-driven models [4,7]. In data-driven approaches the battery is substituted by a black-box model that is able to reproduce its behavior. Nevertheless, the validity of the resulting model relies heavily on the training dataset. This may lead to practical drawbacks, such as a considerably time-consuming training process or overfitting issues [7]. Consequently, electrochemical and ∗Corresponding author. E-mail address: [email protected] (P. Rodríguez-Iturriaga). equivalent circuit models continue to be of great research interest, and we will focus on them in this study. Physics-based approaches model the physical processes and electrochemical reactions that occur in the cell, with the Doyle–Fuller– Newman (DFN) [8,9] being considered the standard battery model. Due to its complex nature, several approximations have been proposed so as to maintain the description of the physical processes with a restrained computational load [10], with the Single Particle Model (SPM) [11,12] being one of the most widely used. The SPM makes the assumptions that each electrode is composed of spherical active material particles with the same physical properties and that electrolyte dynamics are negligible [13], so the current distribution is uniform across the electrode. Consequently, only one representative particle from each electrode needs to be taken into account. This approximation holds true for low current rates [14]; however, at moderate and large currents the effects of concentration and potential gradients in the electrolyte should be accounted for. Therefore, several approaches have been proposed to include electrolyte dynamics with the goal of extending the range of applicability of the SPM [14–18], https://doi.org/10.1016/j.est.2023.107150 Received 22 December 2022; Received in revised form 16 February 2023; Accepted 14 March 2023 Journal of Energy Storage 64 (2023) 107150 2 P. Rodríguez-Iturriaga et al. usually by the names ‘enhanced SPM’ (eSPM) or ‘SPM with electrolyte (SPMe)’, via polynomial approximations of the electrolyte concentration profile [15,16,18] or the numerical resolution of the corresponding PDE [14,17,18], with the former approach being less suitable for pulsed currents [19]. There also exist some alternatives to the direct numerical resolution of the PDEs, in which the transfer functions of the system are obtained from electrochemical models and then transformed into a discrete state-space representation with the Discrete-Time Realization Algorithm (DTRA) [20] or into a SIMULINK [21] model via simplified fractional-order transfer functions [22–24]. Nevertheless, these models still require a precise set of physical parameters, whose determination from non-invasive experimental measurements continues to be an open research topic [25]. On the other hand, ECMs approximate the electrical behavior of a battery cell by that of a specific circuit. This approach is characterized by reduced memory requirements and a low computational load, therefore making them appropriate for online BMSs [26–29]. Moreover, their discrete state-space representation is easily obtained, thus allowing for the implementation of state observers and Kalman filters [30–32] to mitigate the influence of measurement and process noise. The challenge is twofold when developing ECMs: first, determining which circuit topology reproduces battery behavior more accurately, and secondly, obtaining the corresponding parameter values for said topology from experimental tests [33–35]. Most ECMs include a voltage source, which corresponds to the open-circuit voltage (OCV) as a nonlinear function of the cell state of charge (SOC), in series with the internal ohmic resistance and one or several parallel RC elements that model dynamic behavior [30]. A variation, useful for circuit simulators, substitutes the independent OCV source by a voltage source dependent on the state of charge of a capacitor which represents battery SOC [36]. However, these models introduce a strong separation between quasi-static electrode thermodynamics included in the OCV-SOC relationship [37] and the dynamic behavior represented by the RC network, in which different diffusion processes are mixed and a direct equivalence with the internal states of the battery is not often possible. Additionally, physics-based models show that electrode potential is dependent on the lithium concentration at the surface of active material particles rather than their average concentration, which is what SOC stands for. The difference between both magnitudes may be considerably large at high current rates, as predicted by the SPM [38]. Li et al. [39] recently proposed a modification to the model by Chen and Rincon-Mora [36] by which the non-linear voltage source depends on an intermediate variable named 𝑆𝑂𝐶𝑠𝑢𝑟𝑓 , which is obtained as the result of the series connection of the SOC capacitor and a RC-element that represents the non-uniform concentration profile due to lithium diffusion in the active material particles. Nevertheless, the relationship between circuit parameters and their physical counterparts is generally lost, which makes these models inaccurate outside the range in which they have been determined experimentally. For these reasons, several works combine an equivalent circuit model for electrolyte dynamics and ohmic losses with the transfer functions for solid diffusion in physics-informed reduced order models [40–43]. These works make use of Padé’s method in order to obtain rational approximations to these transcendental transfer functions. However, it has to be pointed out that Padé’s approximants are calculated in the neighborhood of a certain point, typically the frequency 𝑠= 0, so the accuracy of the approximation is not guaranteed within the whole frequency range. Furthermore, it has been shown that diffusion processes are more accurately modeled by a continuous distribution of time constants [44], therefore suggesting that fractionalorder transfer functions are more appropriate in this context. This is also corroborated by Electrochemical Impedance Spectroscopy (EIS) tests [6,44–47]. The Nyquist plot of EIS data shows several depressed semicircles and constant-slope tails that may be modeled electrically by Constant-Phase Elements (CPEs) [48] and ZARC elements [49], which are the parallel connection of a resistor and a CPE. The introduction of fractional-order circuit elements entails the computation of fractional-order derivatives, so several approaches have been proposed in the literature for approximating their behavior in the time domain [49–53]. Moreover, fractional-order ECMs have been employed successfully in conjunction with Kalman filters in recent years [54–56] for the concurrent estimation of battery state of charge and electrical parameters. For all the reasons above, the development of a physics-based ECM, which is able to provide information about internal battery states and whose circuit parameters are directly correlated with physical parameters, is an open research topic. Furthermore, we also consider the relationship between the time-domain and frequency-domain behaviors of lithium-ion batteries to be a subject of interest. Therefore, in this article we derive a novel reduced fractional-order ECM from an electrochemical model, whose states and parameters are related to those of the battery cell. Additionally, we have analyzed the correspondence between said circuit parameters and those identified from EIS measurements, thus clarifying the connection between frequency-domain data and physical parameters. For this purpose, we have obtained simplified transfer functions from the SPMe and determined their electrical equivalent with fractional-order elements. Subsequently, the proposed model has been validated first against the SPMe for given a parameter set, and next with experimental data by identifying equivalent circuit parameters from EIS measurements. The major advantage of this model over previous ECMs is the extended insight into the internal battery states and physical processes, whereas it groups physical parameters in resistors and time constants and presents lower computational requirements with respect to the SPMe, thus allowing for a simpler parameterization process as well as the implementation of online estimation algorithms. The main contributions of this article are condensed as follows: 1. Deriving reduced transfer functions directly from the SPMe and establishing their electrical equivalent via ZARC elements, as well as the criteria for the validity of the employed approximations. 2. Constructing an equivalent circuit model whose parameter values are directly correlated with their physical counterparts, and whose states contain information about the internal states of the battery cell. 3. Providing a method to determine said equivalent circuit parameters from EIS data, thus bridging the gap between the time and frequency-domain behaviors of lithium-ion batteries. In consequence, this paper is structured as follows: the analytical derivation of the proposed equivalent circuit model is presented in Section 2. The theoretical validation against an electrochemical model is carried out in Section 3, whereas the experimental parameterization process and results are detailed and discussed in Section 4. Some final remarks are provided in the last section. 2. Equivalent circuit development In this section, the procedure to obtain simplified transfer functions from the SPMe and their electrical equivalence is described with the goal of developing an equivalent circuit model. The main advantage of employing the SPMe instead of the DFN as in [20] lies in the fact that the former removes the coupling between the partial differential equations. Therefore, each diffusion process may be analyzed independently, thus producing considerably simpler transfer functions and allowing for a direct electrical equivalence. Journal of Energy Storage 64 (2023) 107150 3 P. Rodríguez-Iturriaga et al. 2.1. SPMe model description Marquis et al. [14] proposed an SPMe derived as an asymptotic reduction of the standard DFN model, which expresses battery terminal voltage as a function of electrode-averaged quantities. In said model, the output voltage is split into the sum of its components as in Eq. (1): 𝑉=𝑈𝑒𝑞 +𝜂𝑐+𝜂𝑟+𝛥𝛷𝑒+𝛥𝛷𝑠(1) where the term 𝑈𝑒𝑞 stands for the equilibrium potential, 𝜂𝑐and 𝜂𝑟represent the voltage drop due to concentration gradients in the electrolyte and charge-transfer reactions, respectively, and 𝛥𝛷𝑒and 𝛥𝛷𝑠are the ohmic losses in the electrolyte and solid, respectively. The detailed expression for each term is shown in Eq. (2): 𝑈𝑒𝑞 =𝑂𝐶𝑃𝑝(𝜒𝑝|𝑟=𝑅𝑝) − 𝑂𝐶𝑃𝑛(𝜒𝑛|𝑟=𝑅𝑛)(2a) 𝜂𝑐=2𝑅𝑇 𝐹 (1 − 𝑡+) 𝑐𝑒,𝑡𝑦𝑝 (𝑐𝑒,𝑝 −𝑐𝑒,𝑛)(2b) 𝜂𝑟=2𝑅𝑇 𝐹[sinh−1 (𝐼𝑅𝑝 6𝑗0,𝑝𝜖𝑝𝐿𝑝𝐴)+ sinh−1 (𝐼𝑅𝑛 6𝑗0,𝑛𝜖𝑛𝐿𝑛𝐴)] (2c) 𝛥𝛷𝑒=𝐼 𝜅𝐴 (𝐿𝑛 3𝜖𝑏 𝑒,𝑛 +𝐿𝑠 𝜖𝑏 𝑒,𝑠 +𝐿𝑝 3𝜖𝑏 𝑒,𝑝 )(2d) 𝛥𝛷𝑠=𝐼 3𝐴(𝐿𝑛 𝜎𝑛 +𝐿𝑝 𝜎𝑝)(2e) where 𝑂𝐶𝑃𝑝and 𝑂𝐶𝑃𝑛are the open-circuit potentials for the positive and negative electrodes, respectively. 𝜒𝑝|𝑟=𝑅𝑝and 𝜒𝑛|𝑟=𝑅𝑛are the normalized lithium concentrations in the positive and negative active material particles, respectively, evaluated at their surface, with 𝑅𝑝and 𝑅𝑛being their respective particle radius. Regarding the electrolyte, 𝑡+is the cation transference number, 𝜅represents the electrode conductivity and 𝑐𝑒,𝑝 and 𝑐𝑒,𝑛 are the electrode-averaged electrolyte concentrations in the positive and negative electrode, respectively. 𝑗0,𝑝,𝑗0,𝑛,𝜎𝑝,𝜎𝑛𝜖𝑝 and 𝜖𝑛are the exchange current densities, conductivities and active material volume fractions for both electrodes. 𝐿𝑛,𝐿𝑠,𝐿𝑝,𝜖𝑒,𝑝,𝜖𝑒,𝑠 and 𝜖𝑒,𝑛 represent the thickness and porosities of the positive electrode, separator and negative electrode, respectively. Lastly, 𝐴stands for the electrode area, 𝑏is the Bruggeman coefficient and 𝐼is the current applied to the cell. Note that we have changed the current sign criterion to positive while charging in order to obtain impedance expressions with a positive real part. The terms that depend directly on solid or electrolyte concentrations (i.e., 𝑈𝑒𝑞 and 𝜂𝑐) are those that will exhibit a time-domain transient or a frequency response, therefore they will be analyzed in Sections 2.2 and 2.3. The term 𝜂𝑟depends on said concentrations indirectly, so it will be studied subsequently in Section 2.4. The equivalent circuit model along with its discrete state-space representation is presented in Section 2.5 and its correspondence with frequency-domain data is analyzed in Section 2.6. 2.2. Solid diffusion transfer function To obtain the frequency response of the 𝑈𝑒𝑞 term in Eq. (2), the evaluation of the lithium concentration at the surface of the solid particles of the electrodes is necessary. If the particles are assumed to be spherical, the diffusion process taking place within them is described by Eq. (3) [14]: 𝜕𝑐𝑠(𝑟, 𝑡) 𝜕𝑡 =𝐷𝑠 𝑟2 𝜕 𝜕𝑟 (𝑟2𝜕𝑐𝑠(𝑟, 𝑡) 𝜕𝑟 )(3) with the following boundary conditions [14]: 𝜕𝑐𝑠(𝑟, 𝑡) 𝜕𝑟 ||||𝑟=0 = 0,𝜕𝑐𝑠(𝑟, 𝑡) 𝜕𝑟 ||||𝑟=𝑅𝑠 =𝐼𝑅𝑠 3𝐴𝐹 𝐿𝑒𝐷𝑠𝜖𝑠 (4) where 𝑟is the distance along the particle radius, 𝑐𝑠is the lithium concentration in the solid, 𝐷𝑠is the lithium diffusion coefficient in the solid, 𝑅𝑠is the particle radius, 𝐴is the electrode area, 𝐹is Faraday’s constant and 𝐿𝑒is the electrode length. The variation in concentration with respect to the initial value is taken as the equation variable, so as to have zero initial time conditions: 𝑐𝑠=𝑐𝑠−𝑐𝑠,0(5) Upon solving, applying the boundary conditions in Eq. (4) and evaluating at 𝑟=𝑅𝑠, the transfer function from the applied current to the variation in surface concentration is obtained in Eq. (6): 𝑐𝑠,𝑠(𝑠) 𝐼(𝑠)=𝜏𝑠 3𝜖𝑠𝐴𝐹 𝐿𝑒 1 √𝜏𝑠𝑠coth (√𝜏𝑠𝑠)− 1 (6) where 𝜏𝑠is defined as 𝑅2 𝑠 𝐷𝑠. If normalized concentration is considered instead, the transfer function may be rewritten as in Eq. (7): 𝐺𝑠(𝑠) = 𝜒𝑠,𝑠(𝑠) 𝐼(𝑠)=𝐾𝑠 √𝜏𝑠𝑠coth (√𝜏𝑠𝑠)− 1 (7) where 𝜒𝑠,𝑠 =𝑐𝑠,𝑠(𝑠)∕𝑐𝑠,𝑚𝑎𝑥 and 𝐾𝑠=𝜏𝑠 3𝜖𝑠𝐴𝐹 𝐿𝑒𝑐𝑠,𝑚𝑎𝑥 . Taking into account that for small values of 𝑥,√𝑥coth (√𝑥)≈ 1 + 𝑥 3−𝑥2 45 , one can study the behavior of 𝐺𝑠(𝑠)in the limits 𝑠→0 and 𝑠→∞as in Eq. (8): 𝐺𝑠(𝑠) ≈ ⎧ ⎪ ⎨ ⎪ ⎩ 3𝐾𝑠 𝜏𝑠𝑠+𝐾𝑠 5, 𝑠 →0 0, 𝑠 →∞ (8) In order to approximate the frequency response of 𝐺𝑠(𝑠)between these two limits, in this article we propose a transfer function composed of the addition of an integrator and a ZARC element as expressed in Eq. (9): 𝐾𝑠 √𝜏𝑠𝑠coth (√𝜏𝑠𝑠)− 1 ≈3𝐾𝑠 𝜏𝑠𝑠+𝐾𝑠∕5 1 + (𝛽𝑠𝜏𝑠𝑠)𝛼𝑠(9) Therefore, 𝛼𝑠and 𝛽𝑆need to be determined so as to obtain the closest approximation to 𝐺𝑠(𝑠). For this purpose, the frequency response of both transfer functions is plotted in a Nyquist diagram and the weighted impedance error between the exact and the approximate transfer function is minimized within a range of frequencies. In this case, we have considered the interval 𝜔=[1 𝜏𝑠 ,103 𝜏𝑠]due to the fact that at frequencies lower than 2𝜋 𝜏𝑠the integrator behavior is dominant, whereas at frequencies higher than 2𝜋⋅103 𝜏𝑠the amplitude of the frequency response approaches 0. The resulting values are 𝛼𝑠= 0.82 and 𝛽𝑠= 0.0207, and the Nyquist plot of the frequency response of both transfer functions is shown in Fig. 1. Note that Eq. (9) is expressed in terms of the normalized frequency 𝜏𝑠𝑠, so this approximation remains valid regardless of the specific physical parameter values. Furthermore, the fact that the identified order exponent is not equal to 1 proves that a parallel-RC element would not be the most appropriate alternative for approximating the original transfer function. Particularizing for the negative and positive particle, the term 𝑈𝑒𝑞 from Eq. (2a) may be expressed as in Eq. (10), taking into account that the boundary condition for the positive particle is negative if the current applied to the cell 𝐼(𝑠)is considered positive while charging: 𝑈𝑒𝑞(𝑠) = 𝑂𝐶𝑃𝑝⎛⎜⎜⎝⎡⎢⎢⎣ −3𝐾𝑝 𝜏𝑝𝑠−𝐾𝑝∕5 1 + (0.0207𝜏𝑝𝑠)0.82 ⎤⎥⎥⎦ 𝐼(𝑠)⎞⎟⎟⎠ − 𝑂𝐶𝑃𝑛([3𝐾𝑛 𝜏𝑛𝑠+𝐾𝑛∕5 1 + (0.0207𝜏𝑛𝑠)0.82 ]𝐼(𝑠))(10) where 𝜏𝑝=𝑅2 𝑝 𝐷𝑝,𝜏𝑛=𝑅2 𝑛 𝐷𝑛,𝐾𝑝=𝜏𝑝 3𝜖𝑝𝐴𝐹 𝐿𝑝𝑐𝑝,𝑚𝑎𝑥 and 𝐾𝑛=𝜏𝑛 3𝜖𝑛𝐴𝐹 𝐿𝑛𝑐𝑛,𝑚𝑎𝑥 . Journal of Energy Storage 64 (2023) 107150 4 P. Rodríguez-Iturriaga et al. Fig. 1. Nyquist plot of the frequency response of the exact and approximate normalized transfer functions in Eq. (9) for solid diffusion. Note that the real and imaginary axes have different scales. 2.3. Electrolyte diffusion transfer function Next, the frequency response of the 𝜂𝑐term in Eq. (2) must be determined in order to calculate the overpotential due to concentration gradients in the electrolyte. Consequently, the electrode-averaged lithium concentration in the electrolyte has to be evaluated. For this purpose, the diffusion process described in Eq. (11) is analyzed [14]: 𝜖𝑒,𝑛 𝜕𝑐𝑒,𝑛(𝑥, 𝑡) 𝜕𝑡 =𝜖𝑏 𝑒,𝑛𝐷𝑒 𝜕2𝑐𝑒,𝑛(𝑥, 𝑡) 𝜕𝑥2− (1 − 𝑡+)𝐼 𝐴𝐹 𝐿𝑛 ,0< 𝑥 < 𝐿𝑛(11a) 𝜖𝑒,𝑠 𝜕𝑐𝑒,𝑠(𝑥, 𝑡) 𝜕𝑡 =𝜖𝑏 𝑒,𝑠𝐷𝑒 𝜕2𝑐𝑒,𝑠(𝑥, 𝑡) 𝜕𝑥2, 𝐿𝑛< 𝑥 < 𝐿𝑛+𝐿𝑠(11b) 𝜖𝑒,𝑝 𝜕𝑐𝑒,𝑝(𝑥, 𝑡) 𝜕𝑡 =𝜖𝑏 𝑒,𝑝𝐷𝑒 𝜕2𝑐𝑒,𝑝(𝑥, 𝑡) 𝜕𝑥2− (1 − 𝑡+)𝐼 𝐴𝐹 𝐿𝑝 , 𝐿𝑛+𝐿𝑠< 𝑥 < 𝐿 (11c) where 𝐷𝑒is the diffusion coefficient for lithium in the electrolyte, x is the linear distance along the cell length 𝐿=𝐿𝑛+𝐿𝑠+𝐿𝑝, and the reference 𝑥= 0 is placed at the interface between the anode and the current collector. The boundary conditions at the endpoints of the cell are the following [14]: 𝜕𝑐𝑒,𝑛(𝑥, 𝑡) 𝜕𝑥 ||||𝑥=0 = 0,𝜕𝑐𝑒,𝑝(𝑥, 𝑡) 𝜕𝑥 ||||𝑥=𝐿 = 0,(12) as well as continuity in concentration and flux between the different domains [14]: 𝑐𝑒,𝑛(𝑥, 𝑡)|||𝑥=𝐿𝑛 =𝑐𝑒,𝑠(𝑥, 𝑡)|||𝑥=𝐿𝑛 , 𝑐𝑒,𝑠(𝑥, 𝑡)|||𝑥=𝐿𝑛+𝐿𝑠 =𝑐𝑒,𝑝(𝑥, 𝑡)|||𝑥=𝐿𝑛+𝐿𝑠 𝜖𝑏 𝑛 𝜕𝑐𝑒,𝑛(𝑥,𝑡) 𝜕𝑥 ||||𝑥=𝐿𝑛 =𝜖𝑏 𝑠 𝜕𝑐𝑒,𝑠(𝑥,𝑡) 𝜕𝑥 ||||𝑥=𝐿𝑛 , 𝜖𝑏 𝑛 𝜕𝑐𝑒,𝑠(𝑥,𝑡) 𝜕𝑥 ||||𝑥=𝐿𝑛+𝐿𝑠 =𝜖𝑏 𝑝 𝜕𝑐𝑒,𝑝(𝑥,𝑡) 𝜕𝑥 ||||𝑥=𝐿𝑛+𝐿𝑠 (13) Proceeding as in Section 2.2, the variation in concentration with respect to the typical lithium concentration in the electrolyte 𝑐𝑒,𝑡𝑦𝑝 is taken as the equation variable: 𝑐𝑒,𝑘(𝑥, 𝑡) = 𝑐𝑒,𝑘(𝑥, 𝑡) − 𝑐𝑒,𝑡𝑦𝑝 (14) Additionally, the following spatial variables are considered for convenience from now on: 𝑥𝑛=𝑥, 𝑥𝑛∈ [0, 𝐿𝑛], 𝑥𝑠=𝑥−𝐿𝑛, 𝑥𝑠∈ [0, 𝐿𝑠], 𝑥𝑝=𝐿−𝑥, 𝑥𝑝∈ [0, 𝐿𝑝] (15) Taking the Laplace transform of Eq. (11) yields the following set of ordinary differential equations: 𝜕2𝑐𝑒,𝑛(𝑥𝑛, 𝑠) 𝜕𝑥2 𝑛 −𝑠 𝜖𝑏−1 𝑒,𝑛 𝐷𝑒 𝑐𝑒,𝑛(𝑥𝑛, 𝑠) = (1 − 𝑡+) 𝐴𝐹 𝐿𝑛𝐷𝑒𝜖𝑏 𝑒,𝑛 𝐼(𝑠)(16a) 𝜕2𝑐𝑒,𝑠(𝑥𝑠, 𝑠) 𝜕𝑥2 𝑠 −𝑠 𝜖𝑏−1 𝑒,𝑠 𝐷𝑒 𝑐𝑒,𝑠(𝑥𝑠, 𝑠)=0 (16b) 𝜕2𝑐𝑒,𝑝(𝑥𝑝, 𝑠) 𝜕𝑥2 𝑝 −𝑠 𝜖𝑏−1 𝑒,𝑝 𝐷𝑒 𝑐𝑒,𝑝(𝑥𝑝, 𝑠) = (1 − 𝑡+) 𝐴𝐹 𝐿𝑝𝐷𝑒𝜖𝑏 𝑒,𝑝 𝐼(𝑠)(16c) Solving Eq. (16) and applying the boundary conditions at the endpoints of the cell in Eq. (12) yields Eq. (17): 𝑐𝑒,𝑛(𝑥𝑛, 𝑠)=2𝐶𝑛cosh ⎛⎜⎜⎝ 𝑥𝑛√𝑠 𝜖𝑏−1 𝑒,𝑛 𝐷𝑒⎞⎟⎟⎠ −(1 − 𝑡+) 𝐴𝐹 𝐿𝑛𝜖𝑒,𝑛 𝐼(𝑠) 𝑠(17a) 𝑐𝑒,𝑠(𝑥𝑠, 𝑠) = 𝐶𝑠,1exp ⎛⎜⎜⎝ 𝑥𝑠√𝑠 𝜖𝑏−1 𝑒,𝑠 𝐷𝑒⎞⎟⎟⎠ +𝐶𝑠,2exp ⎛⎜⎜⎝ −𝑥𝑠√𝑠 𝜖𝑏−1 𝑒,𝑠 𝐷𝑒⎞⎟⎟⎠ (17b) 𝑐𝑒,𝑝(𝑥𝑝, 𝑠)=2𝐶𝑝cosh ⎛⎜⎜⎝ 𝑥𝑝√𝑠 𝜖𝑏−1 𝑒,𝑝 𝐷𝑒⎞⎟⎟⎠ −(1 − 𝑡+) 𝐴𝐹 𝐿𝑝𝜖𝑒,𝑝 𝐼(𝑠) 𝑠(17c) where 𝐶𝑛,𝐶𝑠,1,𝐶𝑠,2and 𝐶𝑝have to be determined from the boundary conditions in Eq. (13). In this article, instead of solving for the coefficients directly, we propose a method to simplify the system of equations in Eq. (17) by deriving the required conditions to reduce the model. It has to be noted that the consideration of spatial variables as in Eq. (15) allows for the definition of the following timescales: 𝜏𝑒,𝑛 =𝐿2 𝑛 𝜖𝑏−1 𝑒,𝑛 𝐷𝑒 , 𝜏𝑒,𝑠 =𝐿2 𝑠 𝜖𝑏−1 𝑒,𝑠 𝐷𝑒 , 𝜏𝑒,𝑝 =𝐿2 𝑝 𝜖𝑏−1 𝑒,𝑝 𝐷𝑒 (18) The ratio between the lithium migration timescales in the separator and the electrodes is determined by Eq. (19): 𝜏𝑒,𝑠 𝜏𝑒,𝑛,𝑝 = 𝐿2 𝑠𝜖𝑏−1 𝑒,𝑛,𝑝 𝐿2 𝑛,𝑝𝜖𝑏−1 𝑒,𝑠 (19) If this ratio is sufficiently close to 0, the electrolyte concentration in the separator may be considered to be in steady state with respect to that in the electrodes. This results in a linear concentration profile along the separator length, by taking the limit 𝑠→0in Eq. (16b). The slope 𝑚of the concentration in the separator is determined by the boundary condition regarding flux continuity in Eq. (13): 𝑚= 2𝐶𝑛 𝜖𝑏 𝑒,𝑛 𝜖𝑏 𝑒,𝑠 √𝑠 𝜖𝑏−1 𝑒,𝑛 𝐷𝑒 sinh (√𝜏𝑒,𝑛𝑠)(20) A further simplification is carried out by calculating the ratio between the concentration at the electrode-separator interface and the total concentration increment in the negative electrode. The concentration at the electrode-separator interface is calculated assuming that it is equal to half of the total variation in the separator: |||𝑐𝑒,𝑠|𝑥𝑠=0|||=𝑚𝐿𝑠 2=𝐶𝑛 𝜖𝑏 𝑒,𝑛 𝜖𝑏 𝑒,𝑠 √ √ √ √𝑠𝐿2 𝑠 𝜖𝑏−1 𝑒,𝑛 𝐷𝑒 sinh (√𝜏𝑒,𝑛𝑠)(21) The concentration increment in the negative electrode is calculated directly by substituting in Eq. (17a): |||𝛥 𝑐𝑒,𝑛|||= 2𝐶𝑛[cosh (√𝜏𝑒,𝑛𝑠)− 1](22) Therefore, the corresponding ratio is calculated as follows: |||| 𝑐𝑒,𝑠|𝑥𝑠=0 𝛥 𝑐𝑒,𝑠 ||||=1 2 𝜖𝑏 𝑒,𝑛 𝜖𝑏 𝑒,𝑠 sinh (√𝜏𝑒,𝑛𝑠) cosh (√𝜏𝑒,𝑛𝑠)− 1 √ √ √ √𝑠𝐿2 𝑠 𝐷𝑒𝜖𝑏−1 𝑒,𝑛 =1 2 𝜖𝑏 𝑒,𝑛 𝜖𝑏 𝑒,𝑠 coth (√𝜏𝑒,𝑛𝑠 2)√ √ √ √𝑠𝐿2 𝑠 𝐷𝑒𝜖𝑏−1 𝑒,𝑛 (23) Journal of Energy Storage 64 (2023) 107150 5 P. Rodríguez-Iturriaga et al. One can verify that this ratio reduces to 𝐿𝑠𝜖𝑏 𝑒,𝑛 𝐿𝑛𝜖𝑏 𝑒,𝑠 for the frequency range of interest around 𝜔=1 𝜏𝑒,𝑛 . Taking this into account, if 𝐿𝑠𝜖𝑏 𝑒,𝑛,𝑝 𝐿𝑛,𝑝𝜖𝑏 𝑒,𝑠 ≪ 1, the variations in concentration in the separator may be neglected in comparison to those in the electrodes and, as a result, the boundary condition for the concentration in the negative electrode may be approximated by 𝑐𝑒,𝑛(𝑥𝑛, 𝑠)|||𝑥𝑛=𝐿𝑛 = 0. Upon applying said boundary condition, the variation in the electrolyte concentration is obtained as a function of 𝑥𝑛and 𝑠: 𝑐𝑒,𝑛(𝑥𝑛, 𝑠) = (1 − 𝑡+) 𝐴𝐹 𝐿𝑛𝜖𝑒,𝑛 𝐼(𝑠) 𝑠⎡⎢⎢⎢⎢⎣ cosh (𝑥𝑛√𝑠 𝜖𝑏−1 𝑒,𝑛 𝐷𝑒) cosh (√𝜏𝑒,𝑛𝑠)− 1⎤⎥⎥⎥⎥⎦ (24) Next, the electrode-averaged value of the variation in electrolyte concentration is calculated: 𝑐𝑒,𝑛(𝑠) = 1 𝐿𝑛∫𝐿𝑛 0 𝑐𝑒,𝑛(𝑥𝑛, 𝑠)𝑑𝑥𝑛 =(1 − 𝑡+) 𝐴𝐹 𝐿𝑛𝜖𝑒,𝑛 𝐼(𝑠) 𝑠⎡⎢⎢⎢⎣ 1 √𝜏𝑒,𝑛𝑠coth (√𝜏𝑒,𝑛𝑠)− 1⎤⎥⎥⎥⎦ (25) Therefore, the transfer function from the applied current to the electrode-averaged variation in electrolyte concentration may be expressed as: 𝐺𝑒(𝑠) = 𝑐𝑒,𝑛(𝑠) 𝐼(𝑠)=𝐾𝑒,𝑛 ⎡⎢⎢⎢⎣ 1 − √𝜏𝑒,𝑛𝑠coth (√𝜏𝑒,𝑛𝑠) 𝜏𝑒,𝑛𝑠√𝜏𝑒,𝑛𝑠coth (√𝜏𝑒,𝑛𝑠)⎤⎥⎥⎥⎦ (26) where 𝐾𝑒,𝑛 =(1−𝑡+)𝜏𝑒,𝑛 𝐴𝐹 𝐿𝑛𝜖𝑒,𝑛 . Proceeding as in Section 2.2, one can study the behavior of 𝐺𝑒(𝑠)in the limits 𝑠→0and 𝑠→∞as in Eq. (27): 𝐺𝑒(𝑠) ≈ ⎧ ⎪ ⎨ ⎪ ⎩ −𝐾𝑒,𝑛 3, 𝑠 →0 0, 𝑠 →∞ (27) In order to approximate the frequency response of 𝐺𝑒(𝑠)between these two limits, we propose a ZARC element as expressed in Eq. (28): 𝐾𝑒,𝑛 ⎡⎢⎢⎢⎣ 1 − √𝜏𝑒,𝑛𝑠coth (√𝜏𝑒,𝑛𝑠) 𝜏𝑒,𝑛𝑠√𝜏𝑒,𝑛𝑠coth (√𝜏𝑒,𝑛𝑠)⎤⎥⎥⎥⎦ ≈ − 𝐾𝑒,𝑛∕3 1 + (𝛽𝑒𝜏𝑒,𝑛𝑠)𝛼𝑒(28) As in Section 2.2,𝛼𝑒and 𝛽𝑒must be determined so as to obtain the closest approximation to 𝐺𝑒(𝑠). For this purpose, the frequency response of both transfer functions is plotted in a Nyquist diagram and the impedance error between the exact and the approximate transfer function is minimized within a range of frequencies. In this case, we have considered the interval 𝜔=[1 103𝜏𝑒,𝑛 ,103 𝜏𝑒,𝑛 ]and the resulting values are 𝛼𝑒= 0.9936 and 𝛽𝑒= 0.3983. The Nyquist plot of the frequency response of both transfer functions with a positive real part is shown in Fig. 2. A similar result is obtained for the positive electrode with a plus sign. As in Section 2.2, Eq. (28) is expressed in terms of the normalized frequency 𝜏𝑒,𝑛𝑠, so this approximation remains valid regardless of the specific physical parameter values. In the interest of simplicity, a value of 𝛼𝑒= 1 will be used further in this article. Lastly, the overpotential due to concentration gradients in the electrolyte 𝜂𝑐is calculated according to Eq. (2b): 𝜂𝑐(𝑠) = 2𝑅𝑇 𝐹 (1 − 𝑡+) 𝑐𝑒,𝑡𝑦𝑝 (𝑐𝑒,𝑝 −𝑐𝑒,𝑛) =2𝑅𝑇 𝐹 (1 − 𝑡+) 𝑐𝑒,𝑡𝑦𝑝 (𝐾𝑒,𝑝∕3 1 + 𝛽𝜏𝑒,𝑝𝑠+𝐾𝑒,𝑛∕3 1 + 𝛽𝜏𝑒,𝑛𝑠)𝐼(𝑠) (29) Fig. 2. Nyquist plot of the frequency response of the exact and approximate normalized transfer functions in Eq. (28) for diffusion in the electrolyte. If 𝐿𝑛≈𝐿𝑝and 𝜖𝑒,𝑛 ≈𝜖𝑒,𝑝, the terms corresponding to both electrodes in Eq. (29) may be combined into one single RC network. 2.4. Charge transfer overpotential From the expressions for the solid and electrolyte concentrations determined in Sections 2.2 and 2.3, an accurate approximation of the reaction overpotentials may be obtained according to Eq. (2c). For this purpose, the exchange current densities 𝑗0,𝑛 and 𝑗0,𝑝 are calculated as follows [14]: 𝑗0,𝑛 =1 𝐿𝑛∫𝐿𝑛 0𝑚𝑛𝑐𝑛,𝑚𝑎𝑥√𝜒𝑠,𝑛(1 − 𝜒𝑠,𝑛)√𝑐𝑒,𝑛𝑑𝑥𝑛 𝑗0,𝑝 =1 𝐿𝑝∫𝐿𝑝 0𝑚𝑝𝑐𝑝,𝑚𝑎𝑥√𝜒𝑠,𝑝(1 − 𝜒𝑠,𝑝)√𝑐𝑒,𝑝𝑑𝑥𝑝 (30) where 𝑚𝑛and 𝑚𝑝are the reaction rates of the negative and positive electrode, respectively. Note that 𝜒𝑠,𝑛 and 𝜒𝑠,𝑝 have been determined already in Section 2.2 and do not depend on the 𝑥dimension in the SPMe. Therefore, the negative electrode exchange current density may be rewritten as follows: 𝑗0,𝑛 =𝑚𝑛𝑐𝑛,𝑚𝑎𝑥√𝜒𝑠,𝑛(1 − 𝜒𝑠,𝑛)√𝑐𝑒,𝑡𝑦𝑝 1 𝐿𝑛∫𝐿𝑛 0√1 + 𝑐𝑒,𝑛 𝑐𝑒,𝑡𝑦𝑝 𝑑𝑥𝑛(31) Assuming that the variation in the electrolyte concentration is sufficiently smaller than the typical concentration 𝑐𝑒,𝑡𝑦𝑝, the following approximation can be made: 𝑗0,𝑛 ≈𝑚𝑛𝑐𝑛,𝑚𝑎𝑥√𝜒𝑠,𝑛(1 − 𝜒𝑠,𝑛)√𝑐𝑒,𝑡𝑦𝑝 1 𝐿𝑛∫𝐿𝑛 0(1 + 𝑐𝑒,𝑛 2𝑐𝑒,𝑡𝑦𝑝 )𝑑𝑥𝑛(32) Consequently, the negative electrode exchange current density may be rewritten as a function of the electrode-averaged variation in the electrolyte concentration 𝑐𝑒,𝑛, which was determined in the previous Section: 𝑗0,𝑛 ≈𝑚𝑛𝑐𝑛,𝑚𝑎𝑥√𝜒𝑠,𝑛(1 − 𝜒𝑠,𝑛)√𝑐𝑒,𝑡𝑦𝑝 (1 + 𝑐𝑒,𝑛 2𝑐𝑒,𝑡𝑦𝑝 )(33) From this expression of the exchange current density, the chargetransfer overpotential may be calculated according to Eq. (2c): 𝜂𝑟,𝑛 =2𝑅𝑇 𝐹sinh−1 (𝐼𝑅𝑛 6𝑗0,𝑛𝜖𝑛𝐿𝑛𝐴)(34) This may also be expressed via the charge-transfer resistance, including its explicit dependence on the applied current as in Eq. (35) Journal of Energy Storage 64 (2023) 107150 6 P. Rodríguez-Iturriaga et al. Fig. 3. (a) Full equivalent circuit model. (b) Simplified equivalent circuit model. 𝑅𝑐𝑡,𝑛(𝐼) = 2𝑅𝑇 𝐹 𝐼0,𝑛 [𝐼0,𝑛 𝐼sinh−1 (𝐼 𝐼0,𝑛 )], 𝐼0,𝑛 =6𝑗0,𝑛𝜖𝑛𝐿𝑛𝐴 𝑅𝑛 (35) A similar result is obtained for the positive electrode. We consider this to be a major advantage of employing electrode-averaged quantities as in [14], which allows for an accurate estimation of the reaction overpotential from only the electrode-averaged electrolyte concentration instead of having to solve for the concentration profile along the 𝑥dimension. Taking into account the current dependence of the reaction overpotentials is especially critical in applications where a wide range of current rates is expected. 2.5. Equivalent circuit model According to the output voltage expression in Eq. (1) and the resulting equations for the individual terms in Eq. (2), the equivalent circuit model in Fig. 3-(a) may be constructed, where the resistors 𝑅𝑜ℎ𝑚,𝑒 = 1 𝜅𝐴 (𝐿𝑛 3𝜖𝑏 𝑒,𝑛 +𝐿𝑠 𝜖𝑏 𝑒,𝑠 +𝐿𝑝 3𝜖𝑏 𝑒,𝑝 )and 𝑅𝑜ℎ𝑚,𝑠 =1 3𝐴(𝐿𝑛 𝜎𝑛 +𝐿𝑝 𝜎𝑝)correspond to the terms in Eqs. (2d) and (2e). For the time-domain implementation of the circuit, a set of serially connected parallel RC branches will be employed as an approximation Table 1 Parameter equations for the 7-RC approximation as a function of 𝛼[53]. Parameter Expression 𝑟1=𝑟70.14(1 − 𝛼)2 𝑟2=𝑟60.22(1 − 𝛼)−0.08(1 − 𝛼)3 𝑟3=𝑟5(0.12 + 0.057𝑒3.4𝛼)(1 − 𝛼) 𝑟41−2(𝑟1+𝑟2+𝑟3) 𝑡1=1/𝑡71.4⋅10−8𝑒19𝛼(1.6−𝛼) 𝑡2=1/𝑡6 0.078𝛼5.63 0.026+𝛼3.67 𝑡3=1/𝑡5 0.56𝛼2.27 0.4+𝛼1.3 𝑡41 of the ZARC element. In particular, we will use the continuous approximation by 7 RC elements that we presented in [53]. This approach is the most appropriate for the time-domain simulation of the ZARC element as we showed in [56], due to the fact that it avoids the issues caused by the specification of a memory length in the Grünwald– Letnikov approach. Taking the fractional order 𝛼, the resistor 𝑅𝑍𝐴𝑅𝐶 and the time constant 𝜏𝑍𝐴𝑅𝐶 as inputs, the parameter values are directly calculated as 𝑅𝑖=𝑅𝑍𝐴𝑅𝐶 ⋅𝑟𝑖(𝛼)and 𝜏𝑖=𝜏𝑍𝐴𝑅𝐶 ⋅𝑡𝑖(𝛼)for seven-element network according to Table 1 [53]: Journal of Energy Storage 64 (2023) 107150 7 P. Rodríguez-Iturriaga et al. As a result, the approximate discrete state-space representation of a ZARC element is shown in Eq. (36): 𝑥𝑍𝐴𝑅𝐶 (𝑘) = [𝑖1(𝑘)𝑖2(𝑘) … 𝑖7(𝑘)]𝑇 𝑥𝑍𝐴𝑅𝐶 (𝑘) = 𝐴𝑍𝐴𝑅𝐶 𝑥𝑍𝐴𝑅𝐶 (𝑘− 1) + 𝐵𝑍𝐴𝑅𝐶 𝑢(𝑘− 1) 𝑦𝑍𝐴𝑅𝐶 (𝑘) = 𝐶𝑍𝐴𝑅𝐶 𝑥𝑍𝐴𝑅𝐶 (𝑘) (36) where 𝑥𝑍𝐴𝑅𝐶 is the state vector, 𝑢is the input current and 𝑦𝑍𝐴𝑅𝐶 is the total voltage difference in the subcircuit, with the following state-space representation matrices: 𝐴𝑍𝐴𝑅𝐶 =𝑑𝑖𝑎𝑔 [exp (−𝛥𝑡 𝜏1)exp (−𝛥𝑡 𝜏2)… exp (−𝛥𝑡 𝜏7)] 𝐵𝑍𝐴𝑅𝐶 =[1 − exp (−𝛥𝑡 𝜏1)1 − exp (−𝛥𝑡 𝜏2)… 1 − exp (−𝛥𝑡 𝜏7)]𝑇 𝐶𝑍𝐴𝑅𝐶 =[𝑅1𝑅2…𝑅7](37) where 𝛥𝑡 is the sampling time. Consequently, the state vector and the state equation of the ECM are as follows: 𝑥(𝑘) = [𝜒𝑝(𝑘)𝑥𝑍𝐴𝑅𝐶,𝑝(𝑘)𝑐𝑒,𝑝(𝑘)𝜒𝑛(𝑘)𝑥𝑍𝐴𝑅𝐶,𝑛(𝑘)𝑐𝑒,𝑛(𝑘)]𝑇 𝑥(𝑘) = 𝐴𝑥(𝑘− 1) + 𝐵𝑢(𝑘− 1), 𝑢(𝑘) = 𝐼(𝑘) (38) where 𝐴=𝑑𝑖𝑎𝑔 [1𝑑𝑖𝑎𝑔(𝐴𝑍𝐴𝑅𝐶,𝑝) exp (−𝛥𝑡 𝛽𝑒𝜏𝑒,𝑝 )1𝑑𝑖𝑎𝑔(𝐴𝑍𝐴𝑅𝐶,𝑛) exp (−𝛥𝑡 𝛽𝑒𝜏𝑒,𝑛 )] 𝐵= ⎡⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣ −𝐾𝑝𝛥𝑡 3𝜏𝑝 −𝐵𝑍𝐴𝑅𝐶,𝑝 𝐾𝑒,𝑝 3[1 − exp (−𝛥𝑡 𝛽𝑒𝜏𝑒,𝑝 )] 𝐾𝑛𝛥𝑡 3𝜏𝑛 𝐵𝑍𝐴𝑅𝐶,𝑛 −𝐾𝑒,𝑛 3[1 − exp (−𝛥𝑡 𝛽𝑒𝜏𝑒,𝑛 )] ⎤⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦(39) Note the minus sign in the terms corresponding to the positive particle and the electrolyte concentration in the negative electrode in matrix 𝐵. Before evaluating the output voltage expression, it is convenient to calculate explicitly the surface concentration of both particles as well as the exchange current density in both electrodes: 𝜒𝑠,𝑝(𝑘) = 𝜒𝑝(𝑘) + 𝐶𝑍𝐴𝑅𝐶,𝑝𝑥𝑍𝐴𝑅𝐶,𝑝(𝑘) 𝜒𝑠,𝑛(𝑘) = 𝜒𝑛(𝑘) + 𝐶𝑍𝐴𝑅𝐶,𝑛𝑥𝑍𝐴𝑅𝐶,𝑛(𝑘) 𝑗0,𝑝(𝑘) = 𝑚𝑝𝑐𝑝,𝑚𝑎𝑥√𝜒𝑠,𝑝(𝑘)(1 − 𝜒𝑠,𝑝(𝑘))√𝑐𝑒,𝑡𝑦𝑝 (1 + 𝑐𝑒,𝑝(𝑘) 2𝑐𝑒,𝑡𝑦𝑝 ) 𝑗0,𝑛(𝑘) = 𝑚𝑛𝑐𝑛,𝑚𝑎𝑥√𝜒𝑠,𝑛(𝑘)(1 − 𝜒𝑠,𝑛(𝑘))√𝑐𝑒,𝑡𝑦𝑝 (1 + 𝑐𝑒,𝑛(𝑘) 2𝑐𝑒,𝑡𝑦𝑝 ) (40) Taking all the previous equations into account, the output voltage expression is as follows: 𝑣(𝑘) = 𝑂𝐶𝑃𝑝(𝜒𝑠,𝑝(𝑘)) − 𝑂𝐶𝑃𝑛(𝜒𝑠,𝑛(𝑘)) + 2𝑅𝑇 𝐹 (1 − 𝑡+) 𝑐𝑒,𝑡𝑦𝑝 (𝑐𝑒,𝑝(𝑘) − 𝑐𝑒,𝑛(𝑘)) +2𝑅𝑇 𝐹[sinh−1 (𝐼𝑅𝑝 6𝑗0,𝑝(𝑘)𝜖𝑝𝐿𝑝𝐴)+ sinh−1 (𝐼𝑅𝑛 6𝑗0,𝑛(𝑘)𝜖𝑛𝐿𝑛𝐴)] +𝐼 𝜅𝐴 (𝐿𝑛 3𝜖𝑏 𝑒,𝑛 +𝐿𝑠 𝜖𝑏 𝑒,𝑠 +𝐿𝑝 3𝜖𝑏 𝑒,𝑝 ) +𝐼 3𝐴(𝐿𝑛 𝜎𝑛 +𝐿𝑝 𝜎𝑝)(41) If both electrodes are assumed to have similar physical and spatial properties, the ECM can be further simplified to that in Fig. 3-(b), thus reducing in half the number of required states. The developed ECM may be interpreted as the combination of the two discussed approaches in Section 1: on the one hand, the transcendental transfer functions are obtained directly from the electrochemical model and then approximated with fractional-order transfer functions instead of Padé’s approximant [40–43] in our case. On the other hand, the standard RC equivalent circuit model from [36] is split into two sections as in [39], so that the leftmost one accounts for the lithium diffusion within the active material particles. Consequently, the lithium concentration is evaluated at their surface, thus establishing a more accurate correspondence with the manner in which the electrode potential is calculated in electrochemical models. However, a better description of the diffusion process is achieved with a ZARC element instead of a RC network, as indicated by the optimal order exponent 𝛼= 0.82. Therefore, this model, stated in this way, may be useful for most of the applications where ECMs are employed providing additional physical insight, and may also be implemented in standard circuit simulators by substituting the ZARC element by its multiple-RC approximation [53], thus allowing for a simple simulation of constantvoltage phases by substituting the current source in Fig. 3 by a voltage source. 2.6. Frequency domain application From the obtained transfer functions, it is possible to determine the medium and low frequency response of the cell directly given its physical parameters. However, the high-frequency transient effects of the charge transfer processes are usually not considered in the DFN or SPMe due to them being orders of magnitude faster than diffusion dynamics. Consequently, in order to accurately reproduce EIS data, a CPE is usually connected in parallel with the charge-transfer resistance 𝑅𝑐𝑡 in the circuit shown in Fig. 3-(b). The small-signal value of 𝑅𝑐𝑡 is calculated by taking into account that for small values of 𝑥, 𝑥−1 sinh−1 (𝑥)≈ 1: 𝑅𝑐𝑡(𝐼≈ 0) = 2𝑅𝑇 𝐹 𝐼0 , 𝑍𝑐𝑡(𝑠) = 𝑅𝑐𝑡 1 + (𝜏𝑐𝑡𝑠)𝛼𝑐𝑡 (42) Conversely, if the frequency-domain behavior of the battery cell is characterized experimentally in an EIS test, there will be a substantial overlap between the effects of both electrodes and only the effective parameters for the simplified ECM in Fig. 3-(b) will be identifiable. Given that an EIS test consists of small-signal variations around a certain operating point, the ECM is linearized around 𝑆𝑂𝐶 =𝑆𝑂𝐶𝐸𝐼𝑆 . Using the fact that cell SOC remains unchanged, one can write: 𝑂𝐶𝑉 (𝜒𝑠) = 𝑂𝐶𝑉 (𝑆𝑂𝐶𝐸𝐼𝑆 ) + 𝜕𝑂𝐶𝑉 𝜕𝑆𝑂𝐶 ||||𝑆𝑂𝐶𝐸𝐼𝑆 𝜒𝑠(43) where 𝜕𝑂𝐶𝑉 𝜕𝑆𝑂𝐶 |||𝑆𝑂𝐶𝐸𝐼𝑆 is determined empirically from the OCV-SOC relationship and the frequency response of 𝜒𝑠may be approximated by that of a ZARC element. Therefore, a possible transfer function for fitting EIS data is composed of the addition of an ohmic resistor, a highfrequency ZARC element for charge transfer processes, a mid-frequency RC element for the diffusion in the electrolyte and a low-frequency ZARC element for the solid diffusion, as shown in Eq. (44): 𝑍(𝑠) = 𝑅𝑜ℎ𝑚 +𝑅𝑐𝑡 1 + (𝜏𝑐𝑡𝑠)𝛼𝑐𝑡 +𝑅𝑒 1 + 𝜏𝑒𝑠+𝜕𝑂𝐶𝑉 𝜕𝑆𝑂𝐶 ||||𝑆𝑂𝐶𝐸𝐼𝑆 𝐾𝑠 1 + (𝜏𝑠𝑠)𝛼𝑠(44) where and the parameters to be identified are 𝑅𝑜ℎ𝑚,𝑅𝑐𝑡,𝜏𝑐𝑡,𝛼𝑐𝑡,𝑅𝑒, 𝜏𝑒,𝐾𝑠and 𝜏𝑠, although 𝜏𝑐𝑡 and 𝛼𝑐𝑡 will not be used in the ECM. Additionally, the value of the exchange current 𝐼0may be calculated directly from 𝑅𝑐𝑡, so as to take into account the current dependence of the charge transfer resistance for large-signal operation. We believe that the correspondence between the time and frequency-domain behaviors of the cell is an advantage over electrochemical models, as it also provides a consistent characterization method since physical parameters are grouped in resistors and time constants. Journal of Energy Storage 64 (2023) 107150 8 P. Rodríguez-Iturriaga et al. Table 2 Parameter set from [14]. Parameter Units Description Value Negative electrode 𝐿𝑛mThickness 1⋅10−4 𝑅𝑛mParticle radius 1⋅10−5 𝐷𝑛m2∕s Solid diffusivity 3.9⋅10−14 𝜖𝑛– Active material volume fraction 0.6 𝑐𝑛,𝑚𝑎𝑥 mol∕m3Maximum lithium concentration 2.4983 ⋅104 𝜎𝑛S∕m Solid conductivity 100 𝜖𝑒,𝑛 – Porosity 0.3 𝑚𝑛(A∕m2)(m3∕mol)1.5Reaction rate 2⋅10−5 Positive electrode 𝐿𝑝mThickness 1⋅10−4 𝑅𝑝mParticle radius 1⋅10−5 𝐷𝑝m2∕s Solid diffusivity 1⋅10−13 𝜖𝑝– Active material volume fraction 0.5 𝑐𝑝,𝑚𝑎𝑥 mol∕m3Maximum lithium concentration 5.1218 ⋅104 𝜎𝑝S∕m Solid conductivity 10 𝜖𝑒,𝑝 – Porosity 0.3 𝑚𝑝(A∕m2)(m3∕mol)1.5Reaction rate 6⋅10−7 Separator 𝐿𝑠mThickness 2.5⋅10−5 𝜖𝑒,𝑠 – Porosity 1 Overall 𝑐𝑒,𝑡𝑦𝑝 mol∕m3Typical electrolyte concentration 1⋅103 𝐷𝑒m2∕s Typical electrolyte diffusivity 5.34 ⋅10−10 𝜅S∕m Typical electrolyte conductivity 1.1 𝑡+– Transference number 0.4 𝑏– Bruggeman coefficient 1.5 𝐴m2Electrode area 2.8359 ⋅10−2 𝑄Ah Cell capacity 0.68 3. Theoretical results and discussion In this section, the performance of the proposed ECM is validated against the SPMe from [14] for different current profiles. For the simulation of the electrochemical model, we have used PyBaMM (Python Battery Mathematical Modeling) [57]. PyBaMM is a battery modeling software implemented in Python designed to simplify the comparison of standard battery models by providing an interface to discretization methods and numerical solvers. In this case, we have employed a typical discretization consisting of 20 points in each domain as well as 20 points for both particles. The set of physical parameters is shown in Table 2. Before carrying out the simulations, the validity of the approximations derived in Section 2.3 is verified: 𝐿2 𝑠𝜖𝑏−1 𝑒,𝑛,𝑝 𝐿2 𝑛,𝑝𝜖𝑏−1 𝑒,𝑠 = 0.034 ≪1, 𝐿𝑠𝜖𝑏 𝑒,𝑛,𝑝 𝐿𝑛,𝑝𝜖𝑏 𝑒,𝑠 = 0.041 ≪1(45) Next, the performance of the proposed ECM has been validated against the cited SPMe model by comparing their simulation results for distinct operation scenarios. In this case, we have employed four constant-current discharges at 2C, 1C, C/2 and C/5 followed by a 30min rest where C-rate is the measurement of the charge and discharge current with respect to its nominal capacity. Additionally, we have also considered a driving cycle (US06) in order to test our model under dynamic operating conditions. The simulation results for the constant-current discharges are shown in Fig. 4, plotted with respect to normalized time as a function of the Crate. It is observed that the proposed model is able to reproduce the cell voltage accurately during the discharge regardless of the current rate. We attribute this mainly to the precise calculation of the charge-transfer reaction overpotentials as detailed in Section 2.4, which present a highly nonlinear dependency on current. Furthermore, the relaxation profile is also accurately modeled due to the introduction of fractionalorder circuit elements to account for the solid diffusion process. The voltage response of the battery and the US06 current profile are shown in Fig. 5. It can also be observed that the proposed ECM is able to yield a greatly accurate terminal voltage in dynamic conditions. Furthermore, the comparison between the battery internal states and their equivalent in the proposed ECM is shown in Fig. 6. Given that the proposed model has been analytically derived from the SPMe, it is able to provide information about the internal states of the battery. We believe this is a qualitative advantage over standard equivalent circuit models [26], which are constructed to reproduce battery voltage only. Nevertheless, a thorough comparison between the accuracy of different ECMs with respect to experimental data constitutes a separate study that is beyond the scope of this article. There could be many factors involved in the discrepancies between model results and experimental data, such as temperature, hysteresis and current-rate effects, as well as the parameterization process and the operating conditions under which the models are tested. The error results between the SPMe and the proposed ECM are summarized in Table 3. It has to be pointed out that not only does the model provide an accurate approximation of the terminal voltage, but also the internal states of the battery, namely the surface concentration of the solid particles and the electrode-averaged electrolyte concentration, which results in a precise calculation of the charge-transfer reaction overpotentials. The maximum voltage error occurs immediately after the end of the discharge, due to the fact that the proposed approximation for the solid diffusion transfer function is less accurate at higher frequencies, as observed in Fig. 1. This also explains the slightly higher error in the surface concentration of the negative particle in comparison to the positive particle, given its larger solid diffusion time constant. The electrode-averaged electrolyte concentration is also modeled with great accuracy: taking into account that its typical concentration is 𝑐𝑒,𝑡𝑦𝑝 = 1000 mol/m3, the relative RMS and maximum errors are below 0.6% and 3% respectively. Note that the results for the electrolyte concentration are equal in both electrodes due to the fact that their thickness and porosity have the same values in this parameter set; this is not a general result nonetheless. In conclusion, the proposed ECM is able to accurately model the output voltage as well as the internal states with respect to an SPMe. Lastly, a comparison on the computational requirements of both models is carried out. We have determined that the proposed ECM is about 3 to 4 times quicker on average for the same timestep and simulation profile; however, we believe that the main advantage of our model with respect to the SPMe is the number of states required. For a typical discretization, such as the one considered here, 100 states must be stored and updated every timestep for the SPMe, whereas only 18 are necessary for the full ECM and 9 for the simplified version if the 7-RC approximation of the ZARC element is employed. Nevertheless, other continuous approximations by 5 [53] and 3 [49] RC networks have also been reported in recent literature, so should the available memory be tightly constrained, the number of necessary states may be potentially reduced to 14 or 10 for the full ECM and 7 or 5 for the simplified ECM, without a major loss of accuracy. This is a crucial feature given the limitations on memory and computation power in on-board BMSs. Furthermore, although some state observers have been proposed for the SPM and SPMe [15,42], employing an ECM makes it simpler to implement a Dual Fractional-Order Extended Kalman Filter for the concurrent estimation of state of charge and circuit parameters, as we presented in [56], with the advantage that the identified parameters are directly related to their physical counterparts in this model. 4. Experimental application In this section, the simplified version of the proposed equivalent circuit model is compared to experimental data by identifying the equivalent circuit parameters from an EIS test. Journal of Energy Storage 64 (2023) 107150 9 P. Rodríguez-Iturriaga et al. Fig. 4. Terminal voltage for the SPMe and the proposed ECM for a constant-current discharge at (a) 2C, (b) 1C, (c) C/2 and (d) C/5 respectively, followed by a 30-min rest. Curves have been plotted with respect to normalized time, i.e. time divided by the nominal discharge duration according to the C-rate. Fig. 5. Terminal voltage and input current for the SPMe and the proposed ECM for the US06 driving cycle.