An Exact Closed-Form Impedance Model for Porous-Electrode Lithium-Ion Cells

The original DFN P2D model was developed to describe constant-current (not transient) response. It ignores the double-layer of charge that builds up in the electrolyte at the solid–electrolyte interface, so some critical high-frequency cell behaviors can't be obtained. We introduce a more general model for this, as shown in Fig. 1a. In this Figure, ϕ_s is the electric potential in the solid, ϕ_sf is the electric potential in the solid–electrolyte interphase (SEI) film closest to the solid, ϕ_fe is the electric potential in the film closest to the electrolyte, and ϕ_e is the electric potential in the electrolyte. It can be considered that there is a particle/film interface (between ϕ_s and ϕ_sf), a film/electrolyte interface (between ϕ_fe and ϕ_e), and an overall interface (between ϕ_s and ϕ_e). Hence, in addition to the usual faradaic reaction flux, we assume that there may also be non-faradaic (capacitive) fluxes that develop from the charging and discharging of electrical double-layers at these interfaces. Because the double-layer of charge can be modeled as capacitance,¹⁷ we can model non-faradaic flux for any of these three interfaces as

$$\begin{eqnarray}&&{{Fj}}_{\mathrm{dl}}={C}_{\mathrm{dl}}\displaystyle \frac{\partial \left({\phi }_{1}-{\phi }_{2}\right)}{\partial t},\end{eqnarray} \tag{ 1 }$$

where j_dl is the non-faradaic reaction flux, C_dl is the double-layer capacitance, ϕ₁ is the potential of the phase closest to the solid, and ϕ₂ is the potential of the phase closest to the electrolyte.

Zoom In Zoom Out Reset image size

Because the ions moving in the electrolyte also need to overcome resistance, we consider a more realistic double-layer model that includes a series resistance R_dl. In this case,

$$\begin{eqnarray}&&{{Fj}}_{\mathrm{dl}}+{{FR}}_{\mathrm{dl}}{C}_{\mathrm{dl}}\displaystyle \frac{\partial {j}_{\mathrm{dl}}}{\partial t}={C}_{\mathrm{dl}}\displaystyle \frac{\partial ({\phi }_{1}-{\phi }_{2})}{\partial t}.\end{eqnarray} \tag{ 2 }$$

For cells where R_dl is considered to be negligible or nonexistent, we simply set R_dl = 0 in Eq. 2 to obtain the same result as Eq. 1.

After adding the solid-electrolyte interface details as described in Section Solid-electrolyte interface, the P2D model is refined as follows:

1. The governing Equation describing charge conservation in the solid is

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial x}{\sigma }_{\mathrm{eff}}\displaystyle \frac{\partial }{\partial x}{\phi }_{s}={a}_{s}F({j}_{{\rm{f}}}+{j}_{\mathrm{dl}})\end{eqnarray} \tag{ 3 }$$

having boundary conditions

$$\begin{eqnarray}&&-\ {\sigma }_{\mathrm{eff}}{\left.\displaystyle \frac{\partial }{\partial x}{\phi }_{s}\right|}_{x=0,{L}^{n}+{L}^{s}+{L}^{p}}=\displaystyle \frac{{i}_{\mathrm{app}}}{A}\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&-\ {\sigma }_{\mathrm{eff}}{\left.\displaystyle \frac{\partial }{\partial x}{\phi }_{s}\right|}_{x={L}^{n},{L}^{n}+{L}^{s}}=0\end{eqnarray} \tag{ 5 }$$

where σ_eff is the effective electronic conductivity in the solid, ϕ_s is the electric potential in the solid, a_s is the specific interfacial surface area, F is the Faraday's constant, j_f is the reaction flux from solid to electrolyte, A is the surface area, i_app is the cell applied current and x represents the 1-d cross-Section spatial location in the cell. Lⁿ, L^s and L^p are the thickness of the negative-electrode, separator, and positive-electrode, respectively.

2. The governing Equation describing charge conservation in the electrolyte is

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial x}\left({\kappa }_{\mathrm{eff}}\displaystyle \frac{\partial }{\partial x}{\phi }_{e}\right)+\displaystyle \frac{2{RT}\left({t}_{+}^{0}-1\right)}{F}\left(1+\displaystyle \frac{\partial \mathrm{ln}{f}_{\pm }}{\partial \mathrm{ln}{c}_{e}}\right)\\ \ \ \times \ \displaystyle \frac{\partial }{\partial x}\left({\kappa }_{\mathrm{eff}}\displaystyle \frac{\partial }{\partial x}\mathrm{ln}{c}_{e}\right)=-{a}_{s}F({j}_{{\rm{f}}}+{j}_{\mathrm{dl}})\end{array}\end{eqnarray} \tag{ 6 }$$

having boundary conditions,

$$\begin{eqnarray}&&\begin{array}{l}-\ {\kappa }_{\mathrm{eff}}\left[{\left.\displaystyle \frac{\partial }{\partial x}{\phi }_{e}\right|}_{x=0,{L}^{n}+{L}^{s}+{L}^{p}}+\displaystyle \frac{2{RT}\left({t}_{+}^{0}-1\right)}{F}\right.\\ \ \ \left.{\left.\times \left(1+\displaystyle \frac{\partial \mathrm{ln}{f}_{\pm }}{\partial \mathrm{ln}{c}_{e}}\right)\displaystyle \frac{\partial \mathrm{ln}{c}_{e}}{\partial x}\right|}_{x=0,{L}^{n}+{L}^{s}+{L}^{p}}\right]=0\end{array}\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&\begin{array}{l}-\ {\kappa }_{\mathrm{eff}}\left[{\left.\displaystyle \frac{\partial }{\partial x}{\phi }_{e}\right|}_{x={L}^{n},{L}^{n}+{L}^{s}}+\displaystyle \frac{2{RT}\left({t}_{+}^{0}-1\right)}{F}\right.\\ \ \ \left.{\left.\times \left(1+\displaystyle \frac{\partial \mathrm{ln}{f}_{\pm }}{\partial \mathrm{ln}{c}_{e}}\right)\displaystyle \frac{\partial \mathrm{ln}{c}_{e}}{\partial x}\right|}_{x={L}^{n},{L}^{n}+{L}^{s}}\right]=\displaystyle \frac{{i}_{\mathrm{app}}}{A}\end{array}\end{eqnarray} \tag{ 8 }$$

where κ_eff is the effective electrolyte conductivity, ϕ_e is the electric potential in the electrolyte, R is the universal gas constant, T is the temperature, ${t}_{+}^{0}$ is the transference number of the cation in the electrolyte with respect to the solvent, f_± is the mean molar activity coefficient, and c_e is the concentration of lithium in the electrolyte.

3. The governing Equation describing mass conservation in the solid is

$$\begin{eqnarray}&&\displaystyle \frac{\partial {c}_{s}}{\partial t}=\displaystyle \frac{{D}_{s}}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}\displaystyle \frac{\partial {c}_{s}}{\partial r}\right)\end{eqnarray} \tag{ 9 }$$

having boundary conditions

$$\begin{eqnarray}&&{D}_{s}{\left.\displaystyle \frac{\partial {c}_{s}}{\partial r}\right|}_{r=0}=0\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{D}_{s}{\left.\displaystyle \frac{\partial {c}_{s}}{\partial r}\right|}_{r={R}_{s}}=-{j}_{{\rm{f}}}\end{eqnarray} \tag{ 11 }$$

where c_s is the concentration of lithium in a spherical electrode particle, r is the coordinate in the radial pseudo dimension where r = 0 at the particle center and r = R_s at the particle surface, and D_s is the solid diffusivity.

4. The governing Equation describing mass conservation in the electrolyte is

$$\begin{eqnarray}&&{\varepsilon }_{e}\displaystyle \frac{\partial {c}_{e}}{\partial t}=\displaystyle \frac{\partial }{\partial x}\left({D}_{e,\mathrm{eff}}\displaystyle \frac{\partial }{\partial x}{c}_{e}\right)+(1-{t}_{+}^{0}){a}_{s}({j}_{{\rm{f}}}+{j}_{\mathrm{dl}})\end{eqnarray} \tag{ 12 }$$

having gradient boundary conditions,

$$\begin{eqnarray}&&{D}_{e,\mathrm{eff}}{\left.\displaystyle \frac{\partial {c}_{e}}{\partial x}\right|}_{x=0,{L}^{n}+{L}^{s}+{L}^{p}}=0\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{D}_{e,\mathrm{eff}}{\left.\displaystyle \frac{\partial {c}_{e}}{\partial x}\right|}_{x={\left({L}^{n}\right)}^{-}}={D}_{e,\mathrm{eff}}{\left.\displaystyle \frac{\partial {c}_{e}}{\partial x}\right|}_{x={\left({L}^{n}\right)}^{+}}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{D}_{e,\mathrm{eff}}{\left.\displaystyle \frac{\partial {c}_{e}}{\partial x}\right|}_{x={\left({L}^{n}+{L}^{s}\right)}^{-}}={D}_{e,\mathrm{eff}}{\left.\displaystyle \frac{\partial {c}_{e}}{\partial x}\right|}_{x={\left({L}^{n}+{L}^{s}\right)}^{+}}\end{eqnarray} \tag{ 15 }$$

and continuity boundary conditions,

$$\begin{eqnarray}&&{\left.{c}_{e}\right|}_{x={\left({L}^{n}\right)}^{-}}={\left.{c}_{e}\right|}_{x={\left({L}^{n}\right)}^{+}}\end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray}&&{\left.{c}_{e}\right|}_{x={\left({L}^{n}+{L}^{s}\right)}^{-}}={\left.{c}_{e}\right|}_{x={\left({L}^{n}+{L}^{s}\right)}^{+}}\end{eqnarray} \tag{ 17 }$$

where ε_e is the porosity and D_e,eff is the effective diffusivity of the electrolyte.

5. The Butler–Volmer closure term j_f is

$$\begin{eqnarray}&&{j}_{{\rm{f}}}={j}_{0}\left\{\exp \left(\displaystyle \frac{(1-\alpha )F}{{RT}}\eta \right)-\exp \left(\displaystyle \frac{-\alpha F}{{RT}}\eta \right)\right\}\end{eqnarray} \tag{ 18 }$$

where the exchange flux density j₀ and the overpotential η are defined respectively as

$$\begin{eqnarray}&&{j}_{0}={k}_{0}{\left({c}_{e}\right)}^{1-\alpha }{\left({c}_{s,\max }-{c}_{\mathrm{ss}}\right)}^{1-\alpha }{\left({c}_{\mathrm{ss}}\right)}^{\alpha }\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\eta ={\phi }_{s}-{\phi }_{e}-{U}_{\mathrm{ocp}}({c}_{\mathrm{ss}})-{{FR}}_{\mathrm{film}}{j}_{{\rm{f}}}\end{eqnarray} \tag{ 20 }$$

where k₀ is the reaction-rate constant, c_s,max is the maximum lithium concentration in the solid, c_ss is the solid surface concentration, α is the charge-transfer coefficient, U_ocp is the open circuit potential, and R_film is the resistance of the solid surface film.

The set of partial-differential and algebraic equations in Sections Solid-electrolyte interface and Improved P2D model give the full P2D full-order model (FOM) of a lithium-ion cell. These equations will later be used to find an exact closed-form impedance model of the cell. First, however, we elaborate on how to obtain the impedance of the cell and the frequency responses of the internal electrochemical variables using a numeric time-domain simulation of the FOM.

Cell impedance is a function of frequency. It is defined for the case where cell applied current i_app is a sinusoidal excitation:

$$\begin{eqnarray}&&{i}_{\mathrm{app}}(t)=\alpha \cos (\omega t),\end{eqnarray} \tag{ 21 }$$

for $-\infty \lt t\lt \infty$, where ω is the constant angular (radian) frequency of the excitation current, α is the amplitude of the excitation current, and t is time.

As described earlier, the impedance of a cell is equal to the negative of its frequency response. However, frequency responses, by definition, exist only for linear-time-invariant systems. The PDEs describing lithium-ion cell dynamics are nonlinear, but can be very closely approximated as linear perturbations around a constant setpoint v₀ when the input current has small amplitude, which we will assume to be the case here. Under this assumption, the cell's perturbation-voltage response to the input-current stimulus will be a sinusoid having the same frequency as the input, amplified by factor $\left|G(j\omega )\right|$ and phase shifted by arg G(jω). G(jω) is the "frequency response" or "frequency function" of the system, which gives full information to predict the output under the condition $-\infty \lt t\lt \infty$.

In a practical experiment, however, the current cannot be applied starting at time $t=-\infty$. Instead, it is applied at some specific point in time that we reference as t = 0. In this case, the output will include an additional transient term that occurs due to the sudden application of the sinusoid at time zero. For example, the output voltage will be

$$\begin{eqnarray}&&v(t)={v}_{0}+\alpha \left|G(j\omega )\right|\cos (\omega t+\arg \ G(j\omega ))+\mathrm{transient},\end{eqnarray} \tag{ 22 }$$

where v(t) is the cell voltage and v₀ is the at-rest initial open-circuit voltage.

To "measure" a frequency response, we will assume that measurement noises are zero mean and integrate to zero over one period of the input sinusoid, and further assume that the transient has decayed to zero before collecting data. Then, we can obtain the cell impedance from the output voltage according to the following steps, which are based on a Fourier-series analysis method.

First, we define

$$\begin{eqnarray}&&{I}_{\cos }(\omega )={\int }_{0}^{T}(v(t)-{v}_{0})\cos (\omega t)\,{\rm{d}}t\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{I}_{\sin }(\omega )={\int }_{0}^{T}(v(t)-{v}_{0})\sin (\omega t)\,{\rm{d}}t,\end{eqnarray} \tag{ 24 }$$

where T represents the period of the excitation current.

Substituting Eq. 22 into Eqs. 23 and 24, we have

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\cos }(\omega ) & = & {\displaystyle \int }_{0}^{T}\alpha \left|G(j\omega )\right|\cos (\omega t+\arg \ G(j\omega ))\cos (\omega t)\,{\rm{d}}t\\ & = & \ \alpha \left|G(j\omega )\right|{\displaystyle \int }_{0}^{T}\left(\cos (\omega t)\cos (\arg G(j\omega ))\right.\\ & & -\ \left.\sin (\omega t)\sin (\arg G(j\omega ))\right)\cos (\omega t)\,{\rm{d}}t\\ & = & \ \alpha \left|G(j\omega )\right|\cos (\arg G(j\omega )){\displaystyle \int }_{0}^{T}{\cos }^{2}(\omega t)\,{\rm{d}}t\\ & = & \ \displaystyle \frac{\alpha T\left|G(j\omega )\right|\cos (\arg G(j\omega ))}{2}\end{array}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}\begin{array}{rcl}{I}_{\sin }(\omega ) & = & {\displaystyle \int }_{0}^{T}\alpha \left|G(j\omega )\right|\cos (\omega t+\arg \ G(j\omega ))\sin (\omega t)\,{\rm{d}}t\\ & = & \ \alpha \left|G(j\omega )\right|{\displaystyle \int }_{0}^{T}\left(\cos (\omega t)\cos (\arg G(j\omega ))\right.\\ & & -\ \left.\sin (\omega t)\sin (\arg G(j\omega ))\right)\sin (\omega t)\,{\rm{d}}t\\ & = & \ -\alpha \left|G(j\omega )\right|\sin (\arg G(j\omega )){\displaystyle \int }_{0}^{T}{\sin }^{2}(\omega t)\,{\rm{d}}t\\ & = & \ -\displaystyle \frac{\alpha T\left|G(j\omega )\right|\sin (\arg G(j\omega ))}{2}\end{array}\end{eqnarray} \tag{ 26 }$$

We can get the magnitude response $\left|\hat{G}(j\omega )\right|$ and the phase response $\arg \ \hat{G}(j\omega )$ by combining I_cos(ω) and I_sin(ω),

$$\begin{eqnarray}&&\left|\hat{G}(j\omega )\right|=\displaystyle \frac{2\sqrt{{I}_{\cos }^{2}(\omega )+{I}_{\sin }^{2}(\omega )}}{T\alpha }\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&\arg \hat{G}(j\omega )=-{\tan }^{-1}({I}_{\sin }(\omega )/{I}_{\cos }(\omega ))\end{eqnarray} \tag{ 28 }$$

Finally, we obtain the cell impedance Z(jω)

$$\begin{eqnarray}&&Z(j\omega )=-\left|\hat{G}(j\omega )\right|{e}^{j\arg \hat{G}(j\omega )}.\end{eqnarray} \tag{ 29 }$$

In Section Result Comparison, we will use this method to calculate cell impedance numerically. Similarly, the proposed method is also applicable to find the frequency responses of the electrochemical variables inside cells, such as the electric potential in the electrolyte ϕ_e and in the solid ϕ_s, the solid-electrolyte potential difference ${\phi }_{{\rm{s}}{\rm{-}}{\rm{e}}}$, the concentration of lithium in the electrolyte c_e and at the solid surface c_ss, and the reaction flux j. That is, the frequency-response information of internal electrochemical variables also can be obtained simply by replacing the cell voltage v(t) in these equations with these variables.

Impedance at particle/film interface

The primary distinction between the process followed in this Paper to create frequency responses and the process presented by Rodríguez et al. is that we have modified the procedure to allow for a more general model of the solid–electrolyte interface, including the possibility of one or more electrical double layers at the interface. The flowchart presented in Fig. 2 illustrates the revised process, where the most significant change is to the steps required to find this interfacial impedance model and to the resulting expression for this impedance. The position of these new steps in the overall method for finding all transfer functions is illustrated by the top-right box labeled "Interfacial impedance model". Following the determination of this new impedance model, similar steps to those used by Rodríguez et al. are used to derive a transfer function for perturbations to the electrolyte concentration caused by applied electrical current, where the only modifications to the derivation are to allow use of this new general impedance model vs the specific impedance model previously assumed. Once the electrolyte-concentration transfer function is obtained, there are no additional changes in the procedure followed to find the remaining cell transfer functions, which are all dependent on the electrolyte-concentration transfer function. In the Figure, the blue boxes represent time-domain Equations, the yellow boxes represent frequency-domain equations leading up to the electrolyte-concentration transfer function, and the green boxes represent frequency-domain equations resulting from the electrolyte-concentration transfer function.

Zoom In Zoom Out Reset image size

We begin to see how to include the electrical double-layer model by investigating impedance at the particle/film interface in the frequency domain. According to the solid-concentration PDE stated in Eq. 9 and the conversion approach from the time domain to the frequency domain derived by Jacobsen and West,³⁰ we can obtain the transfer function ${\tilde{C}}_{\mathrm{ss}}(s)/{J}_{{\rm{f}}}^{\mathrm{sf}}(s)$ as

$$\begin{eqnarray}&&\displaystyle \frac{{\tilde{C}}_{\mathrm{ss}}(s)}{{J}_{{\rm{f}}}^{\mathrm{sf}}(s)}=\displaystyle \frac{{R}_{{\rm{s}}}}{{D}_{{\rm{s}}}}\left[\displaystyle \frac{1}{1-{R}_{{\rm{s}}}\sqrt{s/{D}_{{\rm{s}}}}\coth \left({R}_{{\rm{s}}}\sqrt{s/{D}_{{\rm{s}}}}\right)}\right].\end{eqnarray} \tag{ 30 }$$

where ${\tilde{c}}_{\mathrm{ss}}(t)={c}_{\mathrm{ss}}(t)-{c}_{{\rm{s}},0}$, c_ss(t) is the concentration of lithium at the surface of the solid and c_s,0 is the initial equilibrium concentration of lithium in the solid. ${\tilde{C}}_{\mathrm{ss}}(s)$ and ${J}_{{\rm{f}}}^{\mathrm{sf}}(s)$ are the representations of ${\tilde{c}}_{\mathrm{ss}}(t)$ and ${j}_{{\rm{f}}}^{\mathrm{sf}}(t)$ in the Laplace frequency domain, respectively, "s" is the Laplace generalized frequency variable, and ${j}_{{\rm{f}}}^{\mathrm{sf}}$ is the faradaic reaction flux from solid to film, ${j}_{{\rm{f}}}^{\mathrm{sf}}={j}_{{\rm{f}}}=j-{j}_{\mathrm{dl}}$.

We can solve for ${J}_{{\rm{f}}}^{\mathrm{sf}}(s)$ by linearizing the Butler–Volmer Equation for this interface and taking the Laplace transform of the result. First, we form a Taylor-series expansion around the linearization setpoint

$$\begin{eqnarray*}&&{p}^{* }=\{{c}_{\mathrm{ss}}={c}_{{\rm{s}},0},{c}_{{\rm{e}}}={c}_{{\rm{e}},0},{\phi }_{{\rm{s}}}-{\phi }_{\mathrm{sf}}={U}_{\mathrm{ocp}}({c}_{{\rm{s}},0}),{j}_{{\rm{f}}}^{\mathrm{sf}}=0\},\end{eqnarray*}$$

where c_{e, 0} is the equilibrium electrolyte concentration, and discard higher-order terms. This gives

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \frac{{j}_{{\rm{f}}}^{\mathrm{sf}}}{{j}_{0}} & \approx & \displaystyle \frac{F}{{RT}}\mathop{\underbrace{\left({\phi }_{{\rm{s}}}-{\phi }_{\mathrm{sf}}-{U}_{\mathrm{ocp}}({c}_{{\rm{s}},0})\right)}}\limits_{{\tilde{\phi }}_{{\rm{s}}}-{\phi }_{\mathrm{sf}}}-\displaystyle \frac{F}{{RT}}{\left[\displaystyle \frac{\partial {U}_{\mathrm{ocp}}({c}_{{\rm{s}}})}{\partial {c}_{{\rm{s}}}}\right|}_{{c}_{{\rm{s}},0}}\\ & & \times \ \mathop{\underbrace{\left({c}_{\mathrm{ss}}-{c}_{{\rm{s}},0}\right)}}\limits_{{\tilde{c}}_{\mathrm{ss}}},\end{array}\end{eqnarray*}$$

where ${\tilde{\phi }}_{{\rm{s}}}={\phi }_{{\rm{s}}}-{U}_{\mathrm{ocp}}({c}_{{\rm{s}},0})$, and ${j}_{0}={k}_{0}{({c}_{{\rm{e}},0})}^{1-\alpha }$ ${({c}_{{\rm{s}},\max }-{c}_{{\rm{s}},0})}^{1-\alpha }{({c}_{{\rm{s}},0})}^{\alpha }$.

Rearranging, we get

$$\begin{eqnarray*}&&{\tilde{\phi }}_{{\rm{s}}}-{\phi }_{\mathrm{sf}}={{FR}}_{\mathrm{ct}}^{\mathrm{sf}}{j}_{{\rm{f}}}^{\mathrm{sf}}+{\left[\displaystyle \frac{\partial {U}_{\mathrm{ocp}}({c}_{{\rm{s}}})}{\partial {c}_{{\rm{s}}}}\right|}_{{c}_{{\rm{s}},0}}{\tilde{c}}_{\mathrm{ss}},\end{eqnarray*}$$

where ${R}_{\mathrm{ct}}^{\mathrm{sf}}={RT}/({j}_{0}{F}^{2})$. Taking the Laplace transform,

$$\begin{eqnarray}&&{\tilde{{\rm{\Phi }}}}_{{\rm{s}}}(s)-{{\rm{\Phi }}}_{\mathrm{sf}}(s)={{FR}}_{\mathrm{ct}}^{\mathrm{sf}}{J}_{{\rm{f}}}^{\mathrm{sf}}(s)+{\left[\displaystyle \frac{\partial {U}_{\mathrm{ocp}}({c}_{{\rm{s}}})}{\partial {c}_{{\rm{s}}}}\right|}_{{c}_{{\rm{s}},0}}{\tilde{C}}_{\mathrm{ss}}(s).\end{eqnarray} \tag{ 31 }$$

We define ideal particle-diffusion impedance Z_s,

$$\begin{eqnarray*}&&{Z}_{{\rm{s}}}={\left[\displaystyle \frac{\partial {U}_{\mathrm{ocp}}({c}_{{\rm{s}}})}{\partial {c}_{{\rm{s}}}}\right|}_{{c}_{{\rm{s}},0}}\displaystyle \frac{{R}_{{\rm{s}}}}{{{FD}}_{{\rm{s}}}}\left[\displaystyle \frac{1}{1-{R}_{{\rm{s}}}\sqrt{s/{D}_{{\rm{s}}}}\coth \left({R}_{{\rm{s}}}\sqrt{s/{D}_{{\rm{s}}}}\right)}\right].\end{eqnarray*}$$

Then, combining Eqs. 31 and 30, we have

$$\begin{eqnarray}&&{\tilde{{\rm{\Phi }}}}_{{\rm{s}}}(s)-{{\rm{\Phi }}}_{\mathrm{sf}}(s)=F({R}_{\mathrm{ct}}^{\mathrm{sf}}+{Z}_{{\rm{s}}}){J}_{{\rm{f}}}^{\mathrm{sf}}(s).\end{eqnarray} \tag{ 32 }$$

In addition to the faradaic flux, the interface supports a non-faradaic flux according to Eq. 2

$$\begin{eqnarray*}\begin{array}{rcl}{C}_{\mathrm{dl}}^{\mathrm{sf}}\displaystyle \frac{\partial ({\phi }_{{\rm{s}}}-{\phi }_{\mathrm{sf}})}{\partial t} & = & {C}_{\mathrm{dl}}^{\mathrm{sf}}\displaystyle \frac{\partial ({\tilde{\phi }}_{{\rm{s}}}+{U}_{\mathrm{ocp}}({c}_{{\rm{s}},0})-{\phi }_{\mathrm{sf}})}{\partial t}\\ & = & \ {{FR}}_{\mathrm{dl}}^{\mathrm{sf}}{C}_{\mathrm{dl}}^{\mathrm{sf}}\displaystyle \frac{\partial {j}_{\mathrm{dl}}^{\mathrm{sf}}}{\partial t}+{{Fj}}_{\mathrm{dl}}^{\mathrm{sf}}{C}_{\mathrm{dl}}^{\mathrm{sf}}\end{array}\end{eqnarray*}$$

Taking the Laplace transform of this relationship, we have

$$\begin{eqnarray*}&&{{sC}}_{\mathrm{dl}}^{\mathrm{sf}}\left({\tilde{{\rm{\Phi }}}}_{{\rm{s}}}(s)-{{\rm{\Phi }}}_{\mathrm{sf}}(s)\right)=F\left({{sR}}_{\mathrm{dl}}^{\mathrm{sf}}{C}_{\mathrm{dl}}^{\mathrm{sf}}+1\right){J}_{\mathrm{dl}}^{\mathrm{sf}}(s).\end{eqnarray*}$$

The interface supports a total flux

$$\begin{eqnarray*}\begin{array}{rcl}{{FJ}}^{\mathrm{sf}}(s) & = & {{FJ}}_{{\rm{f}}}^{\mathrm{sf}}(s)+{{FJ}}_{\mathrm{dl}}^{\mathrm{sf}}(s)\\ & = & \ \mathop{\underbrace{\left(\displaystyle \frac{1}{{R}_{\mathrm{ct}}^{\mathrm{sf}}+{Z}_{{\rm{s}}}}+\displaystyle \frac{{{sC}}_{\mathrm{dl}}^{\mathrm{sf}}}{{{sR}}_{\mathrm{dl}}^{\mathrm{sf}}{C}_{\mathrm{dl}}^{\mathrm{sf}}+1}\right)}}\limits_{{Y}^{\mathrm{sf}}}\left({\tilde{{\rm{\Phi }}}}_{{\rm{s}}}(s)-{{\rm{\Phi }}}_{\mathrm{sf}}(s)\right).\end{array}\end{eqnarray*}$$

where Y^sf is the admittance of the solid/film interface and Z^sf = 1/Y^sf.

Impedance at film/electrolyte interface

In this subsection, we investigate impedance at film/electrolyte interface. Our methods and steps to solve for ${J}_{{\rm{f}}}^{\mathrm{fe}}(s)$ are similar to those used to solve for ${J}_{{\rm{f}}}^{\mathrm{sf}}(s)$. At this interface, we ignore any variation of the open-circuit potential of any reaction occurring at the film/electrolyte interface due to variation in the concentration of lithium in the film or solution (this is the same assumption as made by Meyers et al.¹⁷).

Following the same method as before, we find

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{fe}}(s)-{{\rm{\Phi }}}_{{\rm{e}}}(s)={{FR}}_{\mathrm{ct}}^{\mathrm{fe}}{J}_{{\rm{f}}}^{\mathrm{fe}}(s).\end{eqnarray} \tag{ 33 }$$

where ${R}_{\mathrm{ct}}^{\mathrm{fe}}={RT}/({j}_{0}^{\mathrm{fe}}{F}^{2})$ and ${j}_{{\rm{f}}}^{\mathrm{fe}}$ is the reaction flux from film to electrolyte, ${j}_{{\rm{f}}}^{\mathrm{fe}}={j}_{{\rm{f}}}^{\mathrm{sf}}={j}_{{\rm{f}}}=j-{j}_{\mathrm{dl}}$.

Again, the interface also supports a non-faradaic flux according to Eq. 2

$$\begin{eqnarray*}&&{C}_{\mathrm{dl}}^{\mathrm{fe}}\displaystyle \frac{\partial ({\phi }_{\mathrm{fe}}-{\phi }_{{\rm{e}}})}{\partial t}={{FR}}_{\mathrm{dl}}^{\mathrm{fe}}{C}_{\mathrm{dl}}^{\mathrm{fe}}\displaystyle \frac{\partial {j}_{\mathrm{dl}}^{\mathrm{fe}}}{\partial t}+{{Fj}}_{\mathrm{dl}}^{\mathrm{fe}}\end{eqnarray*}$$

Taking the Laplace transform, we have

$$\begin{eqnarray*}&&{{sC}}_{\mathrm{dl}}^{\mathrm{fe}}\left({{\rm{\Phi }}}_{\mathrm{fe}}(s)-{{\rm{\Phi }}}_{{\rm{e}}}(s)\right)=F\left({{sR}}_{\mathrm{dl}}^{\mathrm{fe}}{C}_{\mathrm{dl}}^{\mathrm{fe}}+1\right){J}_{\mathrm{dl}}^{\mathrm{fe}}(s).\end{eqnarray*}$$

The interface supports a total flux

$$\begin{eqnarray*}\begin{array}{rcl}{{FJ}}^{\mathrm{fe}}(s) & = & {{FJ}}_{{\rm{f}}}^{\mathrm{fe}}(s)+{{FJ}}_{\mathrm{dl}}^{\mathrm{fe}}(s)\\ & = & \ \mathop{\underbrace{\left(\displaystyle \frac{1}{{R}_{\mathrm{ct}}^{\mathrm{fe}}}+\displaystyle \frac{{{sC}}_{\mathrm{dl}}^{\mathrm{fe}}}{{{sR}}_{\mathrm{dl}}^{\mathrm{fe}}{C}_{\mathrm{dl}}^{\mathrm{fe}}+1}\right)}}\limits_{{Y}^{\mathrm{fe}}}\left({{\rm{\Phi }}}_{\mathrm{fe}}(s)-{{\rm{\Phi }}}_{{\rm{e}}}(s)\right)\end{array}\end{eqnarray*}$$

where Y^fe is the admittance of the film/electrolyte interface and Z^fe = 1/Y^fe.

Overall interfacial impedance model

We now combine these results to find an expression for the impedance of the overall interface. First, the resistance of the film is considered to be purely ohmic,

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{sf}}(s)-{{\rm{\Phi }}}_{\mathrm{fe}}(s)={{FR}}_{{\rm{f}}}{J}^{\mathrm{se}}(s).\end{eqnarray} \tag{ 34 }$$

Second, we recognize that the total flux through every one of the interfaces considered so far must be equal. Therefore, J^sf = J^fe = J^se. Therefore, Eqs. 31 and 33 can be written as

$$\begin{eqnarray*}\begin{array}{rcl}{\tilde{{\rm{\Phi }}}}_{{\rm{s}}}(s)-{{\rm{\Phi }}}_{\mathrm{sf}}(s) & = & {{FZ}}^{\mathrm{sf}}{J}^{\mathrm{se}}(s)\\ {{\rm{\Phi }}}_{\mathrm{fe}}(s)-{{\rm{\Phi }}}_{{\rm{e}}}(s) & = & {{FZ}}^{\mathrm{fe}}{J}^{\mathrm{se}}(s).\end{array}\end{eqnarray*}$$

Adding Eqs. 31, 33, and 34, we get

$$\begin{eqnarray*}&&{\tilde{{\rm{\Phi }}}}_{{\rm{s}}}(s)-{{\rm{\Phi }}}_{{\rm{e}}}(s)=F\left({Z}^{\mathrm{sf}}+{R}_{{\rm{f}}}+{Z}^{\mathrm{fe}}\right){J}^{\mathrm{se}}(s).\end{eqnarray*}$$

The interfacial impedance model is nearly complete: all that remains is to consider the double layer between ϕ_s and ϕ_e, which is modeled according to Eq. 2 as

$$\begin{eqnarray*}\begin{array}{rcl}{C}_{\mathrm{dl}}^{\mathrm{se}}\displaystyle \frac{\partial ({\phi }_{{\rm{s}}}-{\phi }_{{\rm{e}}})}{\partial t} & = & {C}_{\mathrm{dl}}^{\mathrm{se}}\displaystyle \frac{\partial ({\tilde{\phi }}_{{\rm{s}}}+{U}_{\mathrm{ocp}}({c}_{{\rm{s}},0})-{\phi }_{{\rm{e}}})}{\partial t}\\ & = & \ {{FR}}_{\mathrm{dl}}^{\mathrm{se}}{C}_{\mathrm{dl}}^{\mathrm{se}}\displaystyle \frac{\partial {j}_{\mathrm{dl}}^{\mathrm{se}}}{\partial t}+{{Fj}}_{\mathrm{dl}}^{\mathrm{se}}.\end{array}\end{eqnarray*}$$

Taking the Laplace transform, we have

$$\begin{eqnarray*}&&{{sC}}_{\mathrm{dl}}^{\mathrm{se}}\left({\tilde{{\rm{\Phi }}}}_{{\rm{s}}}(s)-{{\rm{\Phi }}}_{{\rm{e}}}(s)\right)=F\left({{sR}}_{\mathrm{dl}}^{\mathrm{se}}{C}_{\mathrm{dl}}^{\mathrm{se}}+1\right){J}_{\mathrm{dl}}^{\mathrm{se}}(s).\end{eqnarray*}$$

So the total molar flow rate from solid to electrolyte is

$$\begin{eqnarray*}\begin{array}{rcl}{FJ}(s) & = & {{FJ}}^{\mathrm{se}}(s)+{{FJ}}_{\mathrm{dl}}^{\mathrm{se}}(s)\\ & = & \ \left(\displaystyle \frac{1}{{Z}^{\mathrm{sf}}+{R}_{{\rm{f}}}+{Z}^{\mathrm{fe}}}+\displaystyle \frac{{{sC}}_{\mathrm{dl}}^{\mathrm{se}}}{{{sR}}_{\mathrm{dl}}^{\mathrm{se}}{C}_{\mathrm{dl}}^{\mathrm{se}}+1}\right)\left({\tilde{{\rm{\Phi }}}}_{{\rm{s}}}(s)-{{\rm{\Phi }}}_{{\rm{e}}}(s)\right).\end{array}\end{eqnarray*}$$

Defining debiased solid-electrolyte potential difference ${\tilde{\phi }}_{{\rm{s}}{\rm{-}}{\rm{e}}}\,={\tilde{\phi }}_{{\rm{s}}}-{\phi }_{{\rm{e}}}$, we find the transfer function of the interface to be

$$\begin{eqnarray}&&\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{s}}{\rm{-}}{\rm{e}}}(s)}{J(s)}={{FZ}}_{\mathrm{se}}\end{eqnarray} \tag{ 35 }$$

where

$$\begin{eqnarray*}&&{Z}_{\mathrm{se}}={\left(\displaystyle \frac{1}{{Z}^{\mathrm{sf}}+{R}_{{\rm{f}}}+{Z}^{\mathrm{fe}}}+\displaystyle \frac{{{sC}}_{\mathrm{dl}}^{\mathrm{se}}}{{{sR}}_{\mathrm{dl}}^{\mathrm{se}}{C}_{\mathrm{dl}}^{\mathrm{se}}+1}\right)}^{-1}.\end{eqnarray*}$$

This interfacial model can be visualized using the equivalent circuit shown in Fig. 1b. The reader can verify that linear-circuits principles applied to analyze this circuit will yield the same overall impedance Z_se as found by the rigorous derivation just presented based on the PDEs themselves. Therefore, the circuit is a very helpful tool for visualizing how the individual impedances contribute to the overall impedance of the interface. We also note that this interfacial model is more general than we are likely to need. In particular, it is difficult to identify all the model parameter values corresponding to a physical cell. Instead, in the remainder of this Paper we use the simplified version of the interfacial impedance model shown in Fig. 1(c) which is achieved by setting ${R}_{\mathrm{dl}}^{\mathrm{se}}={R}_{\mathrm{ct}}^{\mathrm{fe}}={R}_{\mathrm{dl}}^{\mathrm{fe}}=0$, and ${C}_{\mathrm{dl}}^{\mathrm{fe}}={C}_{\mathrm{dl}}^{\mathrm{se}}=\infty$. For this model, Z_se = Z_sf + R_f.

Now that the interfacial impedance Z_se has been derived, the next step is to use its value in the development of a transfer function for the electrolyte concentration. The process that we follow here is identical to that of Rodríguez et al.²⁶ except that we must substitute this different interfacial impedance at the correct point in the derivation.

We begin by subtracting Eq. 6 from Eq. 3 to obtain the PDE of solid-electrolyte potential difference

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{\phi }_{{\rm{s}}{\rm{-}}{\rm{e}}}(z,t)}{\partial {z}^{2}}={a}_{{\rm{s}}}{{FL}}^{2}\left(\displaystyle \frac{1}{{\sigma }_{\mathrm{eff}}}+\displaystyle \frac{1}{{\kappa }_{\mathrm{eff}}}\right)j(z,t)+\displaystyle \frac{{\kappa }_{D,\mathrm{eff}}}{{\kappa }_{\mathrm{eff}}}\displaystyle \frac{{\partial }^{2}\mathrm{ln}{c}_{{\rm{e}}}(z,t)}{\partial {z}^{2}}.\end{eqnarray} \tag{ 36 }$$

Until now, we have used natural spatial variable x; we now begin to use normalized spatial variable z to describe the one-dimensional span of the electrode (location z = 0 represents the location at the current collector; location z = 1 represents the location at the electrode/separator boundary).

We define the debiased electrolyte concentration as ${\tilde{c}}_{{\rm{e}}}(z,t)\,={c}_{{\rm{e}}}(z,t)-{c}_{{\rm{e}},0}$. Linearizing using a Taylor-series expansion, we have,

$$\begin{eqnarray*}&&\mathrm{ln}({c}_{{\rm{e}}}(z,t))=\mathrm{ln}({c}_{{\rm{e}},0}+{\tilde{c}}_{{\rm{e}}}(z,t))\approx \mathrm{ln}{c}_{{\rm{e}},0}+\displaystyle \frac{{\tilde{c}}_{{\rm{e}}}(z,t)}{{c}_{{\rm{e}},0}}.\end{eqnarray*}$$

So Eq. 36 becomes,

$$\begin{eqnarray*}&&\displaystyle \frac{{\partial }^{2}{\phi }_{{\rm{s}}{\rm{-}}{\rm{e}}}(z,t)}{\partial {z}^{2}}\approx {a}_{{\rm{s}}}{{FL}}^{2}\left(\displaystyle \frac{1}{{\sigma }_{\mathrm{eff}}}+\displaystyle \frac{1}{{\kappa }_{\mathrm{eff}}}\right)j(z,t)+\displaystyle \frac{{\kappa }_{D,\mathrm{eff}}}{{\kappa }_{\mathrm{eff}}{c}_{{\rm{e}},0}}\displaystyle \frac{{\partial }^{2}{\tilde{c}}_{{\rm{e}}}(z,t)}{\partial {z}^{2}}.\end{eqnarray*}$$

After taking the Laplace transform of both sides of this Equation, we have

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}{{\rm{\Phi }}}_{{\rm{s}}{\rm{-}}{\rm{e}}}(z,s)}{\partial {z}^{2}}={a}_{{\rm{s}}}{{FL}}^{2}\left(\displaystyle \frac{1}{{\sigma }_{\mathrm{eff}}}+\displaystyle \frac{1}{{\kappa }_{\mathrm{eff}}}\right)J(z,s)+\displaystyle \frac{{\kappa }_{D,\mathrm{eff}}}{{\kappa }_{\mathrm{eff}}{c}_{{\rm{e}},0}}\displaystyle \frac{{\partial }^{2}{\tilde{C}}_{{\rm{e}}}(z,s)}{\partial {z}^{2}}.\end{eqnarray} \tag{ 37 }$$

We also take the Laplace transform of Eq. 12,

$$\begin{eqnarray}&&J(z,s)=\displaystyle \frac{-{D}_{{\rm{e}},\mathrm{eff}}}{\left(1-{t}_{+}^{0}\right){a}_{{\rm{s}}}{L}^{2}}\displaystyle \frac{{\partial }^{2}{\tilde{C}}_{{\rm{e}}}(z,s)}{\partial {z}^{2}}+\displaystyle \frac{{\varepsilon }_{{\rm{e}}}}{\left(1-{t}_{+}^{0}\right){a}_{{\rm{s}}}}s{\tilde{C}}_{{\rm{e}}}(z,s).\end{eqnarray} \tag{ 38 }$$

Substitute Eq. 38 into 37, to get

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\partial }^{2}{{\rm{\Phi }}}_{{\rm{s}}{\rm{-}}{\rm{e}}}(z,s)}{\partial {z}^{2}} & = & \left(\displaystyle \frac{-{D}_{{\rm{e}},\mathrm{eff}}F\left(\tfrac{1}{{\sigma }_{\mathrm{eff}}}+\tfrac{1}{{\kappa }_{\mathrm{eff}}}\right)}{(1-{t}_{+}^{0})}+\displaystyle \frac{{\kappa }_{D,\mathrm{eff}}}{{\kappa }_{\mathrm{eff}}{c}_{{\rm{e}},0}}\right)\displaystyle \frac{{\partial }^{2}{\tilde{C}}_{{\rm{e}}}(z,s)}{\partial {z}^{2}}\\ & & +\ \displaystyle \frac{{\varepsilon }_{{\rm{e}}}{{FL}}^{2}\left(\tfrac{1}{{\sigma }_{\mathrm{eff}}}+\tfrac{1}{{\kappa }_{\mathrm{eff}}}\right)}{(1-{t}_{+}^{0})}s{\tilde{C}}_{{\rm{e}}}(z,s).\end{array}\end{eqnarray} \tag{ 39 }$$

At this point, we are ready to substitute the result of the new interfacial impedance Equation. We combine Eqs. 35 and 38, to find

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{{\rm{\Phi }}}}_{{\rm{s}}{\rm{-}}{\rm{e}}}(z,s) & = & {{FZ}}_{\mathrm{se}}\left(\displaystyle \frac{-{D}_{{\rm{e}},\mathrm{eff}}}{(1-{t}_{+}^{0}){a}_{{\rm{s}}}{L}^{2}}\displaystyle \frac{{\partial }^{2}{\tilde{C}}_{{\rm{e}}}(z,s)}{\partial {z}^{2}}\right.\\ & & +\ \left.\displaystyle \frac{{\varepsilon }_{{\rm{e}}}}{(1-{t}_{+}^{0}){a}_{{\rm{s}}}}s{\tilde{C}}_{{\rm{e}}}(z,s)\right).\end{array}\end{eqnarray} \tag{ 40 }$$

Because ${\partial }^{2}{{\rm{\Phi }}}_{{\rm{s}}{\rm{-}}{\rm{e}}}(z,t)/\partial {z}^{2}={\partial }^{2}{\tilde{{\rm{\Phi }}}}_{{\rm{s}}{\rm{-}}{\rm{e}}}(z,t)/\partial {z}^{2}$, we can combine Eqs. 39 and 40

$$\begin{eqnarray*}&&\begin{array}{l}-\displaystyle \frac{{D}_{{\rm{e}},\mathrm{eff}}}{{L}^{2}}\displaystyle \frac{{\partial }^{4}{\tilde{C}}_{{\rm{e}}}(z,s)}{\partial {z}^{4}}+{\varepsilon }_{{\rm{e}}}s\displaystyle \frac{{\partial }^{2}{\tilde{C}}_{{\rm{e}}}(z,s)}{\partial {z}^{2}}\\ \ =\ \left[\displaystyle \frac{-\left(\tfrac{1}{{\sigma }_{\mathrm{eff}}}+\tfrac{1}{{\kappa }_{\mathrm{eff}}}\right){D}_{{\rm{e}},\mathrm{eff}}{a}_{{\rm{s}}}}{{Z}_{\mathrm{se}}}+\displaystyle \frac{(1-{t}_{+}^{0}){a}_{{\rm{s}}}{\kappa }_{D,\mathrm{eff}}}{{\kappa }_{\mathrm{eff}}{c}_{{\rm{e}},0}{{FZ}}_{\mathrm{se}}}\right]\displaystyle \frac{{\partial }^{2}{\tilde{C}}_{{\rm{e}}}(z,s)}{\partial {z}^{2}}\\ \ \ \ +\ \displaystyle \frac{{a}_{{\rm{s}}}{L}^{2}\left(\tfrac{1}{{\sigma }_{\mathrm{eff}}}+\tfrac{1}{{\kappa }_{\mathrm{eff}}}\right){\varepsilon }_{{\rm{e}}}s}{{Z}_{\mathrm{se}}}{\tilde{C}}_{{\rm{e}}}(z,s).\end{array}\end{eqnarray*}$$

Rearranging, and dividing by I_app(s), we have

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{4}}{\partial {z}^{4}}\displaystyle \frac{{\tilde{C}}_{{\rm{e}}}(z,s)}{{I}_{\mathrm{app}}(s)}-{\tau }_{1}(s)\displaystyle \frac{{\partial }^{2}}{\partial {z}^{2}}\displaystyle \frac{{\tilde{C}}_{{\rm{e}}}(z,s)}{{I}_{\mathrm{app}}(s)}+{\tau }_{2}(s)\displaystyle \frac{{\tilde{C}}_{{\rm{e}}}(z,s)}{{I}_{\mathrm{app}}(s)}=0,\end{eqnarray} \tag{ 41 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{\tau }_{1}(s) & = & \left(\displaystyle \frac{1}{{\sigma }_{\mathrm{eff}}}+\displaystyle \frac{1}{{\kappa }_{\mathrm{eff}}}\right)\displaystyle \frac{{a}_{{\rm{s}}}{L}^{2}}{{Z}_{\mathrm{se}}}-\displaystyle \frac{(1-{t}_{+}^{0}){a}_{{\rm{s}}}{L}^{2}{\kappa }_{D,\mathrm{eff}}}{{D}_{{\rm{e}},\mathrm{eff}}{\kappa }_{\mathrm{eff}}{c}_{{\rm{e}},0}{{FZ}}_{\mathrm{se}}}+\displaystyle \frac{{L}^{2}{\varepsilon }_{{\rm{e}}}s}{{D}_{{\rm{e}},\mathrm{eff}}}\\ {\tau }_{2}(s) & = & \displaystyle \frac{{a}_{{\rm{s}}}{L}^{4}\left(\tfrac{1}{{\sigma }_{\mathrm{eff}}}+\tfrac{1}{{\kappa }_{\mathrm{eff}}}\right){\varepsilon }_{{\rm{e}}}s}{{D}_{{\rm{e}},\mathrm{eff}}{Z}_{\mathrm{se}}}.\end{array}\end{eqnarray*}$$

This is exactly the same form of solution as found by Rodríguez et al. except for a slightly different definition of τ₁(s) and τ₂(s) to accommodate the generalized interfacial impedance model. From this point on, the derivation is identical to that in the earlier work, using these revised definitions of τ₁(s) and τ₂(s). In particular, the generic solution is

$$\begin{eqnarray}&&\displaystyle \frac{{\tilde{C}}_{{\rm{e}}}(z,s)}{{I}_{\mathrm{app}}(s)}={c}_{1}{e}^{{{\rm{\Lambda }}}_{1}z}+{c}_{2}{e}^{-{{\rm{\Lambda }}}_{1}z}+{c}_{3}{e}^{{{\rm{\Lambda }}}_{2}z}+{c}_{4}{e}^{-{{\rm{\Lambda }}}_{2}z}\end{eqnarray} \tag{ 42 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{{\rm{\Lambda }}}_{1} & = & \sqrt{\displaystyle \frac{{\tau }_{1}(s)-\sqrt{{\tau }_{1}^{2}(s)-4{\tau }_{2}(s)}}{2}}\\ {{\rm{\Lambda }}}_{2} & = & \sqrt{\displaystyle \frac{{\tau }_{1}(s)+\sqrt{{\tau }_{1}^{2}(s)-4{\tau }_{2}(s)}}{2}}\end{array}\end{eqnarray*}$$

So, the solution in the negative electrode is

$$\begin{eqnarray}&&\displaystyle \frac{{\tilde{C}}_{{\rm{e}}}^{{\rm{n}}}(z,s)}{{I}_{\mathrm{app}}(s)}={c}_{1}^{{\rm{n}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{n}}}z}+{c}_{2}^{{\rm{n}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{n}}}z}+{c}_{3}^{{\rm{n}}}{e}^{{{\rm{\Lambda }}}_{2}^{{\rm{n}}}z}+{c}_{4}^{{\rm{n}}}{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{n}}}z}.\end{eqnarray} \tag{ 43 }$$

Because j^s(z, t) = 0 and ${\phi }_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{s}}}(z,t)$ does not exist in the separator region, the solution in the separator becomes

$$\begin{eqnarray}&&\displaystyle \frac{{\tilde{C}}_{{\rm{e}}}^{{\rm{s}}}(z,s)}{{I}_{\mathrm{app}}(s)}={c}_{1}^{{\rm{s}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{s}}}z}+{c}_{2}^{{\rm{s}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{s}}}z}.\end{eqnarray} \tag{ 44 }$$

The solution in the positive electrode is

$$\begin{eqnarray}&&\displaystyle \frac{{\tilde{C}}_{{\rm{e}}}^{{\rm{p}}}(z,s)}{{I}_{\mathrm{app}}(s)}={c}_{1}^{{\rm{p}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}z}+{c}_{2}^{{\rm{p}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}z}+{c}_{3}^{{\rm{p}}}{e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}z}+{c}_{4}^{{\rm{p}}}{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}z}.\end{eqnarray} \tag{ 45 }$$

In order to find the specific solution, the coefficients ${c}_{1}^{{\rm{n}}}$, ${c}_{2}^{{\rm{n}}}$, ${c}_{3}^{{\rm{n}}}$, ${c}_{4}^{{\rm{n}}}$, ${c}_{1}^{{\rm{s}}}$, ${c}_{2}^{{\rm{s}}}$, ${c}_{1}^{{\rm{p}}}$, ${c}_{2}^{{\rm{p}}}$, ${c}_{3}^{{\rm{p}}}$, and ${c}_{4}^{{\rm{p}}}$, must be found to satisfy the boundary conditions of the problem. The method and even the closed-form equations for doing so are identical to those in,²⁶ except for the different definitions of τ₁(s) and τ₂(s). These equations are listed in the Appendix for completeness. The transfer functions of other electrochemical variables derived based on ${\tilde{C}}_{e}(z,s)/{I}_{\mathrm{app}}(s)$ are also listed in the Appendix. Note that in these Equations, the superscript "n" denotes negative electrode, "s" denotes separator, and "p" denotes positive electrode.

After obtaining the transfer functions of all electrochemical variables, we can build a full-cell impedance model. Using the solid potential at the positive and negative current collectors (where z = 0), we have,

$$\begin{eqnarray*}&&v(t)={\phi }_{{\rm{s}}}^{{\rm{p}}}(0,t)-{\phi }_{{\rm{s}}}^{{\rm{n}}}(0,t)\end{eqnarray*}$$

Because our transfer functions compute ${\tilde{\phi }}_{{\rm{s}}}$ and not ${\phi }_{{\rm{s}}}$ directly and because the biasing term is cell voltage itself in the positive-electrode, so we cannot use the transfer functions for ${\tilde{\phi }}_{{\rm{s}}}$ to find cell voltage. Instead, we write ${\phi }_{{\rm{s}}}={\phi }_{{\rm{s}}{\rm{-}}{\rm{e}}}-{\phi }_{{\rm{e}}}$, so

$$\begin{eqnarray*}\begin{array}{rcl}v(t) & = & \left({\phi }_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{p}}}(0,t)+{\phi }_{{\rm{e}}}^{{\rm{p}}}(0,t)\right)-\left({\phi }_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{n}}}(0,t)+{\phi }_{{\rm{e}}}^{{\rm{n}}}(0,t)\right)\\ & = & \ \left({\tilde{\phi }}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{p}}}(0,t)+{U}_{\mathrm{ocp}}^{{\rm{p}}}({\theta }_{{\rm{s}},0}^{{\rm{p}}})+{\tilde{\phi }}_{{\rm{e}}}^{{\rm{p}}}(0,t)+{\phi }_{{\rm{e}}}^{{\rm{n}}}(0,t)\right)\\ & & -\ \left({\tilde{\phi }}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{n}}}(0,t)+{U}_{\mathrm{ocp}}^{{\rm{n}}}({\theta }_{{\rm{s}},0}^{{\rm{n}}})+{\tilde{\phi }}_{{\rm{e}}}^{{\rm{n}}}(0,t)+{\phi }_{{\rm{e}}}^{{\rm{n}}}(0,t)\right).\end{array}\end{eqnarray*}$$

We define debiased (linearized) voltage

$$\begin{eqnarray*}\begin{array}{rcl}\tilde{v}(t) & = & v(t)-\left({U}_{\mathrm{ocp}}^{{\rm{p}}}({\theta }_{{\rm{s}},0}^{{\rm{p}}})-{U}_{\mathrm{ocp}}^{{\rm{n}}}({\theta }_{{\rm{s}},0}^{{\rm{n}}})\right)\\ & = & \ \left({\tilde{\phi }}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{p}}}(0,t)+{\tilde{\phi }}_{{\rm{e}}}^{{\rm{p}}}(0,t)\right)-\left({\tilde{\phi }}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{n}}}(0,t)+{\tilde{\phi }}_{{\rm{e}}}^{{\rm{n}}}(0,t)\right)\\ & = & \ {\tilde{\phi }}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{p}}}(0,t)+{\tilde{\phi }}_{{\rm{e}}}^{{\rm{p}}}(0,t)-{\tilde{\phi }}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{n}}}(0,t)\end{array}\end{eqnarray*}$$

since, by definition, ${\tilde{\phi }}_{{\rm{e}}}^{{\rm{n}}}(0,t)=0$. Then, we can take the Laplace transform of $\tilde{v}(t)$ to get the transfer function

$$\begin{eqnarray*}&&\displaystyle \frac{\tilde{V}(s)}{{I}_{\mathrm{app}}(s)}=\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{p}}}(0,s)}{{I}_{\mathrm{app}}(s)}+\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{e}}}^{{\rm{p}}}(0,s)}{{I}_{\mathrm{app}}(s)}-\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{n}}}(0,s)}{{I}_{\mathrm{app}}(s)}.\end{eqnarray*}$$

This is an exact closed-form result for the frequency response of a lithium-ion battery cell, based on the assumed PDE model.

Note that an impedance understanding of a cell writes

$$\begin{eqnarray*}\begin{array}{rcl}v(t) & = & \mathrm{OCV}(z(t))-{{Zi}}_{\mathrm{app}}(t)\\ \tilde{v}(t) & = & -{{Zi}}_{\mathrm{app}}(t)\\ \displaystyle \frac{\tilde{V}(s)}{{I}_{\mathrm{app}}(s)} & = & -Z(s).\end{array}\end{eqnarray*}$$

So, we can write the cell impedance using transfer functions as

$$\begin{eqnarray}\begin{array}{rcl}Z(s) & = & -\left[\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{p}}}(0,s)}{{I}_{\mathrm{app}}(s)}+\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{e}}}^{{\rm{p}}}(0,s)}{{I}_{\mathrm{app}}(s)}-\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{n}}}(0,s)}{{I}_{\mathrm{app}}(s)}\right]\\ Z(j\omega ) & = & -\left[\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{p}}}(0,j\omega )}{{I}_{\mathrm{app}}(j\omega )}+\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{e}}}^{{\rm{p}}}(0,j\omega )}{{I}_{\mathrm{app}}(j\omega )}-\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{n}}}(0,j\omega )}{{I}_{\mathrm{app}}(j\omega )}\right]\end{array}\end{eqnarray} \tag{ 46 }$$

since $Z(j\omega )={\left.Z(s)\right|}_{s=j\omega }$.

We begin by examining overall cell impedance response. For the full-order model, we solve the PDEs in Section Improved P2D model using COMSOL commercial software in the time domain, simulating the cell's voltage response to 200 periods of an input-current sinusoid to allow time for the output transient to decay to a negligible level. We refined the FEM mesh until results were mesh independent: a very fine mesh was required for the 2-d domain on which the solid concentration was evaluated (600 elements per electrode), which resulted in slow simulations, with median simulation time of over 20 minutes per frequency on a computer having 16 Gb RAM and a multicore 2.8 GHz Intel i7 processor. The model parameter values were chosen from the widely reported Doyle–Fuller–Newman (DFN) set,¹⁵ and are listed in Table I. The cell SOC is set to 80%. To show that the proposed approach works both with and without an electrical double-layer model, an electrical double layer was included in the negative electrode (by adding the appropriate parameters to the standard DFN set) but not in the positive electrode (by setting those parameter values to zero). The angular frequency range of the cell applied current i_app was set to 10⁻⁴ ∼ 10³ rads⁻¹ (that is 1.59 × 10⁻⁵ ∼ 159 Hz), with six frequency points per decade. The amplitude value of the cell applied current i_app was set to 0.001 A. After the COMSOL simulations, one period of output voltages corresponding to several different exemplary input-current angular frequencies are shown in the magenta curves in Fig. 3a. Then we used the impedance-identification method described in Section Impedance calculation using full-order model to obtain the Nyquist plot of cell impedance. The results are shown as a dotted curve in Fig. 3b. To verify the correctness of the proposed impedance-identification method, the results of $\left|\hat{G}({e}^{j\omega })\right|$ and $\arg \hat{G}({e}^{j\omega })$ are returned to Eq. 22. The resulting voltages are shown in the blue curves in Fig. 3a, and are consistent with the simulation voltages, so correctness of the methodology is verified.

Zoom In Zoom Out Reset image size

We also computed cell impedance using the closed-form frequency responses using the method presented in Section Transfer function of cell impedance. To do so, we used the same cell parameter values that were used for the numeric method. Because the transfer functions have low computational complexity, we computed frequency responses for 200 frequencies over the desired frequency interval. The Nyquist plot of cell impedance obtained from the transfer function are shown as a solid curve in Fig. 3b. We see that the results are consistent with the results of the full-order models. This validates the accuracy of the transfer functions and the effectiveness of the proposed method for simplifying the full-order models.

The overall computation time required by the COMSOL PDE simulations was over 10 h to compute cell impedance and all internal-variable frequency responses for a total of 31 distinct frequencies at the level of numeric accuracy needed for comparing with the transfer-function approach. In contrast, the transfer-function approach found all cell impedances for 200 distinct frequencies in a total of 5 ms using unoptimized MATLAB code on the same computer; it found all electrochemical-variable frequency responses for the same 200 frequencies in 86 ms total. This points out two benefits of the transfer-function approach: (1) it is enormously faster (on the order of 45,000,000 times faster for cell impedance at a single frequency); (2) the transfer-function approach allows computing frequency responses for only the set of variables desired, while the PDE simulation must compute all frequency responses together.

Using the closed-form transfer functions, we obtain Nyquist plots of twelve internal electrochemical variables computed for different locations within the cell electrode. Positive and negative electrode results for electrolyte concentration, solid surface concentration, reaction flux, electrolyte potential, solid potential, and solid-electrolyte potential difference are shown in Fig. 4, through Fig. 9, respectively. The dotted curves denote results obtained from COMSOL numeric simulation of the full-order models; the solid curves are those obtained directly from the transfer functions. Different colors represent different locations within the electrode (again, location z = 0 represents the location at the current collector; location z = 1 represents the location at the electrode/separator boundary). We see here that the results for the internal transfer functions exactly match the results of the full-order models. This further validates the accuracy of the closed-form transfer functions and, importantly, suggest their potential effectiveness in simplifying the full-order models for use in real-time controls. Next we present the frequency responses of each electrochemical variable and briefly describe their characteristics. Note that with the exception of electrolyte potential, the Nyquist plots for variables computed at negative and positive electrodes appear in alternate quadrants of the complex plane. This is due to the reversal of polarity with respect to the applied current for the associated electrode which brings about a 180 degree phase shift.

Frequency responses computed from transfer functions can convey important information about underlying system dynamic behavior. Applied to batteries, we normally consider electrochemical impedance spectroscopy (EIS) and confine our analysis to a spectral characterization of cell output voltage in relation to input current. But this technique examines an external measures only and is thus limited in what it can reveal internal to the cell's electrochemistry. The present work expands the scope of standard analysis by enabling the calculation of individual frequency responses for internal electrochemical variables computed at various locations within the cell.

Frequency response of electrolyte concentration

Nyquist plots for the concentration of lithium in the electrolyte are shown in Fig. 4. The different colors represent frequency responses computed at different locations across the electrode. Since the transfer function model created in Sections Solid-electrolyte interface to Improved P2D model is not a rational function, we are not able to perform meaningful analysis directly (e.g., via pole-zero determination). Nonetheless, the computed frequency responses afford us the ability to observe certain aspects of the underlying behavior.

Zoom In Zoom Out Reset image size

Electrolyte concentration Nyquist plots are clearly dominated by a first-order characteristic for all anode locations; this is representative of a diffusion process, which in the present case reflects both bulk diffusion of lithium ions through the electrolyte as well as through the porous structure of the electrode. The dominant time constant may be found as the reciprocal of the radial frequency associated with the largest magnitude imaginary component of the principle semicircle. Careful analysis of the frequency response indicates that this frequency is consistent for all locations and is found to be ${\omega }_{* }=4.095\times {10}^{-3}\,\tfrac{\mathrm{rad}}{\sec }$ which gives a corresponding diffusion time constant ${\tau }_{* }=244.2\,\sec$. Note that the actual frequency response deviates from the ideal first-order at both high and low frequency ranges indicating processes not modeled by first-order diffusion. Generally speaking, the Nyquist plot for the cathode exhibits the same essential features seen at the anode –i.e. dynamics are dominated by a first-order diffusion process. In this case, measurements show the dominant time constant varies slightly from ${\tau }_{* }=267.9\,\sec$ at the current collector location to ${\tau }_{* }=222.6\,\sec$ at the electrode surface. The key difference noted between the two electrodes is a reversal of sign which introduces a 180 degree phase shift to the plot, as we expect due to the reverse of polarity.

Frequency response of reaction flux

Nyquist plots of reaction flux are shown in Fig. 5. The individual Nyquist plot shapes vary considerably and a direct interpretation of this is unclear although may be due to variations in effective ionic conductivity as functions of both frequency and location. Nonetheless, plots for all locations are limited in their spatial extent and tend toward non-zero values at high frequencies. This essentially represents the nearly direct functional pathway between the applied current and the resultant reaction flux.

Zoom In Zoom Out Reset image size

Frequency response of potential in the electrolyte

Nyquist plots of the electric potential in the electrolyte are shown in Fig. 6. This transfer function represents a proper impedance value as it comprises the ratio of a potential (that of the electrolyte) to the applied current. In the case of both electrodes, the dynamic is dominated by two clearly separated diffusion effects, indicated by the sequential first-order semicircles in each plot. Similar to that observed with the lithium concentration, we may again compute approximate dominant time constants by measuring the radial frequency at which the maximum imaginary part magnitude is attained. For the anode the time constants are well-separated; examination of the case for location 1.0 (electrode surface) gives a dominant time constant of ${\tau }_{* }=293.87\,\sec$ with a secondary (fast) time constant of ${\tau }_{* }=0.040\,63\,\sec .$ Additional information may be gleaned from the real-axis intersection points of the Nyquist plot, which in this case can be associated with effective resistances to electron flow.

Zoom In Zoom Out Reset image size

Frequency response of potential in the solid

Nyquist plots of the electric potential in the solid are shown in Fig. 7. These responses may also be interpreted as true impedances for the corresponding electrodes and behave in similar fashion to cell impedance (to which they contribute). The basic form of the Nyquist plot for the anode again depicts first-order diffusion processes well-separated in time scale. The real axis intersection points further denote appropriate resistance values for electrode and electrolyte materials.

Zoom In Zoom Out Reset image size

The cathode exhibits similar behavior with the expected 180° phase shift. Note that both responses indicate significant variation of relative resistance values with location within the electrode.

Frequency response of the potential difference

Nyquist plots of the solid-electrolyte potential difference are shown in Fig. 8. The frequency response for the anode is very similar in form to the cell impedance obtained through EIS measurements. Regions corresponding to equilibrium capacitance (near-vertical portion), Warburg impedance (approximate 45° slope), and diffusion processes (semicircle) are clearly visible on the plot. It is of special significance that the frequency response of the potential difference (as well as other internal variables) is well-behaved in the sense that it turns out its generating components can be approximated by linear models of relatively low order. This has important implications for reduced-order modeling of these variables²⁹ for use with advanced control strategies.

Zoom In Zoom Out Reset image size

Frequency response of the solid surface concentration

Nyquist plots of the solid-electrolyte potential difference are shown in Fig. 9. The solid surface concentration essentially reflects the charge behavior of the cell and is clearly dominated by the equilibrium capacitance which appears as the vertical line in both electrode plots. Diffusion processes are also captured, but are not easily visible at the scale.

Zoom In Zoom Out Reset image size

After quite a lot of tedious back-substitution, the coefficients ${c}_{1}^{{\rm{n}}}$, ${c}_{2}^{{\rm{n}}}$, ${c}_{3}^{{\rm{n}}}$, ${c}_{4}^{{\rm{n}}}$, ${c}_{1}^{{\rm{s}}}$, ${c}_{2}^{{\rm{s}}}$, ${c}_{1}^{{\rm{p}}}$, ${c}_{2}^{{\rm{p}}}$, ${c}_{3}^{{\rm{p}}}$, and ${c}_{4}^{{\rm{p}}}$, of the electrolyte transfer function can be found exactly in closed form. The result from Rodríguez et al.²⁶ applies exactly here, with the new definition of τ₁(s) and τ₂(s) and consequential modifications to the values of (but not equations for) quantities that depend on τ₁(s) and τ₂(s). The solution is expressed in terms of

$$\begin{eqnarray*}\begin{array}{rcl}{\mathrm{den}}^{{\rm{n}}} & = & {\sigma }_{\mathrm{eff}}^{{\rm{n}}}{\kappa }_{\mathrm{eff}}^{{\rm{n}}}\left({\omega }_{-}^{{\rm{n}}}{\omega }_{-}^{{\rm{p}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}-{\omega }_{+}^{{\rm{n}}}{\omega }_{+}^{{\rm{p}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}\right)\\ {\mathrm{den}}^{{\rm{p}}} & = & {\sigma }_{\mathrm{eff}}^{{\rm{p}}}{\kappa }_{\mathrm{eff}}^{{\rm{p}}}\left({\omega }_{-}^{{\rm{n}}}{\omega }_{-}^{{\rm{p}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}-{\omega }_{+}^{{\rm{n}}}{\omega }_{+}^{{\rm{p}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}\right)\\ {{\rm{\Lambda }}}_{1} & = & \sqrt{\displaystyle \frac{{\tau }_{1}(s)-\sqrt{{\tau }_{1}^{2}(s)-4{\tau }_{2}(s)}}{2}}\\ {{\rm{\Lambda }}}_{2} & = & \sqrt{\displaystyle \frac{{\tau }_{1}(s)+\sqrt{{\tau }_{1}^{2}(s)-4{\tau }_{2}(s)}}{2}}\\ {\lambda }_{1}^{{\rm{r}}} & = & {\left({{\rm{\Lambda }}}_{1}^{{\rm{r}}}\right)}^{3}+{\mu }_{2}^{{\rm{r}}}{{\rm{\Lambda }}}_{1}^{{\rm{r}}}\\ {\lambda }_{2}^{{\rm{r}}} & = & {\left({{\rm{\Lambda }}}_{2}^{{\rm{r}}}\right)}^{3}+{\mu }_{2}^{{\rm{r}}}{{\rm{\Lambda }}}_{2}^{{\rm{r}}}\\ {{\rm{\Upsilon }}}_{1}^{{\rm{r}}} & = & {{\rm{\Lambda }}}_{1}^{{\rm{r}}}\displaystyle \frac{{D}_{e,\mathrm{eff}}^{{\rm{r}}}}{{L}_{{\rm{r}}}}\\ {{\rm{\Upsilon }}}_{2}^{{\rm{r}}} & = & {{\rm{\Lambda }}}_{2}^{{\rm{r}}}\displaystyle \frac{{D}_{e,\mathrm{eff}}^{{\rm{r}}}}{{L}_{{\rm{r}}}}\\ {\mu }_{1} & = & \displaystyle \frac{(1-{t}_{+}^{0}){a}_{s}{L}^{3}}{{{AD}}_{e,\mathrm{eff}}{{FZ}}_{\mathrm{se}}}\\ {\mu }_{2} & = & \displaystyle \frac{{\kappa }_{D,\mathrm{eff}}(1-{t}_{+}^{0}){a}_{s}{L}^{2}}{{D}_{e,\mathrm{eff}}{\kappa }_{\mathrm{eff}}{c}_{e,0}{{FZ}}_{\mathrm{se}}}-\displaystyle \frac{{\varepsilon }_{e}{L}^{2}s}{{D}_{e,\mathrm{eff}}},\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{\omega }_{+}^{{\rm{r}}} & = & \left\{\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{2}^{{\rm{r}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{r}}}}+{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{r}}}}\right)}{\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{r}}}}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{r}}}}\right)}-\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{1}^{{\rm{r}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{r}}}}+{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{r}}}}\right)}{\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{r}}}}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{r}}}}\right)}+\left({\lambda }_{2}^{{\rm{r}}}{{\rm{\Upsilon }}}_{1}^{{\rm{r}}}-{\lambda }_{1}^{{\rm{r}}}{{\rm{\Upsilon }}}_{2}^{{\rm{r}}}\right)\right\}\\ {\omega }_{-}^{{\rm{r}}} & = & \left\{\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{2}^{{\rm{r}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{r}}}}+{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{r}}}}\right)}{\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{r}}}}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{r}}}}\right)}-\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{1}^{{\rm{r}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{r}}}}+{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{r}}}}\right)}{\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{r}}}}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{r}}}}\right)}-\left({\lambda }_{2}^{{\rm{r}}}{{\rm{\Upsilon }}}_{1}^{{\rm{r}}}-{\lambda }_{1}^{{\rm{r}}}{{\rm{\Upsilon }}}_{2}^{{\rm{r}}}\right)\right\}\end{array}\end{eqnarray*}$$

for r ∈ {n, p}.

Then, the electrolyte-concentration transfer-function coefficients for the separator region are:

$$\begin{eqnarray*}\begin{array}{rcl}{c}_{1}^{{\rm{s}}} & = & -\displaystyle \frac{{\mu }_{1}^{{\rm{n}}}{\omega }_{-}^{{\rm{p}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}}{{\mathrm{den}}^{{\rm{n}}}}\left[{\sigma }_{\mathrm{eff}}^{{\rm{n}}}\left\{{{\rm{\Upsilon }}}_{2}^{{\rm{n}}}\coth ({{\rm{\Lambda }}}_{1}^{{\rm{n}}})-{{\rm{\Upsilon }}}_{1}^{{\rm{n}}}\coth ({{\rm{\Lambda }}}_{2}^{{\rm{n}}})\right\}\right.\\ & & +\ \left.{\kappa }_{\mathrm{eff}}^{{\rm{n}}}\left\{{{\rm{\Upsilon }}}_{2}^{{\rm{n}}}\mathrm{csch}({{\rm{\Lambda }}}_{1}^{{\rm{n}}})-{{\rm{\Upsilon }}}_{1}^{{\rm{n}}}\mathrm{csch}({{\rm{\Lambda }}}_{2}^{{\rm{n}}})\right\}\right]\\ & & +\ \displaystyle \frac{{\mu }_{1}^{{\rm{p}}}{\omega }_{+}^{{\rm{n}}}}{{\mathrm{den}}^{{\rm{p}}}}\left[{\sigma }_{\mathrm{eff}}^{{\rm{p}}}\left\{{{\rm{\Upsilon }}}_{2}^{{\rm{p}}}\coth ({{\rm{\Lambda }}}_{1}^{{\rm{p}}})-{{\rm{\Upsilon }}}_{1}^{{\rm{p}}}\coth ({{\rm{\Lambda }}}_{2}^{{\rm{p}}})\right\}\right.\\ & & +\ \left.{\kappa }_{\mathrm{eff}}^{{\rm{p}}}\left\{{{\rm{\Upsilon }}}_{2}^{{\rm{p}}}\mathrm{csch}({{\rm{\Lambda }}}_{1}^{{\rm{p}}})-{{\rm{\Upsilon }}}_{1}^{{\rm{p}}}\mathrm{csch}({{\rm{\Lambda }}}_{2}^{{\rm{p}}})\right\}\right]\\ {c}_{2}^{{\rm{s}}} & = & -\displaystyle \frac{{\mu }_{1}^{{\rm{n}}}{\omega }_{+}^{{\rm{p}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}}{{\mathrm{den}}^{{\rm{n}}}}\left[{\sigma }_{\mathrm{eff}}^{{\rm{n}}}\left\{{{\rm{\Upsilon }}}_{2}^{{\rm{n}}}\coth ({{\rm{\Lambda }}}_{1}^{{\rm{n}}})-{{\rm{\Upsilon }}}_{1}^{{\rm{n}}}\coth ({{\rm{\Lambda }}}_{2}^{{\rm{n}}})\right\}\right.\\ & & +\ \left.{\kappa }_{\mathrm{eff}}^{{\rm{n}}}\left\{{{\rm{\Upsilon }}}_{2}^{{\rm{n}}}\mathrm{csch}({{\rm{\Lambda }}}_{1}^{{\rm{n}}})-{{\rm{\Upsilon }}}_{1}^{{\rm{n}}}\mathrm{csch}({{\rm{\Lambda }}}_{2}^{{\rm{n}}})\right\}\right]\\ & & +\ \displaystyle \frac{{\mu }_{1}^{{\rm{p}}}{\omega }_{-}^{{\rm{n}}}}{{\mathrm{den}}^{{\rm{p}}}}\left[{\sigma }_{\mathrm{eff}}^{{\rm{p}}}\left\{{{\rm{\Upsilon }}}_{2}^{{\rm{p}}}\coth ({{\rm{\Lambda }}}_{1}^{{\rm{p}}})-{{\rm{\Upsilon }}}_{1}^{{\rm{p}}}\coth ({{\rm{\Lambda }}}_{2}^{{\rm{p}}})\right\}\right.\\ & & +\ \left.{\kappa }_{\mathrm{eff}}^{{\rm{p}}}\left\{{{\rm{\Upsilon }}}_{2}^{{\rm{p}}}\mathrm{csch}({{\rm{\Lambda }}}_{1}^{{\rm{p}}})-{{\rm{\Upsilon }}}_{1}^{{\rm{p}}}\mathrm{csch}({{\rm{\Lambda }}}_{2}^{{\rm{p}}})\right\}\right].\end{array}\end{eqnarray*}$$

The coefficients for the negative-electrode region are:

$$\begin{eqnarray*}\begin{array}{rcl}{c}_{1}^{{\rm{n}}} & = & \left({c}_{1}^{{\rm{s}}}-{c}_{2}^{{\rm{s}}}\right)\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{2}^{{\rm{n}}}}{\left({\lambda }_{2}^{{\rm{n}}}{{\rm{\Upsilon }}}_{1}^{{\rm{n}}}-{\lambda }_{1}^{{\rm{n}}}{{\rm{\Upsilon }}}_{2}^{{\rm{n}}}\right)\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}\right)}\\ & & +\ \displaystyle \frac{{\mu }_{1}^{{\rm{n}}}{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}{\left({\lambda }_{2}^{{\rm{n}}}{{\rm{\Lambda }}}_{1}^{{\rm{n}}}-{\lambda }_{1}^{{\rm{n}}}{{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}\left[\displaystyle \frac{{\sigma }_{\mathrm{eff}}^{{\rm{n}}}+{\kappa }_{\mathrm{eff}}^{{\rm{n}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}}{{\sigma }_{\mathrm{eff}}^{{\rm{n}}}{\kappa }_{\mathrm{eff}}^{{\rm{n}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}\right)}\right]\\ {c}_{2}^{{\rm{n}}} & = & \left({c}_{1}^{{\rm{s}}}-{c}_{2}^{{\rm{s}}}\right)\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{2}^{{\rm{n}}}}{\left({\lambda }_{2}^{{\rm{n}}}{{\rm{\Upsilon }}}_{1}^{{\rm{n}}}-{\lambda }_{1}^{{\rm{n}}}{{\rm{\Upsilon }}}_{2}^{{\rm{n}}}\right)\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}\right)}\\ & & +\ \displaystyle \frac{{\mu }_{1}^{{\rm{n}}}{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}{\left({\lambda }_{2}^{{\rm{n}}}{{\rm{\Lambda }}}_{1}^{{\rm{n}}}-{\lambda }_{1}^{{\rm{n}}}{{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}\left[\displaystyle \frac{{\sigma }_{\mathrm{eff}}^{{\rm{n}}}+{\kappa }_{\mathrm{eff}}^{{\rm{n}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}}{{\sigma }_{\mathrm{eff}}^{{\rm{n}}}{\kappa }_{\mathrm{eff}}^{{\rm{n}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}\right)}\right]\\ {c}_{3}^{{\rm{n}}} & = & -\left({c}_{1}^{{\rm{s}}}-{c}_{2}^{{\rm{s}}}\right)\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{1}^{{\rm{n}}}}{\left({\lambda }_{2}^{{\rm{n}}}{{\rm{\Upsilon }}}_{1}^{{\rm{n}}}-{\lambda }_{1}^{{\rm{n}}}{{\rm{\Upsilon }}}_{2}^{{\rm{n}}}\right)\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}\right)}\\ & & -\ \displaystyle \frac{{\mu }_{1}^{{\rm{n}}}{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}{\left({\lambda }_{2}^{{\rm{n}}}{{\rm{\Lambda }}}_{1}^{{\rm{n}}}-{\lambda }_{1}^{{\rm{n}}}{{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}\left[\displaystyle \frac{{\sigma }_{\mathrm{eff}}^{{\rm{n}}}+{\kappa }_{\mathrm{eff}}^{{\rm{n}}}{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}}{{\sigma }_{\mathrm{eff}}^{{\rm{n}}}{\kappa }_{\mathrm{eff}}^{{\rm{n}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}\right)}\right]\\ {c}_{4}^{{\rm{n}}} & = & -\left({c}_{1}^{{\rm{s}}}-{c}_{2}^{{\rm{s}}}\right)\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{1}^{{\rm{n}}}}{\left({\lambda }_{2}^{{\rm{n}}}{{\rm{\Upsilon }}}_{1}^{{\rm{n}}}-{\lambda }_{1}^{{\rm{n}}}{{\rm{\Upsilon }}}_{2}^{{\rm{n}}}\right)\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}\right)}\\ & & -\ \displaystyle \frac{{\mu }_{1}^{{\rm{n}}}{{\rm{\Lambda }}}_{1}^{{\rm{n}}}}{\left({\lambda }_{2}^{{\rm{n}}}{{\rm{\Lambda }}}_{1}^{{\rm{n}}}-{\lambda }_{1}^{{\rm{n}}}{{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}\left[\displaystyle \frac{{\sigma }_{\mathrm{eff}}^{{\rm{n}}}+{\kappa }_{\mathrm{eff}}^{{\rm{n}}}{e}^{{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}}{{\sigma }_{\mathrm{eff}}^{{\rm{n}}}{\kappa }_{\mathrm{eff}}^{{\rm{n}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{n}}}}\right)}\right].\end{array}\end{eqnarray*}$$

The coefficients for the positive-electrode region are:

$$\begin{eqnarray*}\begin{array}{rcl}{c}_{1}^{{\rm{p}}} & = & -\left({c}_{1}^{{\rm{s}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}-{c}_{2}^{{\rm{s}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}\right)\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{2}^{{\rm{p}}}}{\left({\lambda }_{2}^{{\rm{p}}}{{\rm{\Upsilon }}}_{1}^{{\rm{p}}}-{\lambda }_{1}^{{\rm{p}}}{{\rm{\Upsilon }}}_{2}^{{\rm{p}}}\right)\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}\right)}\\ & & -\ \displaystyle \frac{{\mu }_{1}^{{\rm{p}}}{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}{\left({\lambda }_{2}^{{\rm{p}}}{{\rm{\Lambda }}}_{1}^{{\rm{p}}}-{\lambda }_{1}^{{\rm{p}}}{{\rm{\Lambda }}}_{2}^{{\rm{p}}}\right)}\left[\displaystyle \frac{{\sigma }_{\mathrm{eff}}^{{\rm{p}}}+{\kappa }_{\mathrm{eff}}^{{\rm{p}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}}{{\sigma }_{\mathrm{eff}}^{{\rm{p}}}{\kappa }_{\mathrm{eff}}^{{\rm{p}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}\right)}\right]\\ {c}_{2}^{{\rm{p}}} & = & -\left({c}_{1}^{{\rm{s}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}-{c}_{2}^{{\rm{s}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}\right)\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{2}^{{\rm{p}}}}{\left({\lambda }_{2}^{{\rm{p}}}{{\rm{\Upsilon }}}_{1}^{{\rm{p}}}-{\lambda }_{1}^{{\rm{p}}}{{\rm{\Upsilon }}}_{2}^{{\rm{p}}}\right)\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}\right)}\\ & & -\ \displaystyle \frac{{\mu }_{1}^{{\rm{p}}}{I}_{\mathrm{app}}{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}{\left({\lambda }_{2}^{{\rm{p}}}{{\rm{\Lambda }}}_{1}^{{\rm{p}}}-{\lambda }_{1}^{{\rm{p}}}{{\rm{\Lambda }}}_{2}^{{\rm{p}}}\right)}\left[\displaystyle \frac{{\sigma }_{\mathrm{eff}}^{{\rm{p}}}+{\kappa }_{\mathrm{eff}}^{{\rm{p}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}}{{\sigma }_{\mathrm{eff}}^{{\rm{p}}}{\kappa }_{\mathrm{eff}}^{{\rm{p}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}\right)}\right]\\ {c}_{3}^{{\rm{p}}} & = & \left({c}_{1}^{{\rm{s}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}-{c}_{2}^{{\rm{s}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}\right)\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{1}^{{\rm{p}}}}{\left({\lambda }_{2}^{{\rm{p}}}{{\rm{\Upsilon }}}_{1}^{{\rm{p}}}-{\lambda }_{1}^{{\rm{p}}}{{\rm{\Upsilon }}}_{2}^{{\rm{p}}}\right)\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}\right)}\\ & & +\ \displaystyle \frac{{\mu }_{1}^{{\rm{p}}}{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}{\left({\lambda }_{2}^{{\rm{p}}}{{\rm{\Lambda }}}_{1}^{{\rm{p}}}-{\lambda }_{1}^{{\rm{p}}}{{\rm{\Lambda }}}_{2}^{{\rm{p}}}\right)}\left[\displaystyle \frac{{\sigma }_{\mathrm{eff}}^{{\rm{p}}}+{\kappa }_{\mathrm{eff}}^{{\rm{p}}}{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}}{{\sigma }_{\mathrm{eff}}^{{\rm{p}}}{\kappa }_{\mathrm{eff}}^{{\rm{p}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}\right)}\right]\\ {c}_{4}^{{\rm{p}}} & = & \left({c}_{1}^{{\rm{s}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}-{c}_{2}^{{\rm{s}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{s}}}}\right)\displaystyle \frac{{{\rm{\Upsilon }}}_{1}^{{\rm{s}}}{\lambda }_{1}^{{\rm{p}}}}{\left({\lambda }_{2}^{{\rm{p}}}{{\rm{\Upsilon }}}_{1}^{{\rm{p}}}-{\lambda }_{1}^{{\rm{p}}}{{\rm{\Upsilon }}}_{2}^{{\rm{p}}}\right)\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}\right)}\\ & & +\ \displaystyle \frac{{\mu }_{1}^{{\rm{p}}}{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}{\left({\lambda }_{2}^{{\rm{p}}}{{\rm{\Lambda }}}_{1}^{{\rm{p}}}-{\lambda }_{1}^{{\rm{p}}}{{\rm{\Lambda }}}_{2}^{{\rm{p}}}\right)}\left[\displaystyle \frac{{\sigma }_{\mathrm{eff}}^{{\rm{p}}}+{\kappa }_{\mathrm{eff}}^{{\rm{p}}}{e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}}{{\sigma }_{\mathrm{eff}}^{{\rm{p}}}{\kappa }_{\mathrm{eff}}^{{\rm{p}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}\right)}\right].\end{array}\end{eqnarray*}$$

A transfer function for the total reaction flux in region r ∈ {n, p} is

$$\begin{eqnarray}&&\displaystyle \frac{{J}^{{\rm{r}}}(z,s)}{{I}_{\mathrm{app}}(s)}={{\mathfrak{j}}}_{1}^{{\rm{r}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{r}}}z}+{{\mathfrak{j}}}_{2}^{{\rm{r}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{r}}}z}+{{\mathfrak{j}}}_{3}^{{\rm{r}}}{e}^{{{\rm{\Lambda }}}_{2}^{{\rm{r}}}z}+{{\mathfrak{j}}}_{4}^{{\rm{r}}}{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{r}}}z},\end{eqnarray} \tag{ A·47 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl}{{\mathfrak{j}}}_{1}^{{\rm{r}}} & = & {c}_{1}^{{\rm{r}}}\left(\displaystyle \frac{-{D}_{{\rm{e}},\mathrm{eff}}^{{\rm{r}}}{\left({{\rm{\Lambda }}}_{1}^{{\rm{r}}}\right)}^{2}}{(1-{t}_{+}^{0}){a}_{{\rm{s}}}^{{\rm{r}}}{L}_{{\rm{r}}}^{2}}+\displaystyle \frac{{\varepsilon }_{{\rm{e}}}^{{\rm{r}}}}{(1-{t}_{+}^{0}){a}_{{\rm{s}}}^{{\rm{r}}}}s\right)\\ {{\mathfrak{j}}}_{2}^{{\rm{r}}} & = & {c}_{2}^{{\rm{r}}}\left(\displaystyle \frac{-{D}_{{\rm{e}},\mathrm{eff}}^{{\rm{r}}}{\left({{\rm{\Lambda }}}_{1}^{{\rm{r}}}\right)}^{2}}{(1-{t}_{+}^{0}){a}_{{\rm{s}}}^{{\rm{r}}}{L}_{{\rm{r}}}^{2}}+\displaystyle \frac{{\varepsilon }_{{\rm{e}}}^{{\rm{r}}}}{(1-{t}_{+}^{0}){a}_{{\rm{s}}}^{{\rm{r}}}}s\right)\\ {{\mathfrak{j}}}_{3}^{{\rm{r}}} & = & {c}_{3}^{{\rm{r}}}\left(\displaystyle \frac{-{D}_{{\rm{e}},\mathrm{eff}}^{{\rm{r}}}{\left({{\rm{\Lambda }}}_{2}^{{\rm{r}}}\right)}^{2}}{(1-{t}_{+}^{0}){a}_{{\rm{s}}}^{{\rm{r}}}{L}_{{\rm{r}}}^{2}}+\displaystyle \frac{{\varepsilon }_{{\rm{e}}}^{{\rm{r}}}}{(1-{t}_{+}^{0}){a}_{{\rm{s}}}^{{\rm{r}}}}s\right)\\ {{\mathfrak{j}}}_{4}^{{\rm{r}}} & = & {c}_{4}^{{\rm{r}}}\left(\displaystyle \frac{-{D}_{{\rm{e}},\mathrm{eff}}^{{\rm{r}}}{\left({{\rm{\Lambda }}}_{2}^{{\rm{r}}}\right)}^{2}}{(1-{t}_{+}^{0}){a}_{{\rm{s}}}^{{\rm{r}}}{L}_{{\rm{r}}}^{2}}+\displaystyle \frac{{\varepsilon }_{{\rm{e}}}^{{\rm{r}}}}{(1-{t}_{+}^{0}){a}_{{\rm{s}}}^{{\rm{r}}}}s\right).\end{array}\end{eqnarray*}$$

Because there is no solid in the separator region, J^s(z, s) = 0.

According to Eq. 35, we can find ${\tilde{{\rm{\Phi }}}}_{{\rm{s}}{\rm{-}}{\rm{e}}}(z,s)/{I}_{\mathrm{app}}(s)$ in region r ∈ {n, p} as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{s}}{\rm{-}}{\rm{e}}}^{{\rm{r}}}(z,s)}{{I}_{\mathrm{app}}(s)} & = & {{FZ}}_{\mathrm{se}}^{{\rm{r}}}\displaystyle \frac{{J}^{{\rm{r}}}(z,s)}{{I}_{\mathrm{app}}(s)}={\mathfrak{p}}{{}_{{\mathfrak{se}}}}_{1}^{{\rm{r}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{r}}}z}+{\mathfrak{p}}{{}_{{\mathfrak{se}}}}_{2}^{{\rm{r}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{r}}}z}\\ & & +\ {\mathfrak{p}}{{}_{{\mathfrak{se}}}}_{3}^{{\rm{r}}}{e}^{{{\rm{\Lambda }}}_{2}^{{\rm{r}}}z}+{\mathfrak{p}}{{}_{{\mathfrak{se}}}}_{4}^{{\rm{r}}}{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{r}}}z},\end{array}\end{eqnarray} \tag{ A·48 }$$

where ${\mathfrak{p}}{{}_{{\mathfrak{se}}}}_{k}^{{\rm{r}}}={{\mathfrak{j}}}_{k}^{{\rm{r}}}{{FZ}}_{\mathrm{se}}^{{\rm{r}}}$.

Due to the addition of the enhanced solid/electrolyte interface description, the derivation of the transfer function for the solid surface concentration is somewhat different from in Rodríguez et al. We note that solid concentration is affected only by faradaic flux at the particle surface, not the double-layer flux (this is the same assumption made by Meyers et al.). Therefore, we can determine the solid-surface-concentration transfer function in a stepwise fashion for r ∈ {n, p}:

$$\begin{eqnarray}&&\displaystyle \frac{{\tilde{C}}_{\mathrm{ss}}^{{\rm{r}}}(z,s)}{{I}_{\mathrm{app}}(s)}=\displaystyle \frac{{\tilde{C}}_{\mathrm{ss}}^{{\rm{r}}}(z,s)}{{J}_{{\rm{f}}}^{\mathrm{sf},{\rm{r}}}(z,s)}\displaystyle \frac{{J}_{{\rm{f}}}^{\mathrm{sf},{\rm{r}}}(z,s)}{{J}^{{\rm{r}}}(z,s)}\displaystyle \frac{{J}^{{\rm{r}}}(z,s)}{{I}_{\mathrm{app}}(s)}.\end{eqnarray} \tag{ A·49 }$$

The first term is given in Eq. 30 and the third term has just been stated. However, the calculation of ${J}_{{\rm{f}}}^{\mathrm{sf}}$ depends on the exact interfacial model used.

For the default simplified model of Fig. 1c, we can write

$$\begin{eqnarray*}\begin{array}{rcl}{J}_{{\rm{f}}}^{\mathrm{sf},{\rm{r}}}(z,s) & = & \displaystyle \frac{{Z}_{\mathrm{dl}}^{\mathrm{sf},{\rm{r}}}}{{R}_{\mathrm{ct}}^{\mathrm{sf},{\rm{r}}}+{Z}_{{\rm{s}}}^{{\rm{r}}}+{Z}_{\mathrm{dl}}^{\mathrm{sf},{\rm{r}}}}{J}^{{\rm{r}}}(z,s)\\ {Z}_{\mathrm{dl}}^{\mathrm{sf},{\rm{r}}} & = & \displaystyle \frac{1+{{sR}}_{\mathrm{dl}}^{\mathrm{sf},{\rm{r}}}{C}_{\mathrm{dl}}^{\mathrm{sf},{\rm{r}}}}{{{sC}}_{\mathrm{dl}}^{\mathrm{sf},{\rm{r}}}}.\end{array}\end{eqnarray*}$$

In the negative electrode, we define ${\tilde{\phi }}_{s}(z,t)={\phi }_{s}(z,t)-{\phi }_{s}(0,t)$. Then,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{s}}}^{{\rm{n}}}(z,s)}{{I}_{\mathrm{app}}(s)} & = & -\displaystyle \frac{{{zL}}_{{\rm{n}}}}{A{\sigma }_{\mathrm{eff}}^{{\rm{n}}}}+\displaystyle \frac{{a}_{{\rm{s}}}^{{\rm{n}}}{{FL}}_{{\rm{n}}}^{2}}{{\sigma }_{\mathrm{eff}}^{{\rm{n}}}}\left(\displaystyle \frac{{{\mathfrak{j}}}_{1}^{{\rm{n}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{n}}}z}-1-z{{\rm{\Lambda }}}_{1}^{{\rm{n}}}\right)}{{\left({{\rm{\Lambda }}}_{1}^{{\rm{n}}}\right)}^{2}}\right.\\ & & +\ \displaystyle \frac{{{\mathfrak{j}}}_{2}^{{\rm{n}}}\left({e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{n}}}z}-1+z{{\rm{\Lambda }}}_{1}^{{\rm{n}}}\right)}{{\left({{\rm{\Lambda }}}_{1}^{{\rm{n}}}\right)}^{2}}+\ \displaystyle \frac{{{\mathfrak{j}}}_{3}^{{\rm{n}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{n}}}z}-1-z{{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}{{\left({{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}^{2}}\\ & & +\,\left.\displaystyle \frac{{{\mathfrak{j}}}_{4}^{{\rm{n}}}\left({e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{n}}}z}-1+z{{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}{{\left({{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}^{2}}\right).\end{array}\end{eqnarray} \tag{ A·50 }$$

In the positive electrode, we define ${\tilde{\phi }}_{s}(z,t)={\phi }_{s}(z,t)-v(t)$. This gives

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{s}}}^{{\rm{p}}}(z,s)}{{I}_{\mathrm{app}}(s)} & = & \displaystyle \frac{{{zL}}_{{\rm{p}}}}{A{\sigma }_{\mathrm{eff}}^{{\rm{p}}}}+\displaystyle \frac{{a}_{{\rm{s}}}^{{\rm{p}}}{{FL}}_{{\rm{p}}}^{2}}{{\sigma }_{\mathrm{eff}}^{{\rm{p}}}}\left(\displaystyle \frac{{{\mathfrak{j}}}_{1}^{{\rm{p}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}z}-1-z{{\rm{\Lambda }}}_{1}^{{\rm{p}}}\right)}{{\left({{\rm{\Lambda }}}_{1}^{{\rm{p}}}\right)}^{2}}\right.\\ & & +\ \displaystyle \frac{{{\mathfrak{j}}}_{2}^{{\rm{p}}}\left({e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}z}-1+z{{\rm{\Lambda }}}_{1}^{{\rm{p}}}\right)}{{\left({{\rm{\Lambda }}}_{1}^{{\rm{p}}}\right)}^{2}}+\ \displaystyle \frac{{{\mathfrak{j}}}_{3}^{{\rm{p}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}z}-1-z{{\rm{\Lambda }}}_{2}^{{\rm{p}}}\right)}{{\left({{\rm{\Lambda }}}_{2}^{{\rm{p}}}\right)}^{2}}\\ & & +\,\left.\displaystyle \frac{{{\mathfrak{j}}}_{4}^{{\rm{p}}}\left({e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}z}-1+z{{\rm{\Lambda }}}_{2}^{{\rm{p}}}\right)}{{\left({{\rm{\Lambda }}}_{2}^{{\rm{p}}}\right)}^{2}}\right).\end{array}\end{eqnarray} \tag{ A·51 }$$

A.6.1. Negative electrode

We define the unbiased electrolyte potential as ${\tilde{\phi }}_{{\rm{e}}}(z,t)={\phi }_{{\rm{e}}}(z,t)-{\phi }_{{\rm{e}}}(0,t)$. Then,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{e}}}^{{\rm{n}}}(z,s)}{{I}_{\mathrm{app}}(s)} & = & -\displaystyle \frac{{a}_{{\rm{s}}}^{{\rm{n}}}{{FL}}_{{\rm{n}}}^{2}}{{\kappa }_{\mathrm{eff}}^{{\rm{n}}}}\left(\displaystyle \frac{{{\mathfrak{j}}}_{1}^{{\rm{n}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{n}}}z}-1-z{{\rm{\Lambda }}}_{1}^{{\rm{n}}}\right)}{{\left({{\rm{\Lambda }}}_{1}^{{\rm{n}}}\right)}^{2}}+\displaystyle \frac{{{\mathfrak{j}}}_{2}^{{\rm{n}}}\left({e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{n}}}z}-1+z{{\rm{\Lambda }}}_{1}^{{\rm{n}}}\right)}{{\left({{\rm{\Lambda }}}_{1}^{{\rm{n}}}\right)}^{2}}\right.\\ & & \left.+\ \displaystyle \frac{{{\mathfrak{j}}}_{3}^{{\rm{n}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{n}}}z}-1-z{{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}{{\left({{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}^{2}}+\displaystyle \frac{{{\mathfrak{j}}}_{4}^{{\rm{n}}}\left({e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{n}}}z}-1+z{{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}{{\left({{\rm{\Lambda }}}_{2}^{{\rm{n}}}\right)}^{2}}\right)\\ & & -\ \displaystyle \frac{{\kappa }_{D,\mathrm{eff}}^{{\rm{n}}}}{{\kappa }_{\mathrm{eff}}^{{\rm{n}}}{c}_{{\rm{e}},0}}\left({c}_{1}^{{\rm{n}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{n}}}z}-1\right)+{c}_{2}^{{\rm{n}}}\left({e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{n}}}z}-1\right)\right.\\ & & \left.+\ {c}_{3}^{{\rm{n}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{n}}}z}-1\right)+{c}_{4}^{{\rm{n}}}\left({e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{n}}}z}-1\right)\right).\end{array}\end{eqnarray} \tag{ A·52 }$$

A.6.2. Separator region

In the separator,

$$\begin{eqnarray*}\begin{array}{rcl}\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{e}}}^{{\rm{s}}}(z,s)}{{I}_{\mathrm{app}}(s)} & = & -\displaystyle \frac{{{zL}}_{{\rm{s}}}}{{\kappa }_{\mathrm{eff}}^{{\rm{s}}}A}\\ & & -\ \displaystyle \frac{{\kappa }_{D,\mathrm{eff}}^{{\rm{s}}}}{{\kappa }_{\mathrm{eff}}^{{\rm{s}}}{c}_{{\rm{e}},0}}\left({c}_{1}^{{\rm{s}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{s}}}z}-1\right)+{c}_{2}^{{\rm{s}}}\left({e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{s}}}z}-1\right)\right).\end{array}\end{eqnarray*}$$

A.6.3. Positive electrode

In the positive electrode,

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\tilde{{\rm{\Phi }}}}_{{\rm{e}}}^{{\rm{p}}}(z,s)}{{I}_{\mathrm{app}}(s)} & = & \displaystyle \frac{(z-1){L}_{{\rm{p}}}}{{\kappa }_{\mathrm{eff}}^{{\rm{p}}}A}+\displaystyle \frac{{a}_{{\rm{s}}}^{{\rm{p}}}{{FL}}_{{\rm{p}}}^{2}}{{\kappa }_{\mathrm{eff}}^{{\rm{p}}}}\left(\displaystyle \frac{{{\mathfrak{j}}}_{1}^{{\rm{p}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}z}+(1-z){{\rm{\Lambda }}}_{1}^{{\rm{p}}}{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}-{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}\right)}{{\left({{\rm{\Lambda }}}_{1}^{{\rm{p}}}\right)}^{2}}\right.\\ & & +\ \displaystyle \frac{{{\mathfrak{j}}}_{2}^{{\rm{p}}}\left({e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}z}+(z-1){{\rm{\Lambda }}}_{1}^{{\rm{p}}}{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}\right)}{{\left({{\rm{\Lambda }}}_{1}^{{\rm{p}}}\right)}^{2}}\\ & & +\ \displaystyle \frac{{{\mathfrak{j}}}_{3}^{{\rm{p}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}z}+(1-z){{\rm{\Lambda }}}_{2}^{{\rm{p}}}{e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}-{e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}\right)}{{\left({{\rm{\Lambda }}}_{2}^{{\rm{p}}}\right)}^{2}}\\ & & \left.+\ \displaystyle \frac{{{\mathfrak{j}}}_{4}^{{\rm{p}}}\left({e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}z}+(z-1){{\rm{\Lambda }}}_{2}^{{\rm{p}}}{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}\right)}{{\left({{\rm{\Lambda }}}_{2}^{{\rm{p}}}\right)}^{2}}\right)\\ & & -\ \displaystyle \frac{{\kappa }_{D,\mathrm{eff}}^{{\rm{p}}}}{{\kappa }_{\mathrm{eff}}^{{\rm{p}}}{c}_{{\rm{e}},0}}\left({c}_{1}^{{\rm{p}}}\left({e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}z}-{e}^{{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}\right)+{c}_{2}^{{\rm{p}}}\left({e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}z}-{e}^{-{{\rm{\Lambda }}}_{1}^{{\rm{p}}}}\right)\right.\\ & & \left.+\ {c}_{3}^{{\rm{p}}}\left({e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}z}-{e}^{{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}\right)+{c}_{4}^{{\rm{p}}}\left({e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}z}-{e}^{-{{\rm{\Lambda }}}_{2}^{{\rm{p}}}}\right)\right).\end{array}\end{eqnarray} \tag{ A·53 }$$