Journal of Power Sources 283 (2015) 464e477
Contents lists available at ScienceDirect
Journal of Power Sources journal homepage: www.elsevier.com/locate/jpowsour
Reconstruction of relaxation time distribution from linear electrochemical impedance spectroscopy Yanxiang Zhang a, b, *, Yu Chen b, Mufu Yan a, Fanglin Chen b, * a National Key Laboratory for Precision Hot Processing of Metals, School of Materials Science and Engineering, Harbin Institute of Technology, Harbin, Heilongjiang 150001, China b Department of Mechanical Engineering, University of South Carolina, 300 Main Street, Columbia, SC 29208, USA
h i g h l i g h t s
g r a p h i c a l a b s t r a c t
A method for DRT reconstruction from EIS is presented unambiguously. Quality of DRT reconstruction is sensitive to points per frequency decade. Resolution of DRT is studied to distinguish overlapping processes of SOFCs. This method exhibits high robustness to resist noise embedded in EIS data.
a r t i c l e i n f o
a b s t r a c t
Article history: Received 24 September 2014 Received in revised form 27 November 2014 Accepted 18 February 2015 Available online 26 February 2015
Linear electrochemical impedance spectroscopy (EIS), and in particular its representation of distribution of relaxation time (DRT), enables the identification of the number of processes and their nature involved in electrochemical cells. With the advantage of high frequency resolution, DRT has recently drawn increasing attention for applications in solid oxide fuel cells (SOFCs). However, the method of DRT reconstruction is not yet presented clearly in terms of what mathematical treatments and theoretical assumptions have been made. Here we present unambiguously a method to reconstruct DRT function of impedance based on Tikhonov regularization. By using the synthetic impedances and analytic DRT functions of RQ element, generalized finite length Warburg element, and Gerischer element with physical quantities representative to those of SOFC processes, we show that the quality of DRT reconstruction is sensitive to the sampling points per decade (ppd) of frequency from the impedance measurement. The robustness of the DRT reconstruction to resist noise imbedded in impedance data and numerical calculations can be accomplished by optimizing the weighting factor l according to well defined criterion. © 2015 Elsevier B.V. All rights reserved.
Keywords: Relaxation time distribution (DRT) Impedance Solid oxide fuel cell (SOFC) Electrochemical process Resolution Robustness
1. Introduction Electrochemical impedance spectroscopy (EIS) is widely used to
* Corresponding authors. Department of Mechanical Engineering, University of South Carolina, 300 Main Street, Columbia, SC 29208, USA. E-mail addresses:
[email protected] (Y. Zhang),
[email protected] (F. Chen). http://dx.doi.org/10.1016/j.jpowsour.2015.02.107 0378-7753/© 2015 Elsevier B.V. All rights reserved.
understand and resolve the reaction kinetics at the functional interfacial layer and bulk of electrochemical cells, such as solid oxide fuel cells (SOFCs), gas separation membranes and lithium batteries over a wide frequency range. Impedance of an electrochemical cell is a complex resistance encountered when an AC current perturbation, biejut flows through the cell, and an AC voltage b is generated in response [1]. The causality (analytic perturbation, u
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
b and biejut permits the Taylor expansion of u b relation) between u with respect to biejut [2],
b¼ u
∞ b ðkÞ ∞ k X X u biejut ¼ b v k ejkut k! k¼1
(1)
k¼1
k b is the kth derivative, and b b ðkÞbi =k! is the kth harwhere u vk ¼ u b , which can be expressed as a power series [3], monic of u ðkÞ
b vk ¼
∞ . kþ2r2 X bi i* b h k;kþ2r2
(2)
r¼1
where i* is a real current constant that characterizes the nonlinearity of the cell; the harmonic terms b h k;kþ2r2 are univariate functions of angular frequency u. Therefore, the impedance related to the mth order contribution to the kth harmonic b v k through i* is given by Refs. [3],
Zk;m ðuÞ ¼ b h k;m i*
(3)
One can resolve the impedance only when Eq. (1) converges, b ðkÞ = u b ðkþ1Þ . When the module of bi is yielding bi < lim ðk þ 1Þ u k/∞
b is sufficiently small, only the first-order Taylor expansion of u pronounced and the higher harmonics in Eq. (1) are negligible. In addition, the harmonic terms b h k;m for m > 1 in Eq. (2) are also negligible. Consequently, the impedance becomes [2,3],
ZðuÞ ¼ b h 1;1 i* ¼ b v 1 =bi
(4)
This treatment is known as linearization. The validity can be tested by KramerseKronig transform, which demonstrates that the real component Z0 and imaginary component Z00 are interconvertible when Z(u) is analytic in the upper half-plane of u (indicating causal and stable), and juZ(u) uZ0 (∞)j approaches zero as u approaches infinity (indicating no inductive effect at high frequency) [4]. Nowadays, linearization has become the industry standard in EIS instrumentations, and higher harmonics are filtered out during EIS measurements [5]. A well established method for analyzing electrochemical processes is to fit linearized mechanism models with linear EIS measurements [6]. For example, in the application in SOFCs, the following linearized models are oft-used to characterize the cathodic and anodic processes [7]: 1) RQ element model: It consists of a parallel connection of polarization resistor, R, and a constant phase element, 1=ðjuÞn Q . The impedance is accordingly given by [6],
ZðuÞ ¼
R 1 þ ðjuÞn RQ
(5)
where n is a dimensionless constant in the range from 0 to 1, showing non-ideal double-layer capacitor behavior of the constant phase element; Q is the pre-factor of the constant phase element, with unit of secn/[R]. For n ¼ 1, RQ element degrades into RC element, showing ideal double-layer capacitor behavior. The reciprocal of angular frequency corresponding to the maximum -Z00 is known as the pffiffiffiffiffiffiffi characteristic time, given by n RQ . This model is usually used to characterize the impedance response of surface exchange (or charge transfer reactions) near electrode/electrolyte interface [8e17]. 2) Generalized finite length Warburg (G-FLW) element model: It describes the impedance imparted by mass transport, given by [9],
ZðuÞ ¼ Rw
tanh ðjutw Þnw ðjutw Þnw
465
(6)
where Rw is the polarization resistance of mass transport; nw is a constant approaching but not higher than 0.5; tw is defined as l2/ Deff where l is the finite length and Deff is the effective diffusivity. This model is usually used to characterize the impedance response of gas diffusion within porous electrodes [9e12,17,18]. 3) Gerischer element model: It describes the coupling of surface exchange of oxygen and transport of oxygen ions within porous mixed ionic electronic conductor cathodes, where the resistance of electron transport is negligible. The impedance is given by [19],
Rc ZðuÞ ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ jutc
(7)
where Rc is the polarization resistance of the coupling process; tc is the characteristic time. Rc and tc depend on surface exchange coefficient and oxygen diffusion coefficient of the cathode material, electrode microstructure, and the operating conditions. This model is usually used to characterize the impedance response of oxygen surface exchange kinetics and O2 diffusion in the bulk of (La,Sr) (Co,Fe)O3 (LSCF) based cathode [9e12,17]. Depending on SOFC configuration, fuel cell materials, and working conditions, the mechanism model describing the entire picture of electrochemical processes could be a complicated combination of the above element models or their derivatives [6]. However, it is very challenging to identify the optimal model from many analogs, partially due to the following reasons: 1) Many models fit the EIS data equally well in terms of Nyquist, Real-andImaginary and Bode representations [3]; 2) The impedance responses of electrochemical processes might overlap significantly in frequency domain. In other words, the impedance at an angular frequency u contains remarkable contributions from any process with a characteristic time above and below u1 [20]. Thus, it is of great importance to develop new analysis method with improved frequency resolution. By an extension of the Debye relaxation relation with a single relaxation time, the concept of distribution of relaxation time (DRT) with a high frequency resolution (which will be shown in what follows) was proposed and regarded as the origin of constant phase elements [1]. By assuming a linear relaxation process, the principle of superposition allows the impedance to be generalized [1],
ZðuÞ Z 0 ð∞Þ ¼
Z∞ 0
GðtÞ dt 1 þ jut
(8)
where G(t) is defined as the DRT of impedance Z; t is the relaxation time; Z0 (∞) is the limitation of the real part of Z as angular frequency u approaches infinity. Consequently, the impedance can be represented as series connection of infinite number of circuits made of a parallel connection of polarization resistor G(t)dt and a capacitor t=GðtÞdt, such that G(t)dt specifies the fraction of the total polarization resistance associated with relaxation times in the interval between t and t þ dt. It follows that G(t) is subject to,
Z∞
GðtÞdt ¼ Z 0 ð0Þ Z 0 ð∞Þ
(9)
0
where Z0 (0) Z0 (∞) denotes the total polarization resistance; Z0 (0) is the limitation of the real part of impedance Z as u approaches zero. The unique analytic solution of DRT for a generic impedance spectroscopy was first derived by Fuoss and Kirkwood [21],
Condition
LSC//GDC//LSC symmetric cell; 550 C; PO2 > 0.1 atm Ibid., PO2 ¼ 0.01e0.5 atm Ibid., PO2 ¼ 0.01e0.5 atm Ibid., PO2 ¼ 0.2 atm Ibid., PO2 ¼ 0.2 atm
Ibid. Ibid.
Ibid., 800 C; 28% of CO/CO2 are introduced in anode Ibid. Ibid.
LSM-YSZ//YSZ//Ni-YSZ tubular cell; 700 C; PH2O/PH2 ¼ 2/48e24/26 (50% N2); pure O2 Ibid. Ibid. Ibid. Ni-ScYSZ//YSZ//Ni-ScYSZ symmetric cell; 800 C; PH2O ¼ 3% (balance H2) Ibid. Ibid. LSM-YSZ//YSZ//Ni-YSZ planar cell; 850 C; PH2O/PCO2/PH2 ¼ 45/45/10; pure O2 Ibid. Ibid. Ibid.
RQ
G-FLW
fmax, Hz
Ref.
e
1e5
[8]
0.3e1.77 0.009e0.00015
e e
e 300e2000
[8] [8]
0.002
0.0005
e
e
[8]
0.005
0.0013
e
e
[8]
0.002e0.1
e
0.97
0.3e10
[9,10]
0.001e2
e e
0.89e0.95 0.65e0.75
200e3000 3000e50,000
[9,10] [9,10]
4e10 2e500
[9,10] [9,10]
0.6e0.7
[11,12]
9e10
[11,12]
80
[11,12]
R, Ucm2
C, F cm2
n
Gas diffusion in LSC cathode
<0.004
e
O2 diffusion in LSC bulk Oxygen surface exchange in LSC cathode Oxygen surface exchange at LSC/current collector interface Oxygen surface exchange at LSC/GDC electrolyte interface Gas diffusion in the cathode
0.002e0.003 e
Rw, Ucm2
Gerischer
tw, sec
nw
Rc, Ucm2
tc, sec
Gas diffusion coupled with charge transfer and ionic transport in the anode functional layer Gas diffusion in the anode Oxygen surface exchange kinetics and O2 diffusivity in the bulk of the cathode WGS reaction and CO/CO2 diffusion in the anode supporting layer (ASL) gas diffusion (H2/H2O) in the ASL Oxygen surface exchange kinetics and O2 diffusivity in the bulk of the cathode O2 diffusion in YSZ bulk
0.3
e
e
10,000
[13]
O2 reduction Charge transfer in the anode Gas diffusion in the anode Charge transfer in the anode
1 0.2e0.3 0.1e0.9 0.08
e e e 0.0003
e e e 0.94
10 ~1000 ~0.3 e
[13] [13] [13] [18]
10 e
[18] [18]
Gas diffusion in ASL Gas diffusion in stagnant gas layer High frequency O2 transfer in the cathode
Oxidization at TPB in the anode Reduction at TPB in the cathode Gas diffusion in the anode
0.02e0.15
e
0.45e0.49 0.01e2
0.0044e0.0046
e
e
e
0.28
e
e 0.016e0.02
0.16 0.05
e 0.0105
0.42 0.48
e
0.04e0.06
e
e
30,000
[14]
0.03e0.04
e
e
3000e10,000
[14]
0.01e0.04
e
e
1000
[14]
0.04e0.07
e
e
16e21
[14]
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
LSCF/GDC//YSZ//Ni-YSZ cell; 570e870 C; PH2O ¼ 5.5e65% (balance H2); PO2 ¼ 1e21% (balance N2) Ibid. Ibid.
Physics
466
Table 1 Electrochemical processes derived from DRT analysis at open circuit condition and characterized by RQ, G-FLW and Gerischer circuit elements for various SOFC assemblies and working conditions.
[17] [17] 2e20 10e1000
[17] [17] 200e10,000 6000e80,000
[15] [17] 30,000 0.3e10
[15,16] [15] 700 6.2
[15,16] 5000
1.2e1.9
[14]
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
FðtÞ ¼ tGðtÞlnð10Þ i 00 lnð10Þ h 00 ln tþjp=2 ¼ Z e þ Z eln tjp=2 p
467
(10)
where Z00 (u) denotes the imaginary part of Z; F(t) is subject to,
Z∞
FðtÞd log10 t ¼ Z 0 ð0Þ Z 0 ð∞Þ
(11)
e e
e e
e
e e
e e
e
0.1 0.004
0.021e0.026 0.042e0.044
Ibid. Ibid.
Ibid. LSM-YSZ//YSZ//Ni-YSZ cell; 800 C; PH2/PH2O ¼ 34/63; air Ibid. Ibid.
Anode activation Gas diffusion in the anode O2 transfer in the cathode Gas diffusion in the cathode
LSCF-GDC/GDC//YSZ//Ni-YSZ tubular cell; 600 C; PH2O/PN2/PH2 ¼ 3/57/40; air Ibid. Ibid.
Gas diffusion coupled with charge transfer reaction and ionic transport in the anode functional layer Gas diffusion in the anode Oxygen surface exchange kinetics and O2 diffusivity in the bulk of the cathode
e e 0.09e0.1
0.07e0.08 0.16e0.2
e e Ibid.
Gas conversion in the anode Anode activation
0.12
0.023e0.025
e
e
0.056e0.185
e
∞
The contributions of electrochemical processes become visible, characterized by individual peaks in the F(t) vs. elog10t curve. The solution for F(t) exists when the imaginary component Z00 is an analytic and integrable function of Log10t [21]. The linear EIS yielding KramerseKronig relations is sufficient for the two constraints. However, EIS measurements are typically taken at various discrete frequency points in a defined frequency range, and the analytic relation between the impedance and frequency is unknown. Thus, the DRT function F(t) has to be reconstructed numerically from the discrete EIS data set. One of the first attempts for application in SOFCs was reported by Schichlein et al., who reconstructed the DRT for the impedance of a (La,Sr)MnO3 (LSM) jj Zr0.92Y0.08O1.96 (YSZ) jj Ni-YSZ cell by Fourier deconvolution of Eq. (8), and utilized peak fitting and integration of F(t) vs. elog10t curve to calculate the resistances of individual processes [22]. This method was then applied to characterize other SOFCs, such as La0.6Sr0.4CoO3-d (LSC) jj Ce0.9Gd0.1O1.95 (GDC) jj LSC symmetric cells [8], LSCF j GDC jj YSZ jj Ni-YSZ cells [9e12], Ni-ScYSZ jj YSZ jj NiScYSZ symmetric cells [18], and LSM-YSZ jj YSZ jj Ni-YSZ cells [13,14]. By tailoring operating conditions such as cathodic and anodic atmospheres and temperatures, the electrochemical processes can be differentiated from the reconstructed DRT function due to their different dependence on the operating conditions, providing direct evidence for choosing a proper mechanism model. Table 1 summarizes the DRT analysis results for these cell systems, showing the physical origin of the resolved processes and the corresponding mechanism models represented by RQ elements, GFLW elements, and Gerischer elements, as well as their physical quantities. It has been shown that RQ element could be used to model all the processes, while G-FLW element represents the process of gas diffusion within the function layer and supporting layer of the anode, and the stagnant gas layer above the anode. Gerischer element represents the coupling of oxygen surface exchange kinetics and O2 diffusivity in the bulk of the cathode. The characteristic ordinary frequencies for the processes, fmax ¼ 1/ (2ptmax) corresponding to the peaks of F(t), are found to be in range of 0.1 Hz ~ 105 Hz. Depending on test conditions, the peaks for individual processes may overlap with each other, and even merge into one peak. For example, in the LSCF j GDC jj YSZ jj Ni-YSZ cell, the peak for the coupling of oxygen surface exchange and bulk diffusivity in the LSCF cathode and the peak for gas diffusion in the Ni-YSZ anode supporting layer merge into one peak when temperature decreases from 860 C to 730 C [7]. The expected processes of gas diffusion in pores, charge transfer at three phase boundaries and ionic transport in YSZ in the Ni-YSZ anode functional layer can be only detected as two peeks in the range of testing parameters. The process of gas diffusion in the cathode using air as the oxidant for the cell systems listed in Table 1 is not detectable by DRT due to its low resistance [9]. In addition, DRT function is sensitive to the noise imbedded in the EIS measurements [20]. The frequency resolution and the robustness against noise of DRT reconstruction have not been well understood, imposing a critical concern on the capability of resolving electrochemical processes.
468
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
Fig. 1. DRT solutions derived from Eqs. (10), (17) and (22) and quadratic programming (QP) method for the RQ element (a, R ¼ 0.1 Ucm2, Q ¼ 1 secn/Ucm2, n ¼ 0.97), G-FLW element (b, Rw ¼ 0.0218 Ucm2, tw ¼ 0.0783 s, nw ¼ 0.465), and Gerischer element (c, Rc ¼ 0.0633 Ucm2, tc ¼ 0.016 s). ppd ¼ 10, and l ¼ 1 for Eq. (22) and QP solution.
In this work, a feed-forward calculation procedure to reconstruct the DRT function from discrete EIS data set has been developed. We then study its capability to separate two overlapping processes represented by two RQ elements, G-FLW elements, and Gerischer elements, and test the robustness of reconstruction from noise-free and noisy EIS data set. Finally, practical strategies are provided to enhance the frequency resolution and robustness. 2. Methodology of DRT reconstruction Eq. (10) shows that the DRT function F(t) can be recovered by the imaginary component of impedance. Combining Eq. (10) with Eq. (8), the imaginary component of the impedance Z(u) is given by, 00
Z∞
Z ðuÞ ¼
uFðtÞ=lnð10Þ 2
0
1 þ ðutÞ
00
dt
N X um Fðtn Þ=lnð10Þ n¼1
Z
00
(12)
1 þ ðum tn Þ2
dt
(13)
Wm
10
N X 10Wm wn F 10wn ¼ dw 1 þ 102ðWm wn Þ n¼1
(14)
Therein, dw ¼ jWN W1 =N 1j equals to the reciprocal of frequency points per decade (ppd), 1/ppd. Because ppd is an important setting factor for EIS measurements, its effects on the quality of DRT reconstruction will be studied in what follows. Returning to Eq. (11), the normalization condition can be formulized as, N X
It is noted that the impedance data from the experimental measurement is typically not continuous, but recorded at various logarithmically spaced sampling points of frequency. Thus, F(t) has to be reconstructed numerically from the data set of {un, Z00 (un)} for n ¼ 1, 2, …, N. According to the reciprocity between u and t, the discrete structure of F(t) should be {tn, F(tn)}, where tn ¼ 1/un, for n ¼ 1, 2, …, N. The discretization of Eq. (12) gives,
Z ðum Þ ¼
for m ¼ 1, 2, …, N. As un is logarithmically spaced, we introduce W ¼ log10 u, and w ¼ log10 t into Eq. (13). It follows that,
F 10wn dw ¼ Z 0 10WN Z 0 10W1
(15)
n¼1
Eqs. (14) and (15) provide sufficient information to resolve the discrete structure of F(t). For convenience, the equation system is written in the form of matrix product,
GF ¼ Z
(16)
where G is a matrix with N þ 1 rows and N columns, with Gm;n ¼ dw10Wm wn =1 þ 102ðWm wn Þ for m ¼ 1, 2, …, N, and n ¼ 1, 2, …, N; GNþ1;n ¼ dw for n ¼ 1, 2, …, N; F is a vector with N rows, with Fn ¼ Fð10wn Þ for n ¼ 1, 2, …, N; Z is a vector with N þ 1 rows, 00 with Zn ¼ Z ð10Wn Þ for n ¼ 1, 2, …, N; 00 ZNþ1 ¼ Z ð10WN Þ Z 0 ð10W1 Þ. The objective is to solve F and thus the data set of {tn, F(tn)}. As we have N þ 1 equations and N
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
469
Fig. 2. Imaginary part of impedances for the RQ (a), G-FLW (c), and Gerischer (e) circuit elements as a function of frequency. Scatters in (a, c, e) represent the ideal impedance; Solid lines represent the reconstructed impedance by using the QP solutions of DRT for the various number of logarithmically spaced sampling points per decade (ppd) of frequency. (b), (d), and (f) show respectively the DRT solutions of the RQ, G-FLW, and Gerischer circuit elements. Solid lines with enclosed pattern in (b, d, f) show the analytical solutions of DRT of the corresponding circuit elements; Solid lines show the QP solutions of DRT for the various ppd of frequency. RQ: R ¼ 0.1 Ucm2, Q ¼ 1 secn/Ucm2, n ¼ 0.8; G-FLW: Rw ¼ 0.0218 Ucm2, tw ¼ 0.0783 s, nw ¼ 0.465; Gerischer: Rc ¼ 0.0633 Ucm2, tc ¼ 0.016 s; l ¼ 1.
unknowns, the equation system is overdetermined. Thus, the solution of F can be given by ordinary least squares method,
F ¼ ðGT GÞ1 GT Z
(17)
where the superscript T represents transposition. However, the solution is highly sensitive to G and Z, since G is ill conditioned [23]. Mathematically, some regularization term should be introduced to
470
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
Fig. 3. Imaginary impedance of a series circuit consisting of two RQ elements for R1 ¼ 0.1 Ucm2, Q1 ¼ 1 secn/Ucm2, n1 ¼ 0.8 and R2 ¼ 0.1 Ucm2, Q2 ¼ 5 secn/Ucm2, n2 ¼ 0.8 (a). Scatters in (a) represent the ideal impedance; Solid line represents the reconstructed impedance by using the QP solution of DRT for the 10, 50 and 100 ppd of frequency. (b) shows the analytical solution of DRT (black line with enclosed pattern), and the QP solutions of DRT for the 10, 50 and 100 ppd of frequency (l ¼ 1) of the series circuit element.
make the solution robust. Physically, the vector F behaves mostly analytic. In other words, the first-order derivative of F with respect to ew is continuous, indicating the following regularization,
DF ¼ 0
Zerr
(19)
where the first part relates to the standard deviation of impedance,
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðGF ZÞT ðGF ZÞ ¼ N1
(20)
The second part in Eq. (19) relates to the standard deviation of F,
(18)
where D is a N factorial square matrix, with Dn, [n-1, n, nþ1] ¼ {1, 2, 1} for n ¼ 2, 3, …, N 1; D1, [1,2] ¼ {1, 1}; DN, [N-1, N] ¼ {1, 1}; and the other components of D is zero. This is known as Tikhonov regularization, one of the most widely used methods to solve ill conditioned problems [24]. In fact, Sonn et al. have for the first time utilized Tikhonov regularization to reconstruct the DRT function of hydrogen oxidation in a Ni/YSZ anode [25]. However, they did not report more information about the method. This study intends to clarify what mathematical treatments and theoretical assumptions are made to reconstruct the DRT by the linear EIS data. Now the objective is to solve F that minimizes the following residual,
err ¼ ðGF ZÞT ðGF ZÞ þ lðDFÞT ðDFÞ
Fig. 4. Contour maps of the resolution threshold of ppd above which two RQ elements in series are resolvable from the QP solution of DRT for l ¼ 1, in the plane of logarithmical ratios of ohmic resistance and equivalent capacity of the two RQ elements. (a) shows the contour map for baseline R1 ¼ 0.05 Ucm2, Q1 ¼ 0.0386 secn/Ucm2, n1 and 2 n 2 ¼ 0.97; (b) shows the contour map for baseline R1 ¼ 0.5 Ucm , Q1 ¼ 0.0386 sec / Ucm2, n1 and 2 ¼ 0.97; (c) shows the contour map for baseline R1 ¼ 0.05 Ucm2, Q1 ¼ 0.0386 secn/Ucm2, n1 and 2 ¼ 0.65; (d) shows the contour map for baseline R1 ¼ 0.05 Ucm2, Q1 ¼ 0.0386 secn/Ucm2, n1 ¼ 0.97, n2 ¼ 0.65. When the two RQ elements fall into the dark region, they are theoretically not resolvable from the analytical solution of DRT.
Ferr
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðDFÞT ðDFÞ ¼ N1
(21)
l in Eq. (19) is a non-negative weighting factor that balances the goodness-of-fit of impedance and the smoothness of F. The determination of optimal l will be given in what follows, but l ¼ 1 is a good choice. The minimization of Eq. (19) suggests derr/dF ¼ 0, and thus, F ¼ ðGT G þ lDT DÞ1 GT Z
(22)
When l ¼ 0, Eq. (22) degrades into Eq. (17). However, by using Eq. (22), sometimes the elements in the solution of F could be negative for the RQ elements, G-FLW elements, and Gerischer elements. This is not physically rational, since their analytic solutions are non-negative. Although no guarantee exists that F for generic impedance is non-negative, the inductive effect is not usually
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
471
Fig. 5. Contour maps of the resolution threshold of ppd above which two G-FLW elements in series are resolvable from the QP solution of DRT for l ¼ 1, in the plane of logarithmical ratios of ohmic resistance and characteristic time of the two G-FLW elements. (a) shows the contour map for baseline Rw1 ¼ 0.02 Ucm2, tw1 ¼ 0.01 s, nw1 and 2 ¼ 0.49; (b) shows the contour map for baseline Rw1 ¼ 0.02 Ucm2, tw1 ¼ 0.01 s, nw1 and 2 ¼ 0.42; (c) shows the contour map for baseline Rw1 ¼ 0.02 Ucm2, tw1 ¼ 0.01 s, nw1 ¼ 0.49, nw2 ¼ 0.42. When the two G-FLW elements fall into the dark region, they are theoretically not resolvable from the analytical solution of DRT.
3. Results and discussion 3.1. Performance of Eqs. (17) and (22), and quadratic programming method
Fig. 6. Contour map of the resolution threshold of ppd above which two Gerischer elements in series are resolvable from the QP solution of DRT for l ¼ 1, in the plane of logarithmical ratios of ohmic resistance and characteristic time of the two Gerischer elements. The baseline is given as Rc1 ¼ 0.1 Ucm2, tc1 ¼ 0.005 s.
encountered in electrochemistry [5]. Thus, the following constraint may be rational,
F0
(23)
Now the objective is to solve F that minimizes Eq. (19), subjecting to Inequality 23. This is known as a quadratic programming (QP) problem [23]. The solution of the QP problem is unique since GTG þ lDTD is an invertible matrix. Commercial software is available to solve the QP problem. The software used here is the quadprog function built in Matlab package.
To illustrate the performance of Eqs. (17) and (22), and QP method, Fig. 1a, b, and c shows respectively their DRT reconstruction results of a RQ element (a, R ¼ 0.1 Ucm2, Q ¼ 1 secn/Ucm2, n ¼ 0.97), G-FLW element (b, Rw ¼ 0.0218 Ucm2, tw ¼ 0.0783 s, nw ¼ 0.465), and Gerischer element (c, Rc ¼ 0.0633 Ucm2, tc ¼ 0.016 s). We choose the frequency range within which the absolute value of the imaginary impedance is higher than 0.1% of the total resistance of the model element. The ppd is given as 10, typical for EIS measurements. The DRT function F(t), defined as tG(t)ln(10) is plotted versus eLog10(2pt). F(t) has the same unit of polarization resistance, Ucm2; eLog10(2pt) is dimensionless, standing for the exponent of the ordinary frequency, 1=2pt measured in Hz. In this representation, the area enclosed by the curve stands for the polarization resistance. For validation, the analytic DRT solutions for the three elements are plotted in Fig. 1a, b, and c, using Eq. (10). As shown in Fig. 1a, b, and c, Eq. (17) is not capable of reconstructing the DRT function, due to the illconditioned nature of matrix G. Eq. (22) shows better performance in comparison to Eq. (17). The characteristic ordinary frequency corresponding to the peak value of DRT function and the polarization resistance corresponding to the integral area of DRT function by Eq. (22) are identical to that of the analytic solution by
472
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
Fig. 7. Imaginary impedance of a RQ element for R ¼ 0.1 Ucm2, Q ¼ 1 secn/Ucm2, and n ¼ 0.8 (a), showing the ideal impedance (scatter), and the reconstructed impedance by using the QP solutions of 10 ppd DRT for l ¼ 0, 1 and 10. (b) shows the analytical solution of DRT (black line with enclosed pattern), and the QP solutions of 10 ppd DRT for l ¼ 0, 1 and 10 of the RQ element.
Eq. (10). But there exist some negative components in the solution of Eq. (22), which is not physically rational. By introducing the constraint that DRT function is non-negative, the QP solution avoids negative values, showing improved quality compared to Eq. (22). Thus, it is more feasible to use the QP method to reconstruct DRT functions.
element (R ¼ 0.1 Ucm2, Q ¼ 1 secn/Ucm2, n ¼ 0.8), G-FLW element (Rw ¼ 0.0218 Ucm2, tw ¼ 0.0783 s, nw ¼ 0.465), and Gerischer element (Rc ¼ 0.0633 Ucm2, tc ¼ 0.016 s). The physical quantities of these circuit elements approach those of SOFCs in Table 1. The analytic DRT functions for the RQ element, G-FLW element, and Gerischer element are shown as solid lines with enclosed patterns in Fig. 2b, d, and f, respectively. The thick solid lines in Fig. 2b, d, and f show the reconstructed DRT functions by using 10 ppd, 50 ppd, and 100 ppd. As shown in Fig. 2b, the reconstructed DRT approaches the ideal one for the RQ element by increasing ppd value. The broadening effect is eliminated when ppd is above 50. By using Eq. (16), the imaginary component of the RQ impedance can be backward recalculated from the reconstructed DRT function, as the thick solid lines shown in Fig. 2a. It is shown that the recalculated impedances for the various ppd values fit the ideal impedance equally well, indicating that the reconstruction result of DRT function is sensitive to ppd value, but different DRT functions can lead to nearly identical impedance. This result suggests that a test of ppd value may be necessary for EIS measurement to get a satisfactory reconstruction of DRT function. Fig. 2c and d show the results for the G-FLW element. The recalculated impedance for 10 ppd does not show a perfect fitting to the ideal impedance, due to the broadening effect shown in Fig. 2d. The reconstructed DRT function for 10 ppd exhibits only one broad peak. But the analytic DRT function for the G-FLW element exhibits multiple peaks, characterized by a large peak followed by smaller peaks at higher frequencies. By increasing ppd to 50 and 100, the reconstructed DRT function shows two characteristic peaks with narrowed half-width. As shown in Fig. 2c, the recalculated impedances for 50 ppd and 100 ppd agree nicely with the ideal impedance. Fig. 2e and f show similar results for the Gerischer element compared to the G-FLW element. The DRT function for 10 ppd also exhibits one broad peak, and the shape is very similar to that of the G-FLW element. In other words, this may be challenging to distinguish between G-FLW element and Gerischer element when ppd is insufficient. Actually, the analytic DRT function for G-FLW element is very different from that for Gerischer element. As shown in Fig. 2f, there exists a threshold of Log10(2ptc), below which F is zero, and above which F decreases from infinity to zero with Log10(2pt) increasing. By increasing ppd to 50 and 100, the reconstructed DRT function is sharpened, and the threshold at Log10(2ptc) becomes more clear. But a tendency of oscillating declination above Log10(2ptc) is observed. This is induced by the regularization of Eq. (18), because the DRT function at Log10(2ptc) is not analytic. The same problem is encountered by using Fourier deconvolution method [20], where the lowpass filters have the same effects as Eq. (18). The oscillation seems to be unavoidable, because the regularization term is indispensable for DRT reconstruction. 3.3. Frequency resolution of DRT reconstruction
3.2. Points-per-decade dependence of DRT reconstruction It is worth noting that the half-width of the QP solution shown in Fig. 1 is broader than that of the analytic solution. The broadening is caused by the discrete nature of the impedance data, which is also observed by the Fourier deconvolution method [7]. It is natural to consider that increasing the sampling points per decade of frequency makes the impedance data more comprehensive, and thus may suppress the broadening effect. However, the EIS measurements usually employ a typical ppd value such as 10. Effects of ppd on the quality of DRT reconstruction are not clear. This work will verify that the broadening effects can be suppressed or even be eliminated by increasing ppd value, as shown in Fig. 2. The scatters in Fig. 2a, c, and e show respectively the ideal imaginary component as a function of ordinary frequency for a RQ
The frequency resolution of DRT reconstruction is defined by the solution of the following problem: given the impedance of two overlapping processes, find the threshold ppd value above which the two processes can be distinguished from the reconstructed DRT function; otherwise, find the threshold degree to which the two processes are overlapped that the two processes become indistinguishable from the analytic DRT function. A high frequency resolution means a low threshold of ppd value or a high threshold of overlapping degree. Refocusing on the impedance of the RQ element, G-FLW element and Gerischer element shown in Fig. 2a, c and e, and their reconstructed DRT functions shown in Fig. 2b, d and f, the half-width of DRT representation is significantly narrower than that of the Real-and-Imaginary representation, indicating a higher frequency resolution of DRT representation. In
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
473
Fig. 8. Computational errors of 10 ppd impedance (Zerr), DRT (Ferr) and their summation as a function of l for the RQ element (a, R ¼ 0.1 Ucm2, Q ¼ 1 secn/Ucm2, n ¼ 0.8), G-FLW element (b, Rw ¼ 0.0218 Ucm2, tw ¼ 0.0783 s, nw ¼ 0.465), and Gerischer element (c, Rc ¼ 0.0633 Ucm2, tc ¼ 0.016 s).
addition, the reconstructed DRT function approaches the ideal one while being more sharpened by increasing ppd value, suggesting that the frequency resolution can be rationally improved by increasing ppd value. Fig. 3 shows an illustrative example to distinguish two processes made of a series connection of RQ elements for R1 ¼ 0.1 Ucm2, Q1 ¼ 1 secn/Ucm2, n1 ¼ 0.8 and R2 ¼ 0.1 Ucm2, Q2 ¼ 5 secn/Ucm2, n2 ¼ 0.8. Fig. 3a shows the imaginary component of impedance versus frequency, showing that both the recalculated impedances based on DRT reconstruction for 10, 50 and 100 ppd agree nicely with the ideal impedance. The two processes merge into one broad peak in Real-and-Imaginary representation, and thus cannot be separated yet. However, as shown in Fig. 3b, the analytic DRT function for the combined process exhibits two peaks, in other words, the two processes are theoretically distinguishable in the DRT representation. By numerical reconstruction, the DRT function for 10 ppd shows only one peak, while the results for 50 and 100 ppd show two peaks. Thus, the two processes can be separated by DRT reconstruction, and the threshold ppd value is between 10 and 50. In what follows, the capability of separating two overlapping processes with physical quantities approaching those of operating SOFCs will be tested. First, it is necessary to define the two overlapping processes. As listed in Table 1, the physical quantities of the electrochemical processes characterized by RQ elements fall into the range: fmax ¼ 0.3e8 104 Hz; R ¼ 0.001e2 Ucm2;
n ¼ 0.65e0.97. The quantities by G-FLW elements fall into the range: fmax ¼ 2e20 Hz; Rw ¼ 0.02e0.28 Ucm2; nw ¼ 0.42e0.49. The quantities by Gerischer elements fall into the range: fmax ¼ 2e1000 Hz; Rc ¼ 0.01e2 Ucm2. Accordingly, the two overlapping processes can be defined as shown in Figs. 4e6, where one process marked as subscript 1 is the reference process and the other one marked as subscript 2 is the independent process. By tailoring the independent process in the given range, the frequency resolution for a given overlapping degree is calculated. It should be noted that the ppd value for practical EIS measurements should not be extremely high, otherwise the measurements may take a long period of time, and the KramerseKronig relation between the real part and imaginary part of impedance may not be valid. In this regard, the upper limit of ppd in this work is given as 100. Fig. 4 shows the map of frequency resolution for two overlapping RQ processes. Therein, the isolines of ppd threshold are plotted in the coordinates of Log10(R2/R1) and Log10(C2/C1). C is the equivalent capacity, defined as R(1n)/nQ1/n. The dark region represents that the two processes are theoretically indistinguishable if the physical quantities fall into this region. In Fig. 4a, the reference process is given as R1 ¼ 0.05 Ucm2, Q1 ¼ 0.0386 secn/Ucm2, n1 and 2 ¼ 0.97. It is shown that the isolines of threshold ppd are basically in parallel with that of Log10(R2/R1) þ Log10(C2/C1). When the reference R1 changes from 0.05 Ucm2 to 0.5 Ucm2, the resolution map shown in Fig. 4b is not changed. In other words, the threshold
474
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
Fig. 9. Computational errors of 10 ppd impedance (Zerr), DRT (Ferr) and their summation as a function of l for the RQ element (a, R ¼ 1 Ucm2, Q ¼ 1 secn/Ucm2, n ¼ 0.8), G-FLW element (b, Rw ¼ 0.218 Ucm2, tw ¼ 0.0783 s, nw ¼ 0.465), and Gerischer element (c, Rc ¼ 0.633 Ucm2, tc ¼ 0.016 s).
ppd is governed by the difference of characteristic time between the two RQ processes. Increasing ppd can increase frequency resolution. For instance, 4 ppd can resolve two processes within two decades of characteristic time, while 6 ppd increases the resolution to 1.5 decades. But there is not apparent improvement in resolution when ppd increases from 50 to 100. There exists a gap of 0.5 time decade that demands ppd higher than 100. Fig. 4c shows the resolution map when the reference n1 and 2 decreases from 0.97 to 0.65, showing that the dark region expands substantially with the shift and distortion of isolines, because the half-width of DRT function for RQ elements increases with the decrease of factor n. Fig. 4d shows the resolution map when the two overlapping RQ processes have different values of n, given as n1 ¼ 0.97 and n2 ¼ 0.65. It is shown that the dark region is smaller compared to Fig. 4c. Isolines of ppd threshold are not shifted but distorted a little compared to Fig. 4a. Fig. 5 shows the map of frequency resolution for two overlapping G-FLW processes. The isolines of ppd threshold are plotted in the coordinates of Log10(Rw2/Rw1) and Log10(tw2/tw1). Rw1 and tw1 are given as 0.02 Ucm2, and 0.01 s. To study the effects of factor nw, three cases are studied, given as nw1 and 2 ¼ 0.49 for Fig. 5a, nw1 and 2 ¼ 0.42 for Fig. 5b, and nw1 ¼ 0.49, nw2 ¼ 0.42 for Fig. 5c. Similar to RQ elements, the ppd threshold is governed by the difference of characteristic time between the two G-FLW elements. Although the dark region expands by decrease of nw, the isolines of ppd threshold
are not changed apparently. 20 ppd can resolve two G-FLW elements within one decade of characteristic time. But it seems challenging to resolve G-FLW elements within 0.6 decade. Fig. 6 shows the map of frequency resolution for two overlapping Gerischer elements. The isolines of ppd threshold are plotted in the coordinates of Log10(Rc2/Rc1) and Log10(tc2/tc1). The reference process is given as Rc1 ¼ 0.1 Ucm2 and tc1 ¼ 0.005 s. The threshold ppd is also governed by the characteristic time. Because its DRT function is singular at the characteristic time, any two Gerischer elements with different characteristic times are theoretically distinguishable. Thus, the dark region is a vertical line at Log10(tc2/tc1) ¼ 0 as shown in Fig. 6. But there exists a gap of 0.6 decade, that demands ppd higher than 200. 3.4. Robustness of DRT reconstruction The robustness of DRT reconstruction is the stability of the numerical algorithm and its ability to resist the noise imbedded in impedance. This could be well controlled by the weighting factor l in Eq. (19). Fig. 7 shows the effects of l on the quality of DRT reconstruction for the impedance defined in Fig. 2a. It is shown in Fig. 7a that the recalculated impedance by the reconstructed DRT function for l ¼ 0 agrees better with the ideal impedance, compared to that for l ¼ 10. However, as shown in Fig. 7b, the reconstructed DRT function for l ¼ 0 exhibits artificial peaks on
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
475
Fig. 10. Computational errors of 100 ppd impedance (Zerr), DRT (Ferr) and their summation as a function of l for the RQ element (a, R ¼ 0.1 Ucm2, Q ¼ 1 secn/Ucm2, n ¼ 0.8), G-FLW element (b, Rw ¼ 0.0218 Ucm2, tw ¼ 0.0783 s, nw ¼ 0.465), and Gerischer element (c, Rc ¼ 0.0633 Ucm2, tc ¼ 0.016 s).
shoulder of the major peak, which do not exist in the analytic DRT function. The peak artifacts disappear for l ¼ 10, but the DRT function is over smoothed, leading to bad agreement in impedance. Good agreements in both impedance and DRT function are shown for l ¼ 1. Thus, there should be an optimal choice of l. A criterion is necessary to determine the optimal l value. To this end, the dependence of Zerr (Eq. (20)) and Ferr (Eq. (21)) on factors such as the nature of process, absolute value of impedance, value of ppd and value of l has been studied. Fig. 8 shows the effects of the nature of the process, revealing evolutions of standard deviation of impedance, Zerr, standard deviation of DRT function, Ferr, and their summation as a function of l, for the RQ element, G-FLW element, and Gerischer element defined in Fig. 2. It can be seen that with the increase of l, Zerr increases while Ferr decreases. By inspection of Fig. 8, some general criteria can be obtained to estimate the optimal l. One criterion is to use the pattern of Ferr. For example in Fig. 8b and c, there is an oscillation in the Ferr. vs. l curve when l < 0.2, indicating that artificial peaks of DRT function may be generated; when l > 0.2, DRT function tends to be over smoothed because the increase of Zerr. Thus l ¼ 0.2 could be used as the optimal value. However, the threshold value of l above which the oscillation is eliminated depends on nature of process. As shown in Fig. 8a for the RQ element, the optimal l is smaller than 0.01. This criterion may exhibit high frequency resolution due to the significantly low value of Zerr, but particular attention is needed for specific cases. Another
criterion is the summation of Zerr and Ferr. The default viewpoint of this criterion is that Zerr and Ferr contribute equally to the quality of DRT reconstruction. As shown in Fig. 8, there exists a minimal summation of Zerr and Ferr. The corresponding values of l for these elements approach to 1, but not exactly. Without losing generality, we may use the weighted summation, such as Zerr þ cFerr, where c is a weighting factor. For c < 1, the reconstructed DRT function exhibits smaller value of Zerr, with the sacrifice of smoothness. The choice of c depends on one's preference. In this work, Zerr þ Ferr is employed as the criterion. To illustrate the feasibility of this criterion, effects of absolute value of impedance and value of ppd on the evolution of the residual terms are shown respectively in Figs. 9 and 10. Compared to the RQ element, G-FLW element, and Gerischer element in Fig. 8, the absolute values of impedance are magnified ten times. As shown in Fig. 9, both the residual terms, Zerr and Ferr are increased to a level of ten times for the RQ element, G-FLW element, and Gerischer element. The increase of absolute value of impedance shows little effects on the shape of residual vs. l curve, and thus the optimal value of l. The interplay between ppd and l is examined by increase the ppd value from 10 to 100. As shown in Fig. 10, the increase of ppd decreases both the values of Zerr and Ferr to the same degree. For the RQ element, G-FLW element, and Gerischer element, the residuals are decreased to a level of 50 times, 20 times, and 2 times, respectively. This is the reason why the reconstructed DRT function with high ppd value could yield good
476
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
Fig. 11. Ideal and noisy impedance spectrums with 10 ppd of a RQ element for R ¼ 0.1 Ucm2, Q ¼ 1 secn/Ucm2, n ¼ 0.8 (a), and the computational errors of impedance (Zerr), DRT (Ferr) and their summation as a function of l for the noisy impedance (b), and the QP solutions of DRT for the ideal and noisy impedance spectrums (c).
approximation of both the analytical DRT and ideal impedance. The interplay between ppd and l is reflected by the oscillation of residual vs. l curve, however, shows little effects on the optimal l. These results show little dependence of the optimal l on the nature of process, absolute value of impedance and value of ppd. Thus, this criterion permits successful isolation between ppd and l. As derived by the interplay between l and other factors, the optimal value of l ranges from 0.5 to 2. As the change of residuals is not pronounced within this range of l, the optimal value of l is fixed as 1 to study the effects of ppd in the above section. However, if we use some other criterion depending on the oscillation of residual vs. l curve, the interplay between the nature of process, absolute value of impedance, value of ppd, and value of l will be complex. By using this criterion, the DRT reconstruction is found to be very robust when the impedance contains random noise. A representative example is shown in Fig. 11. The ideal impedance and the noisy impedance with 10 ppd of a RQ element for R ¼ 0.1 Ucm2, Q ¼ 1 secn/Ucm2, n ¼ 0.8 are studied as shown in Fig. 11a. For the ideal impedance, the reconstructed DRT function by l ¼ 1 is shown in Fig. 11c. For the noisy impedance, the optimal l value can be given as 2, corresponding to the minimal Zerr þ Ferr as shown in Fig. 11b. For comparison with the reconstructed DRT function of the ideal impedance, Fig. 11c shows the reconstructed DRT functions of
the noisy impedance by l ¼ 0.01, 2 and 10, showing that the one by l ¼ 2 is nearly identical to the DRT function of the ideal impedance. 4. Conclusions The distribution of relaxation time (DRT) can be reconstructed numerically from discrete EIS data set by using quadratic programming method. The feasibility of the method is verified by the impedance and analytic DRT functions of RQ, G-FLW and Gerischer elements. The frequency points per decade (ppd) for the EIS data and the weighting factor (l) in Eq. (19) are found to be important factors governing the frequency resolution and robustness of the DRT reconstruction. It is concluded that: 1) The quality of DRT reconstruction is sensitive to ppd value. By increasing ppd value, the reconstructed DRT function shows improved consistence with the analytic DRT function, while the half-width being more sharpened. 2) Two processes that overlap into one peak in the Real-andImaginary representation of impedance may be distinguishable in the DRT representation. The frequency resolution of DRT representation can be extensively improved by increasing ppd value.
Y. Zhang et al. / Journal of Power Sources 283 (2015) 464e477
3) By tailoring the value of l, the peak artifacts in reconstructed DRT function can be eliminated, and the effects of noise imbedded in impedance can be filtered out. Acknowledgments We gratefully acknowledge the financial support from the U.S. National Science Foundation (DMR-1210792) and Natural Science Foundation of China (51402066). We are grateful to Dr. Kyle N. Grew from the U.S. Army Research Laboratory who has read the manuscript and provided valuable feedback to improve this study. References [1] E. Barsoukov, J. Ross Macdonald (Eds.), Impedance Spectroscopy Theory, Experiment, and Applications, second ed., John Wiley & Sons, Inc., Hoboken, New Jersey, 2005. [2] Q.A. Huang, R. Hui, B. Wang, J. Zhang, Electrochim. Acta 52 (2007) 8144e8164. [3] J.R. Wilson, D.T. Schwartz, S.B. Adler, Electrochim. Acta 51 (2006) 1389e1402. [4] J.M. Esteban, M.E. Orazem, J. Electrochem. Soc. 138 (1991) 67e76. [5] S.M. Park, J.S. Yoo, Anal. Chem. 75 (2003) 455 Ae461 A. [6] B. Tribollet, M.E. Orazem, Electrochemical Impedance Spectroscopy, John Wiley & Sons, Inc., Hoboken, New Jersey, 2008. [7] A. Leonide, SOFC Modelling and Parameter Identification by Means of Impedance Spectroscopy (PhD. dissertation), Karlsruher Institute of Technology, 2010. e, J. Electrochem. Soc. 160 (2013) F1197eF1206. [8] J. Hayd, E. Ivers-Tiffe e, J. Electrochem. Soc. 155 (2008) [9] A. Leonide, V. Sonn, A. Weber, E. Ivers-Tiffe B36eB41.
477
e, [10] A. Leonide, B. Rüger, A. Weber, W.A. Meulenberg, E. Ivers-Tiffe J. Electrochem. Soc. 157 (2010) B234eB239. e, J. Electrochem. Soc. [11] A. Kromp, S. Dierickx, A. Leonide, A. Weber, E. Ivers-Tiffe 159 (2012) B597eB601. e, Electrochim. Acta 106 (2013) [12] A. Kromp, H. Geisler, A. Weber, E. Ivers-Tiffe 418e424. [13] B. Liu, H. Muroyama, T. Matsui, K. Tomida, T. Kabata, K. Eguchi, J. Electrochem. Soc. 157 (2010) B1858eB1864. [14] Y. Tao, S.D. Ebbesen, M.B. Mogensen, J. Electrochem. Soc. 161 (2014) F337eF343. [15] H. Sumi, T. Yamaguchi, K. Hamamoto, T. Suzuki, Y. Fujishiro, J. Power Sources 226 (2013) 354e358. [16] H. Sumi, T. Yamaguchi, K. Hamamoto, T. Suzuki, Y. Fujishiro, T. Matsui, K. Eguchi, Electrochim. Acta 67 (2012) 159e165. e, [17] M. Kornely, A. Neumann, N.H. Menzler, A. Leonide, A. Weber, E. Ivers-Tiffe J. Power Sources 196 (2011) 7203e7208. [18] T. Ramos, M. Søgaard, M.B. Mogensen, J. Electrochem. Soc. 161 (2014) F434eF444. [19] S.B. Adler, J.A. Lane, B.C.H. Steele, J. Electrochem. Soc. 143 (1996) 3554e3564. e, J. Appl. Elec[20] H. Schichlein, A.C. Müller, M. Voigts, A. Krügel, E. Ivers-Tiffe trochem. 32 (2002) 875e882. [21] R.M. Fuoss, J.G. Kirkwood, J. Am. Chem. Soc. 63 (1941) 385e394. e, [22] H. Schichlein, M. Feuerstein, A. Müller, A. Weber, A. Krügel, E. Ivers-Tiffe System identification: a new modelling approach for SOFC single cells, in: S.C. Singhal, M. Dokiya (Eds.), Solid Oxide Fuel Cells (SOFC VI), The Electrochemical Society, Inc., Pennington NJ, 1999, pp. 1069e1077. [23] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes the Art of Scientific Computing, third ed., Cambridge University Press, New York, 2007. [24] A.N. Tikhonov, V.Y. Arsenin, Solutions of Ill-posed Problems, 1977. Winston, Washington, DC. e, J. Electrochem. Soc. 155 (2008) B675eB679. [25] V. Sonn, A. Leonide, E. Ivers-Tiffe