Physics-Based Equivalent Circuit Model Motivated by the Doyle–Fuller–Newman Model

Abstract

This work introduces a sophisticated impedance-based equivalent circuit model of the electrochemical processes inside a lithium-ion battery cell. The influence on the electrical voltage response is derived and merged into a mathematical calculation framework describing all fundamental phenomena inside a battery. The parameters, whose sole influences on the electric behaviour cannot be separated at the cell level, are summarised to derive a model with purely electrical quantities. We significantly reduce the model order compared to a physicochemical model while ensuring a minimal approximation error. Utilising the findings from the model derivation, we develop a parameterisation procedure to separate the individual processes occurring in the battery and to support a hypothesis of the assignment to positive and negative electrodes based on several indicia. For this purpose, electrochemical impedance spectroscopy and correlation analysis are used to calculate the distribution of the time constants. The final parameterised model has physics-based parameter variations, which ensures that the simulation over broad ranges of temperatures and states of charge results in a reasonable voltage response. The model’s physical basis enables extrapolation beyond the measured operation area, and the model verification shows less than a 10 mV root mean square error over a wide range of operations.

1. Introduction

For optimising battery system design for a wide range of applications, it is crucial to model the system’s behaviour in response to electrical load changes. Being highly accurate and as close as possible to the real-world behaviour of battery cells is crucial for simulations of the charging and discharging behaviour using a battery model. Here, an extended electrical equivalent circuit model (ECM) explicitly provides a good trade-off between accuracy, ease of parameterisation, and simulation time [1]. This extended ECM additionally avoids the need to measure all electrochemical parameters in complex analytical experiments while at the same time assigning each equivalent circuit element to an electrochemical process.

Impedance-based analysis of ageing mechanisms regularly uses EIS data in combination with ECMs to gain insights by considering the fitted resistance components. An appropriate model typically has to be selected by an expert in the field [2]. The presented model is derived directly from the DFN model using certain simplifications described in this work. Some of these simplifications are also commonly used in approaches like the single-particle or extended single-particle models. The model can be an appropriate choice for cases that aim to interpret the physical meaning of parameter changes. Besides this, the model is a high-precision model for use cases like development environments for diagnostic algorithms or model-based system design approaches. Physics-based parameterisation allows for the extrapolation of behaviour into unmeasured operating conditions because the parameter variations follow physical constraints. This possibility to extrapolate is a benefit that a mere black-box model may not offer. There are different ECMs like Thévenin ECMs, which consist of OCVs, serial resistance, and a series of n RC circuits. Additionally, there are many approaches with more complex ECMs that usually aim to solve one single problem.

Thèvenin ECMs can represent battery-internal processes that rise or decay over time. These include charge-carrier depletion and recovery processes during battery charging and discharging [3]. Depending on the current applied to the battery, concentration gradients occur in the active material of the electrodes and the electrolyte during the charging and discharging processes. The reason for this is the components’ finite conductivity, leading to overpotentials in the positive electrode during charging and in the negative direction during discharging. When the battery terminals are disconnected and the battery is no longer under load, the overpotential relaxes as the inhomogeneities within the electrodes slowly equalise and concentration gradients are reduced. The potential built up across the double-layer capacitance also relaxes over time [4,5]. Due to the resulting overpotentials, the usable battery capacity decreases with increasing current rates [5]. However, the capacity due to these effects is not lost and can be recovered with sufficient relaxation time.

The chosen number of RC circuits in a Thévenin circuit influences its computational complexity [6]. For different values of n in RC circuits, the models are denoted as “RnRC”. The most popular models in the literature use two or three RC circuits due to their excellent trade-off between computational complexity and model accuracy. Four RC circuits provide only negligible improvement in model accuracy [7]. Other publications also show models with more than two RC circuits [7,8,9,10,11,12,13].

In some approaches in the literature, only one RC circuit is chosen [8,14,15,16,17,18,19,20,21,22]. This RC circuit must represent all the transient behaviour during charging and discharging. The main benefit of this model is the simplicity and ease of implementation [15].

A second RC circuit is used in several publications to separate the frequency-domain response into two parts [8,15,16,17,19,20,21,23,24,25,26,27,28,29,30,31,32,33,34,35]. Here, one RC circuit usually provides the overpotentials caused by both electrodes’ charge-transfer resistance and double-layer capacitance. The other represents the low-frequency behaviour, where diffusion effects are most dominant [12].

Starting with Thévenin models, various other models consisting of different ECEs have been proposed in recent publications. Most likely, in these publications, complex equivalent circuit elements are used to improve the models’ accuracy. Additionally, they try to incorporate further dependencies of ECMs’ parameters. All of the models have in common that either there is a voltage source for the OCV of the full cell or the anode’s and cathode’s open-circuit potentials are modelled separately, which results in two voltage sources. In the case of two voltage sources, the lithium plating process can be described.

In [36], the authors use a series connection of a series resistance and a ZARC element for the charge-transfer behaviour. Others [8,37,38,39] model the influence of the wiring and the current collectors with an inductor. They use a series resistance for the purely ohmic shift of the impedance and two ZARC elements: one for the impact of the solid electrolyte interface (SEI) and one for the charge-transfer behaviour and double-layer capacitance. A Warburg impedance models the diffusion processes that are dominant at low frequencies. The authors of [34,40] take an analogue approach but use a series connection of a resistance and an inductance for high frequencies.

In [41,42], the authors use an ECM consisting of a series resistance, an inductance, and a so-called Randles circuit [3], which is a parallel connection of a CPE and the series connection of a resistance and a Warbug impedance. To integrate a dependency of the electrical behaviour on the load (charging or discharging), the authors of [43,44] propose using ideal diodes in the ECM (so-called “splice-equivalent circuit model”). The authors of [42] add an additional ZARC element to model the impact of the SEI.

Some publications are working to close the gap between physicochemical models and ECMs. The authors of [45,46] describe an ECM implementation of the single-particle model, including electrolyte and thermal dynamics. In contrast to our approach, they aim to describe the separate influence of diffusion, electrolyte, and charge-transfer processes with a network of resistance, capacitance, and voltage source elements. The models are structurally different from commonly used ECMs, which consist of a voltage source resembling an OCV and a serial network of elements, each representing one process. The authors of [47] focus on the numerical approximation of the charge-transfer overpotential, solid diffusion, and electrolyte diffusion transfer functions, reducing the transfer functions to a series connection of RC circuits and ZARC elements to correlate the fitted resistance and time constants with physical parameters.

In this work, a physics-based ECM is used, which divides the battery’s behaviour into three parts depending on the frequency range in which the processes occur. Some of the processes occur only in low-frequency, medium-frequency, or high-frequency ranges, and the presented model can consequently be classified as a grey-box model. Since a simple RC element cannot represent the complex behaviour of all different electrochemical processes, the model consists of advanced equivalent circuit elements, such as ZARC or Warburg impedances.

The presented work provides a more detailed derivation of the impedances resulting from the dominant processes in the electrodes, separator, and electrolyte. During the derivation, commonly used simplifications known from, e.g., the extended single-particle model, are used. The derived model is purposefully not reduced to a series connection of commonly used equivalent circuit elements such as resistances, RC circuits, and ZARC elements before parameterisation. The approximation is only performed after parameterisation, leading to more physically meaningful parameter profiles. Different accuracies can be obtained depending on the selected time increment and by simplifying the model for simulation purposes. The final validation proves that the model and its parameters are valid over a wide operating range concerning the state of charge (SOC) and temperature.

2. Model Derivation

This section describes the derivation of the presented impedance-based model starting from the linearisation of the key processes contributing to the electrical behaviour of lithium-ion batteries. Changes in the electrode potentials are related to the net current to derive the impedance.

2.1. Intra-Particle Diffusion

The derivation of the impedance resulting from solid-state diffusion within particles of different geometries and boundary conditions is described here only for the spherical particles. The derivation for the planar particles follows an analogue approach, and the equations are given in the appendix.

Bard et al. describe Fick’s diffusion as a stochastic process that has its origin in Brownian particle motion [48]. The superposition of these random movements causes a particle flow in the opposite direction of the concentration gradients. This relationship can be expressed as

$$\overrightarrow{J}(\overrightarrow{x},t)=-D\nabla c(\overrightarrow{x},t)$$

(1)

and the continuity equation

$${\displaystyle \frac{\partial c(\overrightarrow{x},t)}{\partial t}}=-\nabla ·\overrightarrow{J}(\overrightarrow{x},t)$$

(2)

where J is the particle flux density of the considered particle, c is the concentration, and D is the diffusion coefficient. These equations can be combined to form Fick’s second law:

$${\displaystyle \frac{\partial c(\overrightarrow{x},t)}{\partial t}}=\nabla (D\nabla \left[c(\overrightarrow{x},t)\right])$$

(3)

Figure 1a shows a discretisation of a spherical particle into a spherical shell. $c_{\mathrm{i}}$, $r_{\mathrm{i}}$, and ${\underline{J}}_{\mathrm{i}}$ denote the lithium concentration, radius, and flux density at the i-th shell. ${\underline{J}}_{r_{\mathrm{P}}}$ and $c_{r_{\mathrm{P}}}$ denote the values at the outer particle radius $r_{\mathrm{P}}$. Figure 1b depicts a circuit representation of Figure 1a for a discretisation into n particles as a Cauer ladder. In the following section, a conversion to a Foster circuit is derived to find a compact formula for the impedance. Additionally, the simulation framework uses the Foster circuit and is therefore preferred here.

Equation (1) multiplied by the surface of the particle $A\left(r\right)$ at radius r can be expressed in spherical coordinates for a symmetrical spherical particle:

$$J(r,t)A\left(r\right)=-D·A\left(r\right){\displaystyle \frac{\partial c(r,t)}{\partial r}}=-D·4\pi r^2{\displaystyle \frac{\partial c(r,t)}{\partial r}}$$

(4)

where $J(r,t)$ is the flux density at radius r and time t. Rearranging this for the derivative of the lithium-ion concentration results in

$${\displaystyle \frac{\partial c(r,t)}{\partial r}}=-{\displaystyle \frac{J(r,t)·A\left(r\right)}{4\pi D}}{\displaystyle \frac{1}{r^2}}$$

(5)

Transferring the diffusion process into an electric equivalent circuit, we can assume that ${\underline{J}}_{r_{\mathrm{P}}}·A_{r_{\mathrm{P}}}$ is analogous to an electrical current driven by the concentration gradient; a resistance describes the diffusion path. The concentration difference between two volume elements resembles a voltage over the resistance, where the volume acts like a capacitance. The related equivalent circuit representation is shown in Figure 1b. Therefore, we assume a constant net flux between two spherical shells with radius $r_1$ and $r_2$:

$$c(r,t)={\displaystyle \frac{J(r,t)·A\left(r\right)}{4\pi D}}{\displaystyle \frac{1}{r}}+c_{\mathrm{const}}$$

(6)

Inserting the resulting boundary condition into Equation (6) leads to the approximation for the resistance element of the diffusion path from $r_1$ to $r_2$:

$${\displaystyle \frac{c_2\left(t\right)-c_1\left(t\right)}{J(r,t)A\left(r\right)}}={\displaystyle \frac{\frac{1}{r_2}-\frac{1}{r_1}}{4\pi D}}$$

(7)

The average radius $r_k$ of the k-th sphere shell of n segments of volume $V_k$ can be written as

$$r_k=r_{\mathrm{P}}\left(1-{\displaystyle \frac{k-1}{n}}\right)$$

(8)

To account for the diffusive losses within each spherical segment, their volumes are considered. The volume $V_k$ of each segment is

$$V_k={\displaystyle \frac{4}{3}}\pi r_{\mathrm{P}}^3\left[{\left(1-{\displaystyle \frac{k-1}{n}}\right)}^3-{\left(1-{\displaystyle \frac{k}{n}}\right)}^3\right]$$

(9)

The next step is the derivation of the transfer function between the particle flow and surface concentration:

$${\displaystyle \frac{{\underline{c}}_{r_{\mathrm{P}}}}{{\underline{J}}_{r_{\mathrm{P}}}·A}}=F\left(s\right)$$

(10)

where $F\left(s\right)$ is the Laplace transfer of the transfer function, which can only be derived with a fixed boundary condition at $r_{\min }$. Jacobsen and West, as well as Barsoukov et al., described the solution for $J_{r_{\min }}=0$ [49,50] as

$${\displaystyle \frac{{\underline{c}}_{r_{\mathrm{P}}}}{{\underline{J}}_{r_{\mathrm{P}}}·A}}={\displaystyle \frac{r_{\mathrm{P}}}{A·D}}·{\displaystyle \frac{tanh\left(\sqrt{j\omega ·\frac{r_{\mathrm{P}}^2}{D}}\right)}{\sqrt{j\omega ·\frac{r_{\mathrm{P}}^2}{D}}-tanh\left(\sqrt{j\omega ·\frac{r_{\mathrm{P}}^2}{D}}\right)}}$$

(11)

Jacobsen and West presented similar derivations for planar and cylindrical diffusion [49]. We assume the diffusivity D to be constant as a function of c. This assumption is a simplification and, therefore, introduces a slight inaccuracy to the model. The concentration at the particle surface gives the measurable potential of the battery:

$$\varphi \left(c_{r_{\mathrm{P}}}\right)\stackrel{!}{=}\varphi \left(c_0\right)+(c_{r_{\mathrm{P}}}-c_0){\displaystyle \frac{\partial \varphi }{\partial c}}$$

(12)

where $c_0$ is the equilibrium concentration and $\varphi \left(c_0\right)$ is the corresponding potential. The relation between the potential and the battery’s current can be described by assuming a piecewise linear relation between the concentration and the potential:

$${\displaystyle \frac{\underline{\varphi }}{{\underline{J}}_{r_{\mathrm{P}}}·A}}={\displaystyle \frac{\partial \varphi }{\partial c}}{\displaystyle \frac{r_{\mathrm{P}}}{A·D}}·{\displaystyle \frac{tanh\left(\sqrt{j\omega ·\frac{r_{\mathrm{P}}^2}{D}}\right)}{\sqrt{j\omega ·\frac{r_{\mathrm{P}}^2}{D}}-tanh\left(\sqrt{j\omega ·\frac{r_{\mathrm{P}}^2}{D}}\right)}}$$

(13)

where ${\displaystyle \frac{\partial \varphi }{\partial c}}\le 0$ during the lithiation of a particle since the potential of an electrode is always given in relation to a lithium reference electrode. The related impedance of a two-particle system can then be written as

$$\underline{Z}={\displaystyle \frac{{\varphi }^{+}}{{\underline{I}}_{\mathrm{Batt}}}}-{\displaystyle \frac{{\varphi }^{-}}{{\underline{I}}_{\mathrm{Batt}}}}={\displaystyle \frac{{\varphi }^{+}}{-{\underline{J}}_{r_{\mathrm{P}}}^{+}{A}^{+}z·F}}-{\displaystyle \frac{{\varphi }^{-}}{{\underline{J}}_{r_{\mathrm{P}}}^{-}{A}^{-}z·F}}={\underline{Z}}_{\mathrm{Diff}}^{+}+{\underline{Z}}_{\mathrm{Diff}}^{-}$$

(14)

where ${\underline{I}}_{\mathrm{Batt}}$ is the net charging current. A corresponding impedance can be derived for the spherical and planar diffusion processes. Making some assumptions to formulate that impedance analogous to a Voigt network is necessary. For planar diffusion, ${\underline{Z}}_{\mathrm{Particle}}^{\mathrm{cc}}$ and ${\underline{Z}}_{\mathrm{Particle}}^{\mathrm{cf}}$ can be analytically expressed by a series expansion. This analytical approach is not possible for spherical diffusion. Therefore, ${\underline{Z}}_{\mathrm{Particle}}^{\mathrm{sph}}$ is approximated using a finite series connection of RC circuits. Then, ${\tau }_{\mathrm{D}}^{\mathrm{sph}}$ and $R_{\mathrm{D}}^{\mathrm{sph}}$ are defined as physically reasonable time constants for the diffusion [51].

The impedance ${\underline{Z}}_{\mathrm{Particle}}^{\mathrm{sph}}$ of a single particle can therefore be written as

$${\underline{Z}}_{\mathrm{Particle}}^{\mathrm{sph}}={\displaystyle \frac{-\partial \varphi }{\partial c}}{\displaystyle \frac{r_{\mathrm{P}}}{A·D·z·F}}·{\displaystyle \frac{tanh\left(\sqrt{j\omega ·\frac{r_{\mathrm{P}}^2}{D}}\right)}{\sqrt{j\omega ·\frac{r_{\mathrm{P}}^2}{D}}-tanh\left(\sqrt{j\omega ·\frac{r_{\mathrm{P}}^2}{D}}\right)}}$$

(15)

Assuming a parallel connection of all particles in an electrode, the associated resistance of the diffusion process can simply be divided by the number of particles. The impedance of diffusion processes in an electrode can then be defined as

$${\underline{Z}}_{\mathrm{Diff}}^{\mathrm{sph}}\stackrel{!}{=}10·R_{\mathrm{D}}·{\displaystyle \frac{tanh\left(\sqrt{j\omega ·{\tau }_{\mathrm{D}}^{\mathrm{sph}}·20}\right)}{\sqrt{j\omega ·{\tau }_{\mathrm{D}}^{\mathrm{sph}}·20}-tanh\left(\sqrt{j\omega {\tau }_{\mathrm{D}}^{\mathrm{sph}}·20}\right)}}$$

(16)

where

$${\tau }_{\mathrm{D}}^{\mathrm{sph}}={\displaystyle \frac{r_{\mathrm{P}}^2}{20·D}}$$

(17)

$$R_{\mathrm{D}}^{\mathrm{sph}}={\displaystyle \frac{1}{10}}{\displaystyle \frac{-\partial \varphi }{\partial c}}{\displaystyle \frac{r_{\mathrm{P}}}{A_{\mathrm{act}},\mathrm{electrode}·D·z·F}}={\displaystyle \frac{1}{10}}{\displaystyle \frac{-\partial \varphi }{\partial c}}{\displaystyle \frac{1}{4\pi r_{\mathrm{P}}·n_{\mathrm{Particle}}·D·z·F}}$$

(18)

For planar diffusion, the derivation using the equations in the appendix given the boundary conditions of the constant concentration ${\underline{Z}}_{\mathrm{Particle}}^{\mathrm{cc}}$ and constant flux ${\underline{Z}}_{\mathrm{Particle}}^{\mathrm{cf}}$, respectively, leads to

$${\underline{Z}}_{\mathrm{Particle}}^{\mathrm{cc}}={\displaystyle \frac{-\partial \varphi }{\partial c}}{\displaystyle \frac{l}{ADzF}}{\displaystyle \frac{tanh\left(\sqrt{j\omega ·\frac{lV}{AD}}\right)}{\sqrt{j\omega ·\frac{lV}{AD}}}}={\displaystyle \frac{-\partial \varphi }{\partial c}}{\displaystyle \frac{l}{ADzF}}{\displaystyle \frac{tanh\left(\sqrt{j\omega ·\frac{l^2}{D}}\right)}{\sqrt{j\omega ·\frac{l^2}{D}}}}$$

(19)

$${\underline{Z}}_{\mathrm{Particle}}^{\mathrm{cf}}={\displaystyle \frac{-\partial \varphi }{\partial c}}{\displaystyle \frac{l}{ADzF}}{\displaystyle \frac{coth\left(\sqrt{j\omega ·\frac{lV}{AD}}\right)}{\sqrt{j\omega ·\frac{lV}{AD}}}}={\displaystyle \frac{-\partial \varphi }{\partial c}}{\displaystyle \frac{l}{ADzF}}{\displaystyle \frac{coth\left(\sqrt{j\omega ·\frac{l^2}{D}}\right)}{\sqrt{j\omega ·\frac{l^2}{D}}}}$$

(20)

The impedance of the planar diffusion process can be derived as

$${\underline{Z}}_{\mathrm{Diff}}^{\mathrm{cc}}\stackrel{!}{=}{\displaystyle \frac{{\pi }^2}{8}}{R}_{\mathrm{D}}^{\mathrm{cc}}{\displaystyle \frac{tanh\left(\sqrt{j\omega \frac{{\pi }^2}{4}{\tau }_{\mathrm{D}}^{\mathrm{cc}}}\right)}{\left.\sqrt{j\omega \frac{{\pi }^2}{4}{\tau }_{\mathrm{D}}^{\mathrm{cf}}}\right)}}$$

(21)

$$R_{\mathrm{D}}^{\mathrm{cc}}={\displaystyle \frac{8}{{\pi }^2}}{\displaystyle \frac{-\partial \varphi }{\partial c}}{\displaystyle \frac{l}{ADzF}},{\tau }_{\mathrm{D}}^{\mathrm{cc}}={\displaystyle \frac{4l^2}{{\pi }^{2}D}}$$

(22)

$${\underline{Z}}_{\mathrm{Diff}}^{\mathrm{cf}}\stackrel{!}{=}{\displaystyle \frac{{\pi }^2}{2}}{R}_{\mathrm{D}}^{\mathrm{cf}}{\displaystyle \frac{coth\left(\sqrt{j\omega {\pi }^2{\tau }_{\mathrm{D}}^{\mathrm{cf}}}\right)}{\sqrt{j\omega {\pi }^2{\tau }_{\mathrm{D}}^{\mathrm{cf}}}}}$$

(23)

$$R_{\mathrm{D}}^{\mathrm{cf}}={\displaystyle \frac{2}{{\pi }^2}}{\displaystyle \frac{-\partial \varphi }{\partial c}}{\displaystyle \frac{l}{ADzF}},{\tau }_{\mathrm{D}}^{\mathrm{cf}}={\displaystyle \frac{l^2}{{\pi }^{2}D}}$$

(24)

For the simulation of the model in the time domain, the described impedance has to be approximated into a series of RC circuits. The associated methodology can be considered in the published implementation of the model (https://git.rwth-aachen.de/isea/framework, https://git.rwth-aachen.de/isea/ecmfitting, accessed on 14 July 2024).

2.2. Electrode–Electrolyte Interface

If a current is applied to a lithium-ion battery, an equal ionic current forms inside the cell under the assumption that no side reactions are taking place [52]. Lithium ions intercalate and de-intercalate from electrode particles, causing overpotentials. Randles derived an equivalent circuit model resembling the impedance resulting from the charge-transfer overpotentials, the double-layer capacitance, and the intra-particle diffusion, as described in Section 2.1 [53].

The overpotentials arising from the charge-transfer reactions are modelled as a current-dependent resistance $R_{\mathrm{ct}}$. The potential difference across the electrolyte causes the negative and positive ions in the electrolyte to orient themselves to the electrode–electrolyte interfaces. This effect can be modelled as a capacitance $C_{\mathrm{dl}}$. The resulting time constant of the combination of charge-transfer resistance and double-layer capacitance is ${\tau }_{\mathrm{dl}}=R_{\mathrm{ct}}·C_{\mathrm{dl}}$. The circuit representation of a Randles circuit is shown in Figure 2a. The associated impedance is

$${\underline{Z}}_{\mathrm{Randles}}={\left({\left(R_{\mathrm{ct}}+{\underline{Z}}_{\mathrm{Diff}}\right)}^{-1}+{\left({\displaystyle \frac{1}{j\omega C_{\mathrm{dl}}}}\right)}^{-1}\right)}^{-1}$$

(25)

2.3. Porous Electrode

In addition to the reaction kinetics and mass transport within particles described in Section 2.1 and Section 2.2, the transport processes of electrons to the reaction interface and lithium ions within the electrolyte and the porous electrode have to be considered.

2.3.1. Ionic and Electric Conductivities

The electrical contact is primarily realised by the current collectors, which have high conductivities compared to the battery’s other components. In the negative electrode, electron conduction is mainly carried out by the conductive graphite and conductive carbon black. In the positive electrode, electron conduction is mainly facilitated by the conductive carbon black. In all cases, only migration caused by electric fields takes place. Overpotentials due to electron conduction can be modelled by an ohmic resistance $R_{\mathrm{ser}}$ due to the lack of concentration gradients.

Ionic conduction in the electrolyte comprises both migration and diffusion, as concentration gradients arise during transport. Both lithium ions and negatively charged ions, e.g., ${\mathrm{PF}}_6^{-}$-anions, participate in the migration current. The flux across the electrode–electrolyte interface can only be realised by lithium cations, causing concentration gradients within the electrolyte. The resulting diffusion current counteracts the polarisation. In an equilibrium state, a constant concentration gradient arises, and the diffusion current of the anions nullifies the anion part of the migration current.

The effective ionic conductivity is defined as ${\kappa }_{\mathrm{eff}}=\kappa ·{\epsilon }^p$, where p is the Bruggeman factor and $\epsilon$ is the porosity. Newman describes the migration overpotential ${\varphi }_{\mathrm{M}}$ in an electrolyte-filled porous structure of length l, surface A, and ${\kappa }_{\mathrm{eff}}$ [54] as

$${\varphi }_{\mathrm{M}}=I·{\displaystyle \frac{l}{A·{\kappa }_{\mathrm{eff}}}}$$

(26)

Therefore, the associated resistance is, e.g., the migration resistance in the separator $R_{\mathrm{M},\mathrm{Sep}}$ and the porous electrode $R_{\mathrm{M},\mathrm{P}}$

$$R_{\mathrm{M},\mathrm{X}}={\displaystyle \frac{l}{A·{\kappa }_{\mathrm{eff}}}}$$

(27)

2.3.2. Electrolyte Overpotential

The potential of an electrode in an electrolyte can be described using the Nernst equation [55]:

$$\varphi ={\varphi }_0+{\displaystyle \frac{RT}{zF}}·ln\left({\displaystyle \frac{a_{{\mathrm{Li}}^{+}}}{a_{\mathrm{Li}}}}\right)$$

(28)

where ${\varphi }_0$ is the potential of the considered metal, $a_{{\mathrm{Li}}^{+}}$ is the electrochemical activity of lithium ions in the electrolyte, and $a_{\mathrm{Li}}$ is the activity of the lithium ions on the electrode surface. The electrochemical activity is directly proportional to the concentration via the activity coefficient $\gamma$:

$$a=\gamma ·c$$

(29)

The equilibrium voltage of an electrode ${\varphi }_{\mathrm{eq}}$ is

$${\varphi }_{\mathrm{eq}}={\varphi }_0+{\displaystyle \frac{RT}{zF}}·ln\left({\displaystyle \frac{{\gamma }_{{\mathrm{Li}}^{+}}·c_{\mathrm{eq},{\mathrm{Li}}^{+}}}{{\gamma }_{\mathrm{Li}}·c_{\mathrm{eq},\mathrm{Li}}}}\right)$$

(30)

The potential of an electrode given a relative change in concentration ${\displaystyle \frac{c_{{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}$ is

$$\begin{array}{cc}\hfill \varphi & ={\varphi }_0+{\displaystyle \frac{RT}{zF}}·ln\left({\displaystyle \frac{{\gamma }_{{\mathrm{Li}}^{+}}·c_{\mathrm{eq},{\mathrm{Li}}^{+}}·\frac{c_{{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}{{\gamma }_{\mathrm{Li}}·c_{\mathrm{eq},\mathrm{Li}}}}\right)\hfill \\ \hfill & ={\varphi }_{\mathrm{eq}}+{\displaystyle \frac{RT}{zF}}·ln\left({\displaystyle \frac{c_{{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)\hfill \end{array}$$

(31)

The change in the local surface potential of an electrode is

$$\mathsf{\Delta }{\varphi }_{\mathrm{Nernst}}={\displaystyle \frac{RT}{zF}}·ln\left({\displaystyle \frac{c_{{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)$$

(32)

The measurable change in potential at the battery tabs, therefore, is

$$\mathsf{\Delta }{\varphi }_{\mathrm{Nernst}}={\displaystyle \frac{RT}{zF}}·ln\left({\displaystyle \frac{c_{{\mathrm{Li}}^{+}}^{\mathrm{pos}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)+{\displaystyle \frac{RT}{zF}}·ln\left({\displaystyle \frac{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}{c_{{\mathrm{Li}}^{+}}^{\mathrm{neg}}}}\right)$$

(33)

Additionally, Hamann and Vetter described an additional voltage change called the fluid diffusion overpotential [55,56]:

$$\mathsf{\Delta }{\varphi }_{\mathrm{Fl},\mathrm{Diff}}=-{\displaystyle \frac{RT}{zF}}\left(t^{+}-t^{-}\right)·ln\left({\displaystyle \frac{c_{{\mathrm{Li}}^{+}}^{\mathrm{pos}}}{c_{{\mathrm{Li}}^{+}}^{\mathrm{neg}}}}\right)={\displaystyle \frac{RT}{zF}}\left(1-2t^{+}\right)·ln\left({\displaystyle \frac{c_{{\mathrm{Li}}^{+}}^{\mathrm{pos}}}{c_{{\mathrm{Li}}^{+}}^{\mathrm{neg}}}}\right)$$

(34)

For two adjacent volumes k and $k+1$, this can be expressed as

$$\mathsf{\Delta }{\varphi }_{k,(k+1),\mathrm{Fl},\mathrm{Diff}}={\displaystyle \frac{RT}{zF}}·\left(1-2t^{+}\right)·ln\left({\displaystyle \frac{c_{k{\mathrm{Li}}^{+}}}{c_{(k+1){\mathrm{Li}}^{+}}}}\right)$$

(35)

Separator

The fluid diffusion potential can be modelled as a Warburg impedance of planar diffusion with a constant concentration in a distance $l/2$, where l is the thickness of the separator (Equation (21)). The derivative of the sum of Equations (35) and (32) is

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \varphi }{\partial c_{{\mathrm{Li}}^{+}}}}& ={\displaystyle \frac{\partial }{\partial c_{{\mathrm{Li}}^{+}}}}\left[{\displaystyle \frac{RT}{zF}}·ln\left({\displaystyle \frac{c_{{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)+{\displaystyle \frac{RT}{zF}}·\left(1-2t^{+}\right)·ln\left({\displaystyle \frac{c_{{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)\right]\hfill \\ \hfill & ={\displaystyle \frac{\partial }{\partial c_{{\mathrm{Li}}^{+}}}}\left[{\displaystyle \frac{RT}{zF}}·2t^{-}·ln\left({\displaystyle \frac{c_{{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)\right]\hfill \\ \hfill & ={\displaystyle \frac{2t^{-}RT}{zF{c}_{{\mathrm{Li}}^{+}}}}\hfill \end{array}$$

(36)

Assuming a symmetrical separator, the impedance in the second half of the separator is equal due to concentration changes of similar sizes and opposing signs. The source term is

$$J·A={\displaystyle \frac{t^{-}·I_{\mathrm{Batt}}}{z·F}}$$

(37)

For a constant flux, the difference in concentration across a volume element in the separator, neglecting any migration, is

$${\displaystyle \frac{\mathsf{\Delta }c}{JA}}={\displaystyle \frac{l}{nA{D}_{\mathrm{eff}}}}$$

(38)

The associated resistance of the planar diffusion model can be calculated as

$$R_{\mathrm{D},\mathrm{Sep}}={\displaystyle \frac{\mathsf{\Delta }U}{I_{\mathrm{Batt}}}}={\displaystyle \frac{\frac{\partial \varphi }{\partial c_{{\mathrm{Li}}^{+}}}·\mathsf{\Delta }c}{JA\frac{zF}{t^{-}}}}={\displaystyle \frac{2t^{-}RT}{zF{c}_{\mathrm{eq},{\mathrm{Li}}^{+}}}}·{\displaystyle \frac{t^{-}·l}{A{D}_{\mathrm{eff}}zF}}={\displaystyle \frac{2{\left(t^{-}\right)}^{2}RTl}{z^2{F}^{2}A{D}_{\mathrm{eff}}{c}_{\mathrm{eq},{\mathrm{Li}}^{+}}}}$$

(39)

Inserting both the effective volume $V_{\mathrm{eff}}=\epsilon V$ and effective diffusion $D_{\mathrm{eff}}=D·{\epsilon }^p$ into Equation (19) results in an adjusted time constant ${\tau }_{\mathrm{D}}$:

$${\tau }_{\mathrm{D},\mathrm{Sep}}={\displaystyle \frac{4\epsilon l^2}{{\pi }^2{D}_{\mathrm{eff}}}}$$

(40)

Therefore, the capacitance is

$$C_{\mathrm{D},\mathrm{Sep}}=\epsilon V·{\displaystyle \frac{z^2{F}^2{c}_{\mathrm{eq},{\mathrm{Li}}^{+}}}{2{\left(t^{-}\right)}^{2}RT}}$$

(41)

which results in the following impedance:

$${\underline{Z}}_{\mathrm{D},\mathrm{Sep}}={\displaystyle \frac{2{\left(t^{-}\right)}^{2}RTl}{z^2{F}^{2}A{D}_{\mathrm{eff}}{c}_{{\mathrm{Li}}^{+}}}}·{\displaystyle \frac{\tanh \left(\sqrt{j\omega ·\frac{\epsilon l^2}{D_{\mathrm{eff}}}}\right)}{\sqrt{j\omega ·\frac{\epsilon l^2}{D_{\mathrm{eff}}}}}}$$

(42)

The impedance related to concentration gradients in the electrolyte is generally small compared to the impedance related to concentration gradients in the electrode due to its smaller thickness. Therefore, it is neglected.

Porous Electrode

Similar equations are valid for both the separator and the electrodes. Käbitz [57] found that the overpotentials due to electrolyte concentration gradients within the porous electrode have a significant impact and should not be neglected. Using a parameter study, he also found that the overpotentials display a sharp time constant and can thus be modelled by a single RC circuit.

The relation between the electrolyte concentration and the respective overpotential has to be linearised to derive an analytical formulation for the impedance of this process. The linearisation $ln\left({\displaystyle \frac{x}{x_0}}\right)\approx (x-x_0)$ is valid for a reasonably large range of the amperage:

$$\mathsf{\Delta }{\varphi }_{k,\mathrm{Nernst}}={\displaystyle \frac{RT}{zF}}·ln\left({\displaystyle \frac{c_{k,{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)\approx {\displaystyle \frac{RT}{zF}}·\left({\displaystyle \frac{c_{k,{\mathrm{Li}}^{+}}-c_{\mathrm{eq},{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)$$

(43)

$$\begin{array}{cc}\hfill \mathsf{\Delta }{\varphi }_{k,(k+1),\mathrm{Fl},\mathrm{Diff}}& ={\displaystyle \frac{RT}{zF}}·\left(1-2t^{+}\right)·ln\left({\displaystyle \frac{c_{k,{\mathrm{Li}}^{+}}}{c_{(k+1),{\mathrm{Li}}^{+}}}}\right)\hfill \\ \hfill & ={\displaystyle \frac{RT}{zF}}·\left(1-2t^{+}\right)·ln\left({\displaystyle \frac{c_{k,{\mathrm{Li}}^{+}}·c_{\mathrm{eq},{\mathrm{Li}}^{+}}}{c_{(k+1),{\mathrm{Li}}^{+}}·c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)\hfill \\ \hfill & ={\displaystyle \frac{RT}{zF}}·\left(1-2t^{+}\right)·\left(ln\left({\displaystyle \frac{c_{k,{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)-ln\left({\displaystyle \frac{c_{(k+1),{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)\right)\hfill \\ \hfill & \approx {\displaystyle \frac{RT}{zF}}·\left(1-2t^{+}\right)·\left({\displaystyle \frac{c_{k,{\mathrm{Li}}^{+}}-c_{(k+1),{\mathrm{Li}}^{+}}}{c_{\mathrm{eq},{\mathrm{Li}}^{+}}}}\right)\hfill \end{array}$$

(44)

The parameters of the associated RC circuit can be calculated analogously to Equations (37) and (40):

$$\begin{array}{cc}\hfill R_{\mathrm{D},\mathrm{P}}& ={\displaystyle \frac{\mathsf{\Delta }{\varphi }_{k,\mathrm{Nernst}}}{I_{\mathrm{Batt}}}}={\displaystyle \frac{\frac{RT}{zF{c}_{\mathrm{eq},{\mathrm{Li}}^{+}}}·\mathsf{\Delta }c}{\frac{JAzF}{t^{-}}}}\hfill \\ \hfill & ={\displaystyle \frac{RT}{z^2{F}^2{c}_{\mathrm{eq},{\mathrm{Li}}^{+}}}}·\mathsf{\Delta }c·{\displaystyle \frac{t^{-}}{JA}}={\displaystyle \frac{RT{t}^{-}}{z^2{F}^2{c}_{\mathrm{eq},{\mathrm{Li}}^{+}}}}·{\displaystyle \frac{l}{A{D}_{\mathrm{eff}}}}\hfill \end{array}$$

(45)

Again, inserting effective quantities into the time constant for planar diffusion with a constant flux at the boundary and solving for the capacitance provides

$$C_{\mathrm{D},\mathrm{P}}={\displaystyle \frac{z^2{F}^2{c}_{\mathrm{eq},{\mathrm{Li}}^{+}}}{RT{t}^{-}}}·\epsilon V$$

(46)

2.3.3. Current Distribution

When a constant current is applied, a concentration distribution is established across the electrodes within a few seconds. A discharge current causes pores within the anode close to the separator to have a lower salt concentration than those close to the current collectors. This affects both the equilibrium potential of the electrode segments and the reaction kinetics. Both effects, in turn, influence the distribution of the current density across the electrode thickness.

If a concentration gradient is established during discharge, the particles near the separator are increasingly involved in the discharge. An equilibrium is formed, which accounts for the potentials due to changes in concentration, the reaction overpotentials, and the degree of loading of the individual particles.

The circuit diagram in Figure 2b can be used to calculate the equilibrium that develops given a concentration gradient within a porous electrode. Strictly speaking, an electrode cannot be considered individually, as the entire electrolyte volume has to be taken into account in the two-electrode setup. A similar structure in both electrodes is assumed, allowing for the depicted symmetry. Barsoukov [58] found that the overpotentials in the electrolyte can be approximated as one process if the electrolyte overpotentials within the separator are small compared to those in the porous electrode.

The resulting circuit can be solved analytically for an arbitrarily fine discretisation (see [51]), but the resulting computational complexity is impractical for an ECM. Witzenhausen chose an approximation neglecting the Nernst overpotentials and fluid diffusion overpotentials due to their minor influence [51]

$${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}=\sqrt{R_{\mathrm{M},\mathrm{P}}·{\underline{Z}}_{\mathrm{Randles}}}·coth\left(\sqrt{{\displaystyle \frac{R_{\mathrm{M},\mathrm{P}}}{{\underline{Z}}_{\mathrm{Randles}}}}}\right)$$

(47)

After subtracting ${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}^{+}$, ${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}^{-}$, and ${\underline{Z}}_{\mathrm{Sep}}$, only the previously neglected electrolyte overpotentials remain. Their shape in the impedance spectrum can be approximated with an RC circuit:

$${\underline{Z}}_{\mathrm{RC},\mathrm{Por},\mathrm{Diff}}={\displaystyle \frac{R_{\mathrm{D},\mathrm{P},\mathrm{eff}}}{1+j\omega {\tau }_{\mathrm{D},\mathrm{P}}}}$$

(48)

2.3.4. Particle Size Distribution

The influence on the electrical behaviour of the distribution of particle properties, such as the particle radius, is closely related to the lithium diffusion inside the particles. The voltage behaviour represents a superposition of the dynamic and capacitive elements of the diffusion impedance and the balancing processes between the different particles. This superposition creates a fundamental problem for parameter identification concerning diffusion and particle distributions based solely on the battery voltage. Both effects cannot be separated from each other since they merge smoothly. Similar to what was explained in Section 2.3.3 for the porous structure, the current is also distributed non-homogeneously to different active mass particles. At first, the current distribution is almost exclusively proportional to the active surface of the particles, as described in the model of Meyers [59]. As a result, smaller particles are loaded more heavily in relation to their volume than larger ones, and the change in the SOC is greater. This, in turn, induces a charge mismatch between the particles.

Meyers developed a complex model framework for describing this interplay, based on the Newman model equations developed in his research group [59]. From an electrotechnical point of view, the different particles with their charge-dependent surface potentials are connected in parallel, and the current is distributed according to Kirchhoff’s rules. To calculate the small-signal impedance of a particle distribution, it is sufficient to connect the equivalent circuit elements developed up to this point for individual particles and their surface potentials in parallel to calculate the resulting impedance.

To account for interactions due to the distribution of particle radii and interactions between the transfer reaction and diffusion when the respective time constants are of similar size, the ZARC element in the Randles circuit is replaced with an HN element, named after the researchers Havriliak and Negami [60]:

$${\underline{Z}}_{\mathrm{HN}}={\displaystyle \frac{R}{{\left(1+{\left(j\omega \tau \right)}^{\theta }\right)}^{\phi }}}$$

(49)

Using the HN element, the circuit diagram accounts for cases with asymmetrical distributions of time constants.

In order to further improve the model’s accuracy, an additional RC circuit with the impedance ${\underline{Z}}_{\mathrm{LF}}$ is added. This RC circuit accounts for other slow processes not already covered by the presented equation system. It can be calculated using

$${\underline{Z}}_{\mathrm{LF}}={\displaystyle \frac{R_{\mathrm{LF}}}{1+j\omega {\tau }_{\mathrm{LF}}}}$$

(50)

where $R_{\mathrm{LF}}$ and ${\tau }_{\mathrm{LF}}=R_{\mathrm{LF}}·C_{\mathrm{LF}}$ represent the resistance and time constant of these slow processes.

The resulting series connection of all derived equivalent circuit elements is presented in Figure 3. The capacitances in the ECM are correlated with their time constants in the specific equations describing the equivalent circuit elements. A collection of all parameters to be identified is shown in

Table 1. ECM parameters and their dependencies.

Parameter	Part of	Equation	Physical Effect
$V_{\mathrm{OCV}}$	-	-	Open-circuit voltage
$R_{\mathrm{ser}}$	-	-	Electron conduction
$R_{\mathrm{ct}}^{+}$	${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}^{+}$	Equation (25)	Electrode–electrolyte
$C_{\mathrm{dl}}^{+}$	${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}^{+}$	Equation (25)	interface
${\underline{Z}}_{\mathrm{Diff}}^{+}$	${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}^{+}$	Equation (25)
$R_{\mathrm{M},\mathrm{P}}^{+}$	${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}^{+/-}$	Equation (47)
$R_{\mathrm{M},\mathrm{Sep}}$	${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}^{+/-}$	Equation (47)
$R_{\mathrm{M},\mathrm{P}}^{-}$	${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}^{+/-}$	Equation (47)
$R_{\mathrm{D},\mathrm{P}}$	${\underline{Z}}_{\mathrm{RC},\mathrm{Por},\mathrm{Diff}}$	Equation (48)	Porous electrode
$C_{\mathrm{D},\mathrm{P}}$	${\underline{Z}}_{\mathrm{RC},\mathrm{Por},\mathrm{Diff}}$	Equation (48)
${\underline{Z}}_{\mathrm{Diff}}^{+}$	${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}^{-}$	Equation (25)	Electrode–electrolyte
$R_{\mathrm{ct}}^{-}$	${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}^{-}$	Equation (25)	interface
$C_{\mathrm{dl}}^{-}$	${\underline{Z}}_{\mathrm{Por},\mathrm{El},\mathrm{M}}^{-}$	Equation (25)
$R_{\mathrm{LF}}$	${\underline{Z}}_{\mathrm{LF}}$	Equation (50)	Equalisation processes
$C_{\mathrm{LF}}$	${\underline{Z}}_{\mathrm{LF}}$	Equation (50)

.

3. Validation

For the validation, we present the measurements used to extract the raw data for our parameter identification process. As a result of this process, the extracted parameters exhibit steady behaviour over the SOC and temperature, and simulations using these parameters show promising results for verifying the model.

3.1. Measurements

We conducted measurements at seven different temperatures for the holistic characterisation of the battery. The operating temperature ranged from −15 °C to +35 °C in 10 K steps. The measurement device recorded current, voltage, and temperature values at every temperature for the same measurement profile, as shown in Figure 4a.

For the measurements, the following equipment was used. The current load for the battery was charged and discharged using a cycler of the MCT 100-06-12 ME type by Digatron, which has an accuracy of 0.1% of the end value. It was calibrated before the test to 0.14 mV and 0.41 mA for voltage and current measurements, respectively. A temperature chamber of the MK 240 type by Binder ensured a correct constant ambient temperature, with an accuracy of 0.1 K to 0.5 K depending on the target temperature, and 0.1 K to 1.2 K over time and space, respectively. We performed electrochemical impedance spectroscopy (EIS) using an EISMeter device by Digatron. A temperature sensor attached to the surface of the prismatic battery housing, also by Digatron, measured the cell temperature.

The measurement procedure and the model’s resulting parameterisation are described for a prismatic automotive battery cell of the first-generation Mitsubishi iMiEV, labelled LEV50 by GS Yuasa. According to the datasheet, this cell can be discharged with a maximum current of 300 A. The capacity at the time of parameterisation, measured with a discharge current of 50 A (1 C), was 48 Ah. The end-of-charge and end-of-discharge voltages for a constant current (CC) phase were defined as 4.2 V and 2.5 V, respectively. The charging process’s constant voltage (CV) phase terminated if the charging current decreased to less than C/20.

We fully charged the battery cell at 25 °C and performed a capacity test. Here, the capacity of the battery cell is the amount of charge discharged. A resting phase followed, during which the environmental temperature changed from 25 °C to the specified test temperature. This process took several hours and ended when the cell’s target temperature was stable for more than one hour.

At the target temperature, the characterisation procedures were carried out, including a relaxation phase, EIS, and pulse tests (see Figure 4b,c), performed at different SOCs in 5% steps. For the initial relaxation phase at every SOC, which began after the discharge ended, a $\mathsf{\Delta }U$ criterion was applied. Experience from relaxation tests at the institute indicated that the phase ended when the voltage changed by less than 1 mV in 1 h. To ensure a completed relaxation process and a valid criterion for a linear time-invariant system [61], the EIS measurement took place after the criterion was met [62]. The EIS excited the battery cell with a sinusoidal current at various frequencies from 6 kHz down to 1 mHz, and the division of voltage and current resulted in the frequency-dependent impedance. The measured frequencies were distributed logarithmically over the range, as shown in Figure 4d. The maximum amplitude of the sinusoidal excitation was 2 A, and the required amplitude of the voltage response was 3 mV. Figure 4b shows the results for the EIS measurements at 50% SOC and 15 °C in the Nyquist plot.

Afterwards, we performed pulse tests. For this purpose, we charged and discharged the battery cell with current rates of 1 C and 2 C for a duration of 20 s. At SOCs of 100% and 0%, the charge or discharge pulses were not performed due to the battery cell’s voltage limits. After every pulse and the subsequent relaxation phase, we reset the SOC by charging or discharging to the previous value with 0.5 C. This reset ensured that all pulses started from the same SOC value. There was a relaxation phase before each pulse to ensure that all cell-internal electrochemical processes had stabilised. For this, a slightly different $\mathsf{\Delta }U$ criterion had to be met, where the voltage decrease within 15 min had to be smaller than 10 mV. The last pulse represented the end of the measurements at the specific SOC level, and the cell was discharged by a further 5% SOC.

The characterisation of one temperature finished with a quasi-open-circuit voltage (qOCV) measurement, where the battery cell was discharged with a minor current of C/20. The results could be used for the open-circuit voltage (OCV) element of the ECM and provided a good trade-off between measurement accuracy and capturing all dynamics of the OCV curve.

3.2. Model Parameters

After measuring all the supporting points of the test matrix regarding the SOC and temperature, we performed the parameter identification process. For the implementation of the process, we provide the source in the git repository. Starting with the EIS measurement data, a subsequent parameter fitting, which reduces the residual between measurement and calculation, was performed for parameters with fast to slow time constants. We started with $R_{\mathrm{ser}}$, usually defined as $\Re \left\{\underline{Z}\right\}(\Im \left\{\underline{Z}\right\}=0)$. After that, using the Distribution of Relaxation Times (DRT) [63] method, we determined the time constants for the double-layer capacitances ${\tau }_{\mathrm{dl}}^{-}$ and ${\tau }_{\mathrm{dl}}^{+}$ at the anode and cathode. The parameter identification of the other parameters of both Randles circuits followed by fitting them to the impedance spectrum. Subsequently, we transferred the fitted ECM from the frequency to the time domain. We identified the parameters for $R_{\mathrm{D},\mathrm{P}}$, ${\tau }_{\mathrm{D},\mathrm{P}}$, $R_{\mathrm{LF}}$, and ${\tau }_{\mathrm{LF}}$ by minimising the error of the ECM calculation to the initial relaxation measurement at every SOC before the EIS measurement was conducted. After finishing all local parameter fittings, manual post-processing to address outliers and physically unreasonable behaviour was performed.

The parameter identification process resulted in a look-up table for each parameter consisting of values for all $SOC$-T combinations. For clarity, the parameter variations are shown in SOC dependency for all temperatures, but the temperature dependency is only shown for four SOCs. As the values for the elements describing the slower processes become very high at −15 °C, in the figures showing the profiles of $R_{\mathrm{D}}^{+}$, ${\tau }_{\mathrm{D}}^{+}$, $R_{\mathrm{D},\mathrm{P}}$, ${\tau }_{\mathrm{D},\mathrm{P}}$, $R_{\mathrm{slow}}$, and ${\tau }_{\mathrm{slow}}$, this temperature is neglected.

3.2.1. Ohmic Resistance

The model’s series resistance $R_{\mathrm{ser}}$ includes all overpotentials that occur instantaneously after changing the current. These include all electrical connection resistances, electrical conductors, and ionic migration resistances, such as those at the solid electrolyte interface (SEI). A further separation of the individual causes is impossible without disassembling and analysing the battery cell. The series resistance $R_{\mathrm{ser}}$ showed hardly any dependence on the SOC and increased slightly with decreasing temperature. Despite the temperature range of 50 K, the resistance varied only by half a decade (see Figure 5a). This slight variance indicates that this fastest possible measurement does not include any electrochemical components. Experience has shown that the losses from these reactions would increase much more with decreasing temperature. Even if we look at the change in $R_{\mathrm{ser}}$ in the logarithmic representation (Figure 5d), it should not be assumed that $R_{\mathrm{ser}}$ has an Arrhenius relation because the change of only half a decade is too small.

3.2.2. Reaction of the Negative Electrode

First, we assumed ${\phi }_{\mathrm{HN}}^{-}=1$ over the entire measurement range, which means that the behaviour of the HN element was like a ZARC element (see Equation (49)). The underlying assumption was that the interaction between the charge transfer and the diffusion process would result in normally distributed time constants. We neglected the influences on the high-frequency uncertainty of the time constant due to a distribution of the particles’ binding since they were not in the measurable range, especially for the faster of the two reaction processes. The current dependencies of both $R_{\mathrm{ct}}$ elements were neglected due to the fact that the investigated high-energy battery was not capable of high currents.

The faster of the two electrode reaction processes was assigned to the anode. An electrode’s double-layer capacitance decreased during the SEI’s growth, and the interaction between charge-transfer resistance and double-layer capacitance accelerated. Ageing tests also confirmed that the anode was the faster electrode for the considered battery cell [51].

Accordingly, it is plausible that the negative electrode also formed a faster time constant ${\tau }_{\mathrm{dl}}^{-}$ than the positive electrode. Another indication is the lower SOC dependence of ${\tau }_{\mathrm{dl}}^{-}$ (see Figure 5e,f) compared to ${\tau }_{\mathrm{dl}}^{+}$ (see Figure 5h,i) since the surface concentration was independent of the battery’s SOC over long ranges and remained constant due to the two-phase material behaviour of graphite. Only a slight increase in the time constant can be seen towards the end of the discharge when a pure phase was formed. The time constant ${\tau }_{\mathrm{dl}}^{-}$ shows this dependence clearly, indicating that a single reaction process could be identified here.

The time constant ${\tau }_{\mathrm{dl}}^{-}$ was determined using DRT process fitting [63] (see Figure 4d). The associated resistance $R_{\mathrm{ct}}^{-}$ was subsequently identified by parameter variation in the impedance spectrum. Although steady curves over the SOC were identified here, the dependence on temperature did not exhibit a clear trend (see Figure A2). Comparing $R_{\mathrm{ct}}^{-}$ (see Figure 5g) with $R_{\mathrm{ser}}$, it is noticeable that the ohmic resistance showed a piecewise linear behaviour with one region from −15 °C to 15 °C and one region from 15 °C to 35 °C. The deviation was of approximately the same order of magnitude as $R_{\mathrm{ct}}^{-}$ in these temperature ranges, while $R_{\mathrm{ct}}^{-}$ with exponential extrapolation was below the measurement accuracy. Due to the high dynamics of the reaction process, the impedance fitting algorithm did not seem to be able to distinguish clearly between the pure resistive resistance $R_{\mathrm{ser}}$, the migration resistance $R_{\mathrm{M},\mathrm{P}}$ in the pores (see Figure 6a,d), and the charge-transfer resistance $R_{\mathrm{ct}}^{-}$. Therefore, measuring both model parts in the low-temperature range is sufficient, extrapolating values above 15 °C, because $R_{\mathrm{ct}}^{-}$ would then fall to minimal values.

3.2.3. Reaction of the Positive Electrode

The reaction process of the positive electrode was in a lower frequency range and was much easier to isolate and characterise. This reaction process also slowed down exponentially with decreasing temperature. The approximately linear shape of the curve of ${\tau }_{\mathrm{dl}}^{+}$ in the Arrhenius plot in Figure 5i indicates the splendid isolation of this process by parameter identification.

The time constant ${\tau }_{\mathrm{dl}}^{+}$ of the reaction process increased significantly as the battery’s SOC decreased. The double-layer capacitance’s relaxation and recharging became very slow at low temperatures. Here, the DRT method was again used for identification. Since the time constant of diffusion inside the particles showed much less dependence on temperature and even tended to accelerate at low SOCs, the separation of these states would not be valid without further consideration. For a more straightforward calculation of the voltage in the time domain, this separation is necessary to split up the parameters of the Randles circuit.

The charge-transfer resistance of the positive electrode $R_{\mathrm{ct}}^{+}$, presented in Figure 5b,c, shows similar behaviour to the associated time constant. This behaviour suggests that the capacitance of the double layer remained comparatively stable over SOCs and temperatures.

Additionally, $R_{\mathrm{ct}}^{+}$ increased significantly at lower SOCs. There are several explanations for this. For example, the exchange current density—and thus directly the charge-transfer resistance—was influenced by the concentrations of the individual reaction reductants, as described by Newman’s model [64]. Towards the end of the battery’s discharge, the lithium concentration in the positive electrode particles slowly approached the maximum. The free lattice spaces necessary for the intercalation of lithium became scarce, and the reaction required higher overpotentials to continue at the same rate.

The second explanation for such an increase supports the hypothesis that $R_{\mathrm{ct}}^{+}$ represents the positive electrode reaction process. According to Whittingham [65] and Huggins [66], at very high degrees of lithiation, nickel-containing mixed oxide cathode materials, such as NMC, exhibit low electrical conductivity. Since the reaction also requires electrons, this poorer conductivity results in an increase in the reaction resistance, among other effects.

3.2.4. Intra-Particle Diffusion

For fitting the low-frequency region, we neglected the graphite electrode’s influence due to the flat course of the OCV. Since the measurement accuracy was very low in the low-frequency range and only a few measurement points were recorded, separating the two electrodes in the impedance spectrum was impossible with the available means. Accordingly, a pure cathode model was assumed for the small-signal behaviour. With this limitation, we tended to overestimate the parameter $R_{\mathrm{D}}^{+}$. Due to the limited frequency range, the parameter value for ${\tau }_{\mathrm{D}}^{+}$ can probably only be interpreted qualitatively. Considering these limitations in the interpretation, conclusions can be drawn from the parameters obtained in the impedance spectrum. Parameters characterising the large-signal behaviour were identified from the time-domain measurements to implement diffusion in the interior of the active mass particles into a time-domain model. These included ${\tau }_{\mathrm{D},\mathrm{P}}$, $R_{\mathrm{D},\mathrm{P}}$, ${\tau }_{\mathrm{L},\mathrm{F}}$, $R_{\mathrm{L},\mathrm{F}}$, as described in Section 3.2.6 and Section 3.2.7.

Although we measured the small-signal impedance at all mentioned temperatures and performed parameter fitting, the diffusion time constant below 15 °C was outside the measured frequency range, and only its initial part could be analysed. A stable extraction of the necessary parameters was not possible. For this reason, we only used the three upper measured temperatures for the parametrisation of the diffusion.

The time constant ${\tau }_{\mathrm{D}}^{+}$ can be seen in Figure 6c,f. This parameter is divided into two ranges: above 70% SOC, a much slower relaxation occurred compared to below this threshold. Also, the voltage swing, proportional to the diffusion resistance $R_{\mathrm{D}}^{+}$, as shown in Figure 6b,e, was much higher above this limit. It is likely that this limit coincidentally fell near the same SOC at which the change from the 80 mV to the 120 mV plateau occurred on the negative electrode. The very flat characteristic of graphite in the highly lithiated state rules out a significant effect on $R_{\mathrm{D}}^{+}$ in this simple model. Accordingly, the effect must have been located on the positive electrode alone.

For the same intercalation material, which did not exhibit significantly different intercalation mechanisms and did not form phases of different concentrations, it is implausible that such a sharp transition would occur. Due to the high amplitude of the overpotential, overlapping effects from the negative electrode can be ruled out.

3.2.5. Porous Electrode

In this work, it is assumed that the reaction overpotential of an electrode can never be considered isolated. The impedance behaviour is always an interaction between the reaction overpotential over $R_{\mathrm{ct}}$, the double-layer capacitance $C_{\mathrm{dl}}$, and the migration resistance $R_{\mathrm{M},\mathrm{P}}$ of the porous electrode structure. It is an arc in the Nyquist diagram that exhibits a 45° angle at high frequencies and the behaviour of an RC circuit ($R_{\mathrm{ct}}\left|\right|C_{\mathrm{dl}}$) at low frequencies.

The positive and negative electrodes have different migration resistances $R_{\mathrm{M},\mathrm{P}}$. The negative electrode’s high-frequency part was not entirely recorded over a wide temperature range. Therefore, a separate fitting of $R_{\mathrm{M},\mathrm{P}}^{-}$ with the available data was not possible. However, since the geometric structure of the two electrodes should not differ significantly, and the same electrolyte is present in both electrodes, it is assumed that both migration resistances are identical:

$$R_{\mathrm{M},\mathrm{P}}=R_{\mathrm{M},\mathrm{P}}^{-}=R_{\mathrm{M},\mathrm{P}}^{+}$$

Figure 6a,d show the variation in this parameter over the battery’s temperature and SOC. The linear Arrhenius plot again indicates an excellent identification of a single physical process. Notable for this parameter is the substantial, abrupt increase in resistance towards low SOCs, especially at low temperatures. The model concept underlying the equations cannot entirely explain this increase. The battery’s SOC should not significantly influence the migration resistance in the pores of the electrodes. Particularly on the positive electrode, a change in the volume of the particles is hard to expect, and on the negative electrode, the particles should be smallest in this state, making the pore volume as a whole the largest. With respect to usability, we extended the model complexity here, but the SOC dependencies of the parameters compensated for the shortcomings. However, when interpreting the results, this has to be considered.

3.2.6. Diffusion

The diffusive elements of the electrolyte concentration overpotential can be seen in Figure 6g,h for $R_{\mathrm{D},\mathrm{P}}$ and Figure 7a,d for ${\tau }_{\mathrm{D},\mathrm{P}}$. Except for the measurements at temperatures below 0 °C, the parameter identification of this RC element yielded physically reasonable results. At this point, it should be mentioned that in some publications [67,68], this separable RC element with time constants in the two-digit or small three-digit second range is attributed to the equalisation of two different particle groups or electrode layers. However, this cannot be separated from electrolyte diffusion in a porous electrode. As described, an interaction between the equalisation of the individual electrode layers and the electrolyte salt diffusion must always be considered. Meanwhile, the existence of the separable RC element can only be postulated with the inclusion of both electrodes.

The resulting time constant ${\tau }_{\mathrm{D},\mathrm{P}}$ and the effective resistance $R_{\mathrm{D},\mathrm{P}}$ exhibited behaviour that was qualitatively aligned with the gradient of the open-circuit voltage characteristic of the battery. In the small-signal behaviour, this was the inverse of the limiting capacitance of the particle diffusion impedance ${\underline{Z}}_{\mathrm{Diff}}$. In regions with low gradients and correspondingly high limiting capacitances, the compensation of the electrolyte concentration differences between individual layers occurred more slowly than in regions with significant gradients. At the same time, both parameters ${\tau }_{\mathrm{D},\mathrm{P}}$ and $R_{\mathrm{D},\mathrm{P}}$ showed a slight dependence on temperature. This dependency indicates the influence of diffusion on the electrolyte since its diffusion constant is affected by temperature.

3.2.7. Time-Domain Parameters

In the large-signal relaxation behaviour (see Figure 4e), the corresponding equivalent circuit elements were also parameterised for the slow homogenisation process with different particles. Therefore, we transformed the elements with already identified parameters into the time domain and calculated the voltage response of the model using all the elements identified in the frequency domain. We used the elements that had not been parameterised so far to minimise the voltage difference from the voltage course of the relaxation measurement (see Figure 4a). The time constant of this prolonged exchange process ${\tau }_{\mathrm{LF}}$ can be seen in Figure 7c,f. ${\tau }_{\mathrm{LF}}$ behaved qualitatively similarly to the time constant of homogenisation inside the particles ${\tau }_{\mathrm{D}}^{+}$, where ${\tau }_{\mathrm{LF}}$ was always significantly slower than ${\tau }_{\mathrm{D}}^{+}$. This correlation substantiates the hypothesis that this represents the exchange between particles. The homogenisation across different particle groups was closely related to particle diffusion and was always slower.

Figure 7e clearly shows an Arrhenius dependence of the resistance $R_{\mathrm{LF}}$, indicating that a single process was again identified. However, it is unclear whether the balancing mechanism was due to the diffusion of the intercalated material across grain boundaries or via recharge involving ionic and electrical exchange. In the region above 70% SOC, extremely slow voltage relaxations occurred. Extremely large resistances $R_{\mathrm{LF}}$ and time constants ${\tau }_{\mathrm{LF}}$ were required to describe the slow voltage relaxation at certain high SOCs (see Figure 7b). Implementing these values in a time-domain model would lead to extreme overpotentials at precisely these points. However, the battery is a parallel circuit of countless particles. Therefore, when the battery is discharged, the entire electrode never has the same degree of lithiation. In the dynamic voltage response, such a particular SOC with extremely high overpotential and prolonged relaxation does not occur unless explicitly set in the parameterisation. It is, therefore, entirely legitimate to discard extreme parameter values occurring at a single, selective SOC to obtain a continuous voltage while exceeding this SOC. Alternatively, one would have to simulate the parallel connection of different particles, which would require considerable additional computational effort.

3.2.8. Voltage Source

As the battery has an SOC-dependent OCV, it is important to parameterise it as an SOC-dependent voltage source $V_{\mathrm{OCV}}$, as shown in Figure 6i. The overpotentials across all equivalent circuit elements caused by an electrical load add to the OCV.

3.3. Model Verification

In the last subsection of this section, we show that the simulated voltage response to a specified current profile follows the expected course. The verification of the behaviour for different operating points is discussed below.

3.3.1. Simulation

We performed the simulation using the ISEAFrame tool [69]. This tool can be used to simulate variable ECMs. A current or power profile can serve as its input. Here, we used a current profile as input. The model input for the simulation tool consisted of SOC/T look-up tables for each parameter, with the measured SOCs and temperatures as supporting points. For operation between two values in the look-up tables, linear interpolation was applied.

The current profiles used for verification were measured electrically on the test bench, and their voltage curves can be considered for comparison. The verification was performed on current profiles measured at mostly constant SOCs. Thus, the simulated behaviour of the parameterised model can be verified locally. A simulation with a dynamic driving profile for validation was also conducted; however, the analysis is beyond the scope of this work.

3.3.2. Comparison

First, the voltage curve for a pulse profile at 25 °C and 80% SOC is considered (see Figure 8a). It is immediately apparent that the voltage curve is very close to the measured voltage curve, both in peak values and in dynamics. The model reacted as expected to different current strengths, whether in the charging or discharging direction.

Figure 8b clearly shows that the error peaks occurred almost exclusively at the moments of current changes. This behaviour is examined in more detail in Figure 8e. When only the moments of current changes are considered, an instantaneous jump in the voltage response due to the ohmic components of the internal resistance can be seen in the voltage curve of the simulation, as expected. In the measurement, a voltage increase in the form of a ramp is more apparent. This ramp resulted from the control of the cyclers, whose rising steepness was not infinitely large; therefore, a current ramp, instead of a current jump, was applied. Nevertheless, we take these errors into account.

The voltage error between the simulation and measurement is shown in Figure 8c. Here, an average error of less than 10 mV can be seen. Only during load changes do peak values of up to 100 mV appear. Looking at the different simulated current profiles under varying operating conditions, we observe that the simulation yielded excellent results across the entire temperature and SOC ranges. The figure shows current profiles simulated for three different temperatures (−5 °C, 5 °C, and 25 °C) and three different SOCs (20%, 50%, and 80%), with the voltage error between the measurement and simulation plotted over time. Maximum peak deviations of 300 mV can be explained by the offset in the current rise between the measurement and the simulation. When considering the mean difference, the simulation has an error of less than 10 mV over the wide operating range (see Figure 8f). Only for negative temperatures do we see higher deviations, which might be due to effects at low temperatures that have not been modelled or parameterised. Another possible explanation is that while the validation profile is experimentally set at an ambient temperature of −5 °C, the temperature of the battery cell might be different.

Additionally, Figure 8d shows a comparison of the voltage curves of the simulation and measurement regarding a 1C CC discharge. We can see quite an accurate representation of the real-world measured behaviour, with only slight inaccuracies at very small SOCs around 3400 s.

This overview shows that a physics-based parameter set does not only achieve optimal accuracy for one measuring point. It also highlights the trade-off between local accuracy and global physical reasonability in parameterisation. While local results might be improved, the ability to extrapolate could be lost, and the deviation from physical reasonability could become more significant.

4. Conclusions

We comprehensively investigated what information about the battery cell can be extracted from the current and voltage behaviour alone. We started with the basic physical equations of diffusion and migration and the equilibrium conditions of electrochemistry to develop a model that describes the influence of each relevant process on the battery voltage.

Based on the concept of complex impedance, the respective transfer functions from a current excitation to a voltage response were derived for all relevant phenomena. At this point, it became clear that there are natural limits to extracting individual physical quantities from the voltage response alone. For example, whenever two quantities occur exclusively in superposition, one can only be inferred if the other is known from additional sources. For further interpretation, the parameter values were combined so that each impedance element has a minimum number of degrees of freedom.

The presented model employs some simplifications that limit its usability. The model does not take hysteresis effects into account. The resulting offset of the equilibrium voltage can be addressed by replacing the SOC-OCV curve with a hysteresis model. The effects of hysteresis on the cell’s impedance should be addressed in future work. The parameterisation technique neglects the current dependence of the transfer resistance $R_{\mathrm{ct}}$, thereby limiting its applicability to low C rates. This shortcoming should be addressed in future work to enable the model for HP cells. Setting $\phi$ = 1 (Equation (49)) reduces the model’s runtime by approximating ZARC elements as RC circuits but reduces accuracy and physical interpretability concerning particle distributions in the electrodes.

Correspondingly, electrical parameter quantities for the impedances are derived from resistances, time constants, and quantities without units. In a porous electrode, the individual reaction and transport processes are interconnected in chain ladder-like circuits or linked via source terms. The derived impedances were transferred into a mathematical framework based on these circuits, which allows for describing the voltage behaviour of the entire cell.

In the first application, the possibilities of an ECM derived from electrical measurements were demonstrated using a complete parameterisation of the first-generation Mitsubishi iMiEV’s battery cell. In contrast to simpler ECMs, this model is valid over a wide operating range, allowing for the separation of the overpotentials of both electrodes and for extrapolation of the parameters and, therefore, the battery cell’s electrical behaviour.

The verification results show that even for a load profile with varying dynamics, the time constants of the battery’s processes are sufficiently met, and the voltage error is low across the entire operational range. For a single local operation point, this model can achieve even more precise simulation results compared to the measurements, but with the drawback that the processes are no longer accurately represented and the physical reasonability is lost. Therefore, this ECM parameter set represents the best trade-off between accuracy and physical realism.