Journal of Power Sources 449 (2020) 227499

Contents lists available at ScienceDirect

Journal of Power Sources

journal homepage: www.elsevier.com/locate/jpowsour

Prediction of overpotential and concentration profiles in solid oxide fuel

cell based on improved analytical model of charge and mass transfer
Daili Feng, Cheng Bao *, Ting Gao
School of Energy and Environmental Engineering, University of Science and Technology Beijing, Beijing, 100083, PR China

H I G H L I G H T S

� A linear translation (LT) greatly improves the conventional linear approach.

� Improved power-law (PL) nonlinear approximation combined with the LT approach.
� Analytical solutions of decoupled charge and mass transfer.
� Validation of LT þ PL approach in a wide range of SOFC operations.
� The closed-form approaches are applicable for general electrochemical systems.

A R T I C L E I N F O A B S T R A C T

Keywords: The accuracy of prediction of overpotential and relevant concentration distributions is important for under­
Solid oxide fuel cell standing fundamentals and electrochemical reaction engineering. The conventional linear approximation (CLA)
Butler-Volmer kinetics of the exponent items in the Butler-Volmer kinetics causes significant discrepancies in the estimation of the
Analytical modeling
overpotential profile at moderate and high electric loads. Based on a one-dimension dimensionless charge
Overpotential
Linear translation
transfer model, an approach of linear translation (LT) is presented to solve the problem of the CLA and improve
Mass transfer the power-law (PL) approach of nonlinear approximation (J Power Sources, 210 (2012) 67–80), which is ac­
curate at high current densities but is still limited to cases of thick electrodes. The hybrid LT þ PL approach is
then convincingly validated in a wide range of fuel cell operations, including moderate and high operating
current densities, electrolyte- and electrode-supported electrodes, and asymmetric anodic and cathodic transfer
coefficients. Then the analytical species concentration profiles are obtained from a decoupled charge/mass
transfer model. The closed-form electrode-level model significantly improves the computational efficiency in
thermo-electrochemical cell-level simulation. Although discussed in the context of solid oxide fuel cells, the
approaches in this work are applicable to general electrochemical systems which hold Butler-Volmer kinetics.

1. Introduction while in the introduction of the book, it makes a declaration that all
analytical models presented can be equally applied to solid oxide fuel
An accurate solution of the electrochemical reaction kinetics cells (SOFCs) since the similar working principles. Moreover, Kulikov­
involved in transport phenomena is essential to the effective prediction sky’s analytical solution was also adopted to SOFC by Prokop et al. [2],
of the profiles of current densities, overpotential and concentration and who also modified the solution to include the microstructure of the
of the design of fuel cells and system. Compared to time-cost numerical electrode. Bao developed an analytical modeling framework of elec­
computations, exact solutions or analytical approximations can provide trochemistry and transport in solid oxide fuel cells (SOFCs) with his
a better understanding of the native of physicochemical phenomena in specific work [3,4].
closed form and more efficient calculations, which are valuable for The Butler-Volmer (BV) equation is one of the standard formulations
system-level optimization and real-time control. An overview of describing the global or/and elementary reaction kinetics in electro­
analytical models of fuel cells was seen in Kulikovsky’s book [1], which chemical systems. The BV equation was originally derived for electro­
is mainly focused on polymer electrolyte membrane fuel cells (PEMFCs), chemical reactions with a single charge transfer as the single rate-

* Corresponding author.
E-mail address: [email protected] (C. Bao).

https://doi.org/10.1016/j.jpowsour.2019.227499
Received 24 July 2019; Received in revised form 18 November 2019; Accepted 22 November 2019
Available online 29 November 2019
0378-7753/© 2019 Elsevier B.V. All rights reserved.
D. Feng et al. Journal of Power Sources 449 (2020) 227499

limiting step at metal-solution interfaces with the variation of potential

across a double layer [5]. The electrochemical kinetics is generally
believed to be in BV formulation in PEMFCs in which the metal catalysts
are effectively immersed in water to facilitate the transfer of proton into
the electrolyte. Although there are some arguments on the applicability
of BV kinetics in solid-oxide systems [6] and some non-BV kinetics
presented in SOFC modeling [7–11], the BV equation was proved,
especially by the fundamental research work performed by the group of
Prof. Ivers-Tiff�ee, to be a very valid approximation in SOFC, by elegantly
extracting the correct kinetic and physically meaningful parameters
from electrochemical measurements (e.g. by using impedance spec­
troscopy measurements) [12–15]. As extensively used in SOFC
modeling, the simplicity and generality of the BV equation provides a
good starting formulation for both fundamental analysis and fitting of
electrode kinetics in engineering applications [3,16–19]. The elemen­
tary electrochemistry modeling (not limited to SOFCs) could also
generate a global BV kinetics, although in some cases there might be
more complex mathematical formulations, e.g. a more nonlinear
dependence of the exchange current density on species concentrations
[5,20].
However, the nonlinearity of exponent items in the BV equation
brings difficulties to obtain the explicitly analytical overpotential as a
function of current density, which is normally required in the prediction
of a voltage-current performance of fuel cells. The BV equation expresses
the electrochemical reaction rate as j ¼ j0[exp(αneFη/RT) exp( βneFη/
RT)], in which η is the overpotential, j0 is the volume-specific (with unit
A m 3) or area-specific (A m 2) exchange current density, ne is the
number of electrons involved in the reaction, and α and β are the anodic
and cathodic transfer coefficient, respectively, and can be normally
taken to be symmetric as α ¼ β ¼ 0.5 for an ideal elementary charge
transfer step. Chan et al. [21] introduced an algebraic calculation of a
total overpotential ηt ¼ RT/αneF.sinh 1(j/2j0) for area-specific modeling
when α ¼ β. For volume-specific charge transfer modeling, differential
equations are needed (cf. Section 2), the overpotential distribution is
significantly influenced by the magnitudes of j0, operating current
density, electrode thickness and also the symmetry of transfer co­
efficients. In the viewpoint of mathematics for a linear approximation of
exponent items exp(x) � 1 þ x, j ¼ j0(αþβ)neFη/RT is thought to be Fig. 1. Overpotential profile in a typical Ni/YSZ cermet anode of SOFC oper­
accurate at small overpotentials. Costamagna et al. [22] used such a ations at 1073.15 K. (a) Parameters are: electrode thickness l ¼ 10 5 m, ex­
linear approximation and made an analytical simulation model of change current density j0 ¼ 4.94 � 108 A m 3, operating current density I ¼
overpotential distribution in a SOFC electrode. In their simulation cases, 1000 A m 2, anodic and cathodic transfer coefficients α ¼ 1 and β ¼ 0.5. (b)
symmetric transfer coefficients α ¼ β ¼ 0.5 and ne ¼ 1 were used. For Parameters are: l ¼ 10 3 m, j0 ¼ 4.94 � 108 A m 3, I ¼ 4 � 104 A m 2, α ¼ 1
practical SOFC operation, multiple charge transfer steps are involved in and β ¼ 0.5.
and the global kinetics in the BV formulation generally leads to asym­
metric transfer coefficients [3,23]. As shown in Fig. 1a, in the case of α the anodic and cathodic transfer coefficients, etc. After analyzing the
¼ 1, β ¼ 0.5 and ne ¼ 2 and a thin electrode, there is an obvious relative shortcomings of the conventionally used linear approximation, we
error (although the absolute error is not big) between the fully numer­ introduce an approach of linear translation to improve the linear
ical and linearly approximated overpotential profiles, even if at a very approximation and use it to extend the scope of application of the
low operating current density or total overpotential. Bao et al. [24] power-law approach. The analytical profile of an activation over­
obtained an approximate analytical solution of the transport model in potential is therefore introduced for a multi-component mass transfer
electrodes for anode-supported SOFCs by using the perturbation method model to get an analytical concentration distribution and therefore a
on the basis of the linear approximation. Not surprisingly, the linear concentration overpotential. The algorithm is discussed in the context of
approximation makes a significant error in predicting the total over­ SOFCs, but applicable for the general electrochemical systems when the
potential, as shown in Fig. 1b. By keeping the nonlinear features of the BV kinetics is applicable.
system, Bao et al. [25] developed a hybrid algorithm of perturbation
method and power-law approach and realized an accurate predicting of 2. Mathematical model of charge transfer
the overpotential profile and the total overpotential. However, this
approach is currently just applicable for thick electrodes as it requires a A SOFC electrode (or the part of catalyst layer) is formed by a
minimum overpotential of zero (i.e. η ¼ 0 when dη/dz ¼ 0, z denotes the mixture of ionic (ion) conductor/electronic (el) conductor materials,
coordinate along the direction of the electrode thickness), which does and porosities are present in the structure for the gas diffusion, ion/
not hold for thin electrodes, as obviously shown in Fig. 1a. electron transport to make the three-phase boundary (TPB). At TPBs,
This paper is proposed to present an algorithm for accurate pre­ electrochemical reactions occur and current exchanges between ionic
dictions of an overpotential profile in SOFC electrodes, which is appli­ and electronic conducting phases. With the assumptions that all species
cable for a wide range of operations including cases of moderate or large are ideal gases, temperature and pressure are uniform throughout the
operating current densities (or total overpotentials), thin (electrolyte- electrode, both of the two conducting phases are continuous and ho­
supported) or thick (electrode-supported) electrodes, and asymmetry of mogeneous, the model here is steady-state and one-dimensional along z

2
D. Feng et al. Journal of Power Sources 449 (2020) 227499

coordinate, the direction of the electrode thickness. Only the Cartesian � � � � ��

αne F βne F
coordinate system (planar geometry) is discussed in the following j ¼ j0 exp ηact exp ηact
mathematics. RT RT
� � �� � �γi (5)
The charge transfer in the ionic and electronic conducting phases is Eact 1 1 pi;b
j0 ¼ j0;ref exp Π
governed by Ohm’s law, R T Tref p0

d2 φion d2 φ where the concentration correction factors (ci/ci,b) in Eq. (3) are
σion;eff ¼ σel;eff 2el ¼ j (1)
removed and j0 is relevant to the bulk species concentrations, which is
dz 2 dz
thought to be more consistent with the physical nature of exchange
where φ is the potential, and σk,eff ¼ σkfk (k ¼ el, ion) is the effective
current density (no concentration gradient at zero-overpotential or
ionic/electronic conductivity which should be corrected by the volume
equilibrium condition). In the next, we will firstly focus on the system of
fraction (fk) of the corresponding conducting phase in the electrode.
the activation overpotential and come up with a wide-range closed-form
The boundary conditions are: at the electrode/gas channel (E/C)
analysis. With the predetermined profile of an activation overpotential,
interface (z ¼ 0), the ionic current density completely transfers into the
we will then extend to get an analytical solution of multi-component
electronic phase, while at the electrode/electrolyte layer (E/E) interface
mass transfer and subsequently the concentration overpotential and
(z ¼ l), the electronic current density completely transfers into the ionic
the total overpotential in electrode, i.e. η ¼ ηact þ ηconc.
phase,
With the assumption of homogeneous ion/electron conducting
� �
dφ � dφ � phase, i.e. constant ionic and electronic conductivity at a given tem­
σion;eff ion �� ¼ σel;eff el �� ¼ 0
dz z¼0 dz z¼l perature, we can get the system of ηact as follows.
� � 8 2 � �
> d ηact 1 1 1 1
dφel �� dφion �� > ¼ j where ​ ​ ¼ þ
σel;eff ¼ σ ion;eff ¼I (2) >
< dz2 σt σt σel;eff σ ion;eff
dz �z¼0 dz �z¼l � � (6)
>
> dη � dη �
>
: σel;eff act �� ¼ σ ion;eff act �� ¼ I
where I is the operating current density, l is the electrode thickness. dz z¼0 dz z¼l
The electrochemical reaction rate j, can be described by the general
current-overpotential equation, Define the following variables,
�
�
creact
�
αne F
�
cprod
�
βne F
�� ζ ¼ z=l; y ¼ f ηact ; k0 ¼ j0 fl2 σ t (7)
j ¼ j0 exp η exp η (3)
creact;b RT cprod;b RT
where f ¼ neF/RT, we obtain the dimensionless system as follows.
where R is the universal gas constant, F is the Faraday constant, T is the � ’’
y ¼ k0 ½expðαyÞ expð� βyÞ�
operating temperature, ne is the number of electrons involved in the
� (8)
y’ ðζ ¼ 0Þ ¼ N0 ¼ Ifl σ el;eff ; ​ ​ ​ ​ ​ ​ y’ ðζ ¼ 1Þ ¼ N1 ¼ Ifl σ ion;eff
electrochemical reaction, ci,b and ci are bulk concentration and local
concentration at the reaction site of reactant or product species i, Because k0 is proportional to l2, the nonlinearity of the system is
respectively, α and β are the anodic and cathodic charge transfer co­ strongly related to the electrode thickness. In cases of thin electrodes, k0
efficients, β ¼ 1 – α holds for an elementary charge-transfer step but is small enough (k0 → 0), Eq. (8) approaches to be linear (y’’ � 0) and the
generally does not for an overall multi-step kinetics, j0 is the volume- difference between the absolute values of N0 and N1 is small, the over­
specific exchange current density, which is often semi-empirically potential presents a smooth distribution in the electrode, as shown in
taken a function of temperature and species concentrations following Fig. 1a. However, in mathematics, the linear system does not satisfy the
the framework of Arrhenius’s kinetics two Newman-type boundaries at the same time. In cases of thick elec­
� � �� � �γi trodes, k0 is big enough (k0 → ∞), by multiplying 1/k0 at both sides of
j0 ¼ j0;ref exp
Eact 1 1
Π
pi
(4) Eq. (8), it reduces from a differential system to an algebraic equation, of
R T Tref p0 which the algebraic solution does not satisfy the boundary conditions at
the two ends. Larger the value of l, larger the difference between the
where j0,ref is the reference exchange current density at the reference absolute values of N0 and N1. This denotes a strong boundary effect in
temperature Tref, Eact is the activation energy, pi ¼ ciRT and γ i are the physics, i.e. a sharp variation of the overpotential profile (or current
partial pressure and reaction order of species i, respectively, and p0 ¼ variation) close to the boundary, as shown in Fig. 1b. The thickness of
101325 Pa is the standard pressure. The fitting parameters of j0,ref, Eact, the concentration boundary layer (similar to the boundary layer in hy­
etc. can be fitted from the characteristics of Tafel slope rate, polarization drodynamics) has been discussed in Ref. [24].
resistance and temperature dependence of V–I and electrochemical For a linear approximation of exponent items, i.e. exp(αy) � 1 þ αy,
impedance spectra (EIS) measurements. exp( βy) � 1 – βy, an exact solution of Eq. (8) is obtained as follows.
The overpotential η is defined as the potential difference between the
two phases subtracts the potential difference in equilibrium, i.e. η ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
λ0
N1 N0 e N1 N0 eλ0
ylin ¼ e λ0 ζ þ e λ0 ζ
; λ0 ¼ k0 ðα þ βÞ (9)
φel φion ðφel φion Þeq . There are generally two ways in fuel cell λ0 ðeλ0 e λ0 Þ λ0 ðeλ0 e λ0 Þ
modeling to calculate the total overpotential [26]. The first method is The variable λ0 is similar to the Thiele modulus for a classical
based on the general current-overpotential equation with the effect of problem of diffusion and reaction of gases throughout a porous catalyst
species concentrations, i.e. adding concentration correction factor (ci/ci, pellet in isothermal conditions.
b) before the exponent terms in the BV equation, cf. Eq. (3). In this case, η
has been the summation of activation overpotential and concentration 3. Analytical solution of overpotential profile
overpotential. Because of coupled modeling of charge and mass transfer,
it is generally difficult to get an analytical solution in closed form, Let us first analyze the overpotential profiles in Fig. 1, which show
especially for multi-component fuel system. the advantages and disadvantages of the linear approximation and then
Another widely used method is to decouple charge and mass transfer guide us on how to improve.
by separately calculating activation overpotential (ηact) and concentra­
tion overpotential (ηconc), based on the standard the BV equation and (1) The linear approximation normally works well at low over­
linear model of Fick’s diffusion, respectively. For charge transfer only, potentials, except for some cases like that in Fig. 1a.

3
D. Feng et al. Journal of Power Sources 449 (2020) 227499

(2) Even if in the case of that in Fig. 1a, the linear approximation Using Y ¼ 0 when dY/dζ ¼ 0, Eq. (14) yields,
makes a good prediction of the shape of the overpotential profile, sffi�
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi��
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi�
ffiffiffi
i.e. the main error lies in a constant deviation from the exact 2N1 1þ
σion;eff λ σion;eff λ
σ el;eff e 1þ σel;eff e
solution. Such an error comes from the intrinsic structure of Eq. eαA e βA
¼ (16)
(8), which lacks Dirichlet-type boundary conditions for this λðeλ e λÞ αeαA þ βe βA

boundary-value problem.
(3) In the case of that in Fig. 1b, most of the electrode (the region out which is used to determine the magnitude of translation, A.
of the reaction boundary layer) is inactive to electrochemical
reactions. The linear approximation works therefore well in the 3.1.2. Power-law approach
inactive region, but naturally makes a higher error in predicting The idea of the power-law approach is to construct a power function
the overpotential at E/E interface, ηE/E, at the higher operating to approach the exponent items in BV equation in the region close to the
current density. On the other hand, the accuracy of the magni­ E/E interface, i.e.,
tude of ηE/E are more important than that of the overpotential
expðαyÞ expð βyÞ ¼ kðy þ aÞn (17)
profile in the active region (reaction boundary layer), as that the
active region is thin in a thick electrode and ηE/E is normally taken where k, a and n need to be determined.
as the total overpotential because σ ion ≪ σel holds generally in Substituted into Eq. (8), it yields,
fuel cell operations. � ’’ �
y ¼ k0 kðy þ aÞn t¼yþa t’’ ¼ k0 kt n
�
�! (18)
For a summary, the linear approximation can predict well the com­ y’ ð1Þ ¼ N1 t’ ð1Þ ¼ N1
plete shape of the overpotential distribution, while a linear translation
Considering a power-law solution of Eq. (18), i.e. t ¼ t(1)[1 þ b
should be introduced to alleviate the direct-current deviation at the low
(1 ζ)]m, we obtain m ¼ 2/(1 n) and the following expression.
overpotentials. At moderate and high operating current densities, the
system nonlinearity must be kept in the active region, which can be b2 k0 ð1 nÞ2
¼ ½yð1Þ þ a�n 1 ; 2b½yð1Þ þ a� ¼ ð1 nÞN1 (19)
solved by the power-law approximation to the exponent items. k 2ð1 þ nÞ
To get a smooth approximation of the constructed power function to
3.1. Algorithm
the exponent function, we require both of zero- and first-order differ­
entials of the two sides of Eq. (17) at the E/E interface.
3.1.1. Linear translation
With the signs of the 2nd-order derivative and the 1st-order deriv­ f1 ¼ eαyð1Þ e βyð1Þ ¼ k½yð1Þ þ a�n
ative at the both boundaries (d2y > 0, dy(0) < 0, dy(1) > 0), the solution (20)
f2 ¼ αeαyð1Þ þ βe βyð1Þ ¼ nk½yð1Þ þ a�n 1

of Eq. (8) presents a concave profile. Mark the minimum of the dimen­
sionless overpotential, Then k, a, n, and b can be explicitly solved from Eqs. (19) and (20).

ymin ¼ A (10) f2 N 21 f1 N 21
n¼ 2
; ​ ​ ​ ​ ​ ​ ​ ​a¼ yð1Þ;
2k0 f 1 f2 N 1 2
2k0 f 21 f2 N 21
and define, �� �n � (21)
f1 N 21 2N1 2k0 f 21 f2 N 21
yðζÞ ¼ YðζÞ þ A (11) k ¼ f1 2 2
; ​ ​ ​ ​b¼ 2
2k0 f 1 f2 N 1 ð1 nÞf1 N 1

where Y is the linearly translated dimensionless overpotential, we can Note that k, a, n and b are all dependent on y(1) ¼ ηE/Ef, i.e. the total
get the following system. overpotential is used here as a parameter instead of a performance
� ’’ � � (outcome) variable in a conventional fuel cell modeling. Notice the
Y ¼ k0 eαA expðαYÞ e βA expð βYÞ implicit exact solution of Eq. (8), i.e.,
’ ’ (12)
Y ðζ ¼ 0Þ ¼ N0 ; ​ ​ ​ ​ ​ ​ Y ðζ ¼ 1Þ ¼ N1 � �
1 1
The minimum of Y is zero, i.e. Y ¼ 0 when dY/dζ ¼ 0. In mathe­ ðy’ Þ2 ¼ 2k0 eαy þ e βy þ C (22)
α β
matics, we can therefore safely approximate the exponent items, exp
(αY) � 1 þ αY, exp( βY) � 1 – βY, in the region close to the minimum of and Eq. (10), i.e.,
Y. In physics, the first-order differential of the overpotential is corre­
sponding to the electrochemical reaction (or ionic/electronic current y’ ¼ 0; y¼A (23)
densities), the minimum region (dY/dζ ¼ dy/dζ ¼ 0) represents an
inactive region which is safe for linear approximation. Compared to the we get,
conventional linear approximation, the idea of the linear translation is to �
1 αA 1
�
linearly approximate the exponent items exp(αY) and exp( βY) instead C¼ e þ e βA
(24)
α β
of exp(αy) and exp( βy). As a result, it yields,
� ’’ � � and
Y ¼ k0 αeαA þ βe βA Y þ k0 eαA e βA
’ ’ (13) � �
Y ðζ ¼ 0Þ ¼ N0 ; ​ ​ ​ ​ ​ ​ Y ðζ ¼ 1Þ ¼ N1 1 1 1 1
2k0 eαyð1Þ þ e βyð1Þ
eαA e βA
¼ N 21 (25)
α β α β
which has an exact solution,
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi When α ¼ β, there is an exact solution as follows.
eαA e βA
YðζÞ ¼ C3 eλζ þ C4 e λζ
α
; ​ ​ ​ ​ ​ ​ ​ λ ¼ k0 ðαeαA þ βe βA Þ 0 2 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 31
αe þ βe
A βA � 2 �2
1 @1 4βN 21 βN 1
(14) yð1Þ ¼ ln þ eαA þ e βA þ þ eαA þ e βA 4 5A (26)
β 2 2k0 2k0
where,
Because Eqs. (25) and (26) come from the exact solution of Eq. (8), y
N1 N0 e λ N1 N0 eλ (1) will be significantly accurate as long as the linear translation, cf. Eq.
C3 ¼ λ ; C4 ¼ λ (15)
λðe e λÞ λðe e λÞ (23) holds well. Then the power-law solution in the region close to the

4
D. Feng et al. Journal of Power Sources 449 (2020) 227499

E/E interface can be obtained.

yp1 ¼ ½yð1Þ þ a�½1 þ bð1 ζÞ�2=ð1 nÞ

a (27)

Up to now, we got the linear approximation and power-law

approximation which are valid in the inactive region and active (close
to the E/E interface), respectively. Based on the analysis of boundary-
layer effect [24], we construct a hybrid function to make a smooth
transition between them to denote the overpotential profile in the whole
electrode,
�
y ¼ ylt þ yp1 ylt e λð1 ζÞ þ g1 ðζÞ (28)

where ylt(ζ) ¼ Y(ζ) þ A, cf. Eq. (11), denotes the above-mentioned lin­
early translated solution.
From the following requirements,
� �
g1 ð0Þ ¼ 0; g’1 ð0Þ ¼ 0; g1 ð1Þ ¼ 0; g’1 ð1Þ ¼ λ yp1 ð1Þ ylt ð1Þ (29)

the complementary function of g1 can be constructed as follows.

� �
g1 ðζÞ ¼ λ yp1 ð1Þ ylt ð1Þ ζð1 ζÞe 1:3λð1 ζÞ (30)

The power-law approximation can be obtained similarly in the re­

gion close to the E/C boundary, whatever it is normally unnecessary
because of the negligible boundary effect at the E/C interface as σion ≪
σ el in practical SOFC operation.

3.2. Validation of the algorithm

The algorithm is first tested for a thin Ni/YSZ (Yttria-stabilized zir­

conia) cermet anode in an electrolyte-supported SOFC. The anode
thickness is l ¼ 10 5 m. Fig. 2a and b shows the cases at the moderate
and high operating current densities, respectively. The numerical solu­
tion of Eq. (8) is solved by using the bvp4c function in MATLAB toolbox
for boundary value problems for ordinary differential equations. In both
cases, the conventional linear approximation (CLA) makes significant
overestimations of the overpotential, while the linear translation (LT)
solution works well singly and the hybrid solution of the linear trans­
lation and power-law (PL) approach further improves the performance.
Because of the small electrode thickness, even the power function at the
E/E interface works not badly in the whole region of the electrode. Fig. 2. Validation of the overpotential profile in a thin Ni/YSZ anode at the
operating current densities (a) I ¼ 1 � 104 A m 2 (b) I ¼ 4 � 104 A m 2. Pa­
Fig. 3 shows the cases at high operating current densities, respec­
rameters are: l ¼ 10 5 m, j0 ¼ 4.94 � 108 A m 3, α ¼ 1, β ¼ 0.5, ne ¼ 2, T ¼
tively, for a thick Ni/YSZ cermet anode in an anode-supported SOFC.
1073.15 K, σion ¼ 3.34 � 104 exp( 10300/T) S m 1, σ el ¼ 9.5 � 107/T.exp
The anode thickness is l ¼ 10 3 m. Because of the large electrode ( 1150/T) S m 1, f ¼ 21.6286, k0 ¼ 0.4714.
thickness, the electrochemical reaction (or the ionic current density)
occurs intensively in the region close to the E/E interface (i.e. the re­
4. Extensions of the approach
action boundary layer) and most of the anode is inactive. As a result, the
conventional linear approximation (CLA) works well in the inactive
4.1. Analytical solution to mass transfer
region (predicting almost zero overpotentials), but not surprisingly
overestimates significantly in the active region. In these cases, the linear
The mass transfer of multi-component gas species in a porous elec­
translation (LT) solution is coincident with the CLA, as the linear
trode is governed by the Stefan-Maxwell model (SMM), or the Dusty gas
translation magnitude A ¼ 0. The power function achieves a good
model (DGM). The DGM includes the Knudsen diffusion term in
approximation to the BV-type kinetics in the reaction boundary layer.
Maxwell-Stefan framework, as a result, the total pressure drop holds
Again, the hybrid solution (LT þ PL) works excellent in the whole region
when the effects of Knudsen diffusion become significant. In our previ­
of the electrode. Fig. 4 further shows the index of the power-law func­
ous work [27], we have discussed detailly the pressure gradient under
tion, n, with respect to the variation of the electric load. At the small
the framework of the DGM. Although there is pressure gradient when
operating current densities, a linear approximation (n ¼ 1) works well,
accounting for the effects of Knudsen diffusion, the pressure-driven
and in our expectation, the value of n increases sharply at the high
viscous flow contributes the convective flux and has a little influence
operating current densities.
on the profile of gas concentrations (which are dominant by diffusion).
Table 1 compared the above-calculated values of the minimum
With this regard, for an easy close-form analysis, we used the extended
dimensionless overpotential (ymin) and the dimensionless overpotential
SMM (eSMM) approach, which follows the framework of the SMM and
at the E/E interface (yE/E) between the numerical solution, conventional
takes the Knudsen diffusion effect in the effective Maxwell-Stefan dif­
linear approximation (CLA), linearly translated (LT) solution, and
fusivities [27]. The eSMM model has been validated in predicting the
hybrid solution combined with LT and power-law (PL) approach. The
concentration loss [27].
values of the linear translation magnitude (A) and the power index (n)
For 1D modeling, it is,
are also listed.

5
D. Feng et al. Journal of Power Sources 449 (2020) 227499

dxi Xxi Nj xj Ni
¼ (31)
dz j6¼i
cDij;eff

where c is the concentration of the gas mixture, xi ¼ ci/c and Ni are the
molar fraction and flux of species i, respectively, Dij,eff is the effective
diffusivity, which includes the binary Maxwell-Stefan diffusivity be­
tween species i and j (Dij) and their Knudsen diffusion coefficients (Di,K,
Dj,K) and is corrected by the electrode porosity/tortuosity ratio (ε/τ) as
[27],
� � rffiffiffiffiffiffiffiffi
ε 1 1 2 8RT
Dij;eff ¼ � � þ � � ; Di;K ¼ rP (32)
2τ 1 Dij þ 1 Di;K 1 Dij þ 1 Dj;K 3 πMi

where rP is the average particle size of electrode and Mi is the molecular

weight of species i.
The mass balance of gas species is,
dNi X
¼ υi j=ne F þ υi;k rk (33)
dz k

Fig. 3. Validation of the overpotential profile in a thick Ni/YSZ anode at the

where υi, υi,k are the stoichiometric coefficients of species i in the elec­
operating current densities I ¼ 4 � 104 A m 2. Parameters are identical to those
in Fig. 2 except for l ¼ 10 3 m, k0 ¼ 4713.7. trochemical reaction and the kth chemical reaction, e.g. reforming re­
action and water gas shift reaction.
The species concentration at E/C interface is the bulk concentration,
and species can not penetrate the dense electrolyte layer i.e. the
boundary condition are: xi(z ¼ 0) ¼ xi,b, Ni(z ¼ l) ¼ 0.

4.1.1. H2(1) H2O(2) inert(3) or CO(1) CO2(2) inert(3) system

This is typical for SOFC anodes using single fuel (H2 or CO).
Considering the inert species (N3 ¼ 0) and the equimolar counter
diffusion between H2–H2O or CO–CO2 (N1 ¼ N2) and no chemical
reactions (rk ¼ 0), substituting it into Eqs. (31) and (33), one gets
d2(lnx3)/dz2 ¼ υ1j(1/D13,eff 1/D23,eff)/cneF. Combined with d2ηact/dz2
¼ j/σt (cf. Eq. (6)), it yields d2(lnx3)/dz2 ¼ υ1σt/cneF(1/D13,eff 1/D23,eff)
∙d2ηact/dz2. Integrating it and considering the boundary conditions, i.e.
xi(z ¼ 0) ¼ xi,b, dηact/dz(z ¼ l) ¼ I/σ ion,eff and dxi/dz(z ¼ l) ¼ 0, we can
get the concentration profiles.
� � �� ��
σ t υ1 1 1 I
x3 ðzÞ ¼ x3;b exp ηact ðzÞ ηact ð0Þ z (34)
cne F D13;eff D23;eff σion;eff

D13;eff D23;eff x
x1 ðzÞ ¼ x1;b þ x3;b x3 ðzÞ þ � ln 3;b (35)
D12;eff D23;eff D13;eff x3 ðzÞ

Fig. 4. Power index (n) versus the variable operating current density. Param­
For the binary species, i.e. H2(1) H2O(2) or CO(1) CO2(2) system,
eters are identical to those in Fig. 3. x3 ¼ 0, it becomes
� �
υ1 σ t I
x1 ðzÞ ¼ x1;b ηact ðzÞ ηact ð0Þ z (36)
cD12;eff ne F σion
Table 1
Some results in Fig. 2a (case 1), Fig. 2b (case 2) and Fig. 3 (case 3). Note that, the profile of the activation overpotential ηact(z) and ηact(z
Numerical CLA LT LT þ PL ¼ 0) has been predetermined by using the algorithm in Section 3. The
Case 1 ymin 0.8276 1.2027 0.8395 0.8226
profile of the concentration overpotential is then calculated by the linear
yE/E 1.2551 1.6535 1.2724 1.2615 model of Fick’s diffusion, i.e.
REa 31.75% 1.38% 0.52%
A 0.8395 0.8395 RT x1;b x2 ðzÞ
ηconc ðzÞ ¼ ln (37)
n 0.6885 ne F x1 ðzÞx2;b
Case 2 ymin 1.6259 4.8109 1.7405 1.6838
yE/E 3.0448 6.6140 3.3019 3.0707
4.1.2. O2(1) N2(2) system
RE 117.23% 8.44% 0.85%
A 1.7405 1.7405 It is typical operation in SOFC cathodes. With substitution of N2 ¼
n 2.7721 0 into Eq. (31), it yields d2(lnx2)/dz2 ¼ υ1j/cD12,effneF and then the ni­
Case 3 ymin 0 0 0 0 trogen concentration profile.
yE/E 2.8888 4.5387 4.5387 2.8892 � � ��
RE 57.12% 57.12% 0.01% υ1 σ t I
A 0 0 x2 ðzÞ ¼ x2;b exp ηact ðzÞ ηact ð0Þ z (38)
cD12;eff ne F σ ion
n 7.941
a
RE denotes the relative error (RE) of the approximated yE/E referred to the 4.1.3. H2(1) H2O(2)–CO(3) CO2(4) inert(5) system
numerical one. This is general in SOFC anodes using syngas fuel, reformatted natural
gas and hydrocarbon fuel. The hybrid-fuel operation brings two main

6
D. Feng et al. Journal of Power Sources 449 (2020) 227499

different phenomena: the first is the chemical reactions (here is the

water gas shift reaction (WGSR), i.e. CO þ H2O → H2 þ CO2), and the
second is the electrochemical co-oxidation of H2 and CO, which can be The above expression includes the effects of the temperature and
modeled as, species partial pressures, where j0,ref is the reference exchange current
� �� �αn F � � �� density at the reference temperature of Tref, p and p0 are the total gas
βne F
j ¼ ϑH2 j0;H2 þ ΓϑCO j0;CO exp
e
ηact exp ηact (39) pressure and standard pressure, respectively, γi is the reaction order of
RT RT
species i, Eact is the activation energy, and ωA/E is the hydrogen current
fraction at the A/E interface. As Γ is a function of the operating current
where ϑH2 ¼ x1,b/(x1,b þ x3,b) and ϑCO ¼ 1 ϑH2 is the weighting factor
density and fuel compositions, Γj0,CO can be taken as the effective ex­
accounting for the H2/CO bulk compositions, Γ is defined as the
change current density of CO, the above algorithm establishes an
enhancement factor to denote the enhancement effect of the WGSR in
adaptive superposition mechanism of single H2-fuel and CO-fuel elec­
the hybrid-fuel atmosphere on the CO electrochemical rate compared to
trochemistry for hybrid-fuel operation.
that in the single CO-fuel atmosphere [3,28–30]. Instead of a constant, Γ
Fig. 5a compares the simulated V–I performance of a syngas-feed
here varies with respect to fuel compositions and current density. As Γ ¼
anode-supported SOFC button cell with the experiments operated by
1, the above model correlating hybrid-fuel kinetics with nonlinear su­
Jiang and Virkar [31]. The algorithm works well in the region of small
perposition of single-fuel rates, recovers the state-of-the-art line­
and medium operating current densities. The discrepancy at high cur­
ar-superposition rate [2].
rent densities and fuel utilization, i.e. in the region of concentration
To get an analytical solution to the five-species mass transfer, we
overpotential of V–I curves, could be alleviated with a more detailed
used here the mixture-average model (i.e. Fickian model) to govern the
mass-transfer model, not discussed here, by correcting the species TPB
mass transfer,
concentrations with consideration of surface diffusion [29,30]. To our
dxi X5 knowledge, there is generally a significant degradation of single-fuel V–I
Ni ¼ cDi þ xi Nj ði ¼ 1; 2; 3; 4Þ (40) performance as the inlet fuel molar fraction decreases, and the electro­
dz j¼1
chemical oxidation rate of CO in Ni-YSZ cermet is several times slower
where Di is the Fick diffusivity of species i, which can be estimated by the than that of H2 [32]. The hybrid-fuel operations show somehow
Wilke formulation and taken as constant. different behaviors with comparable V–I performance between 86%
The species mass balance in this case is, H2 14% CO and 20% H2 80% CO fuel system. The outcome can be
explained by the enhancement effect of CO kinetics within the
dN1
¼
dN2
¼
ωj
þ rWGS ;
dN3
¼
dN4
¼
ð1 ωÞj
rWGS (41) syngas-fuel atmosphere. As shown in Fig. 5b, the value of Γ is always
dz dz ne F dz dz ne F higher than unity, that is why we call it as ‘enhancement factor’, and can
be up to order of 102–103 at low CO concentrations. It is interesting that
where rWGS is the reaction rate of the WGS, ω is defined as the hydrogen the value of Γ � 4 in the case of 20% H2 80% CO is consistent with the
current fraction which denotes the contribution percent of the hydrogen ratio of electrochemical oxidation rates of H2 to CO in Ni-YSZ anode
electrochemical rate in the total. [32]. Besides to the comparable performance between CO-rich and
Accounting for the boundary conditions at the anode/electrolyte (A/ H2-rich fuel system, the modeling framework of the enhancement effect
P
E) interface (Ni|z¼l ¼ 0), which leads to N1 ¼ N2, N3 ¼ N4, Ni ¼ 0, also provides a perspective to explain the experiments by O’Brien and
2 2 2 2
dN1/dz þ dN3/dz ¼ cD1,effd x1/dz cD3,effd x3/dz ¼ j/neF, and Giorgi [33], in which they found 5–8 times exchange current density of
assuming the WGSR in equilibrium, we can get the H2 concentration 25% CO 75% H2 higher than that of pure H2 in the Ni0.7Co0.3/YSZ
profile [3,28,29], cermet anode. As shown in Fig. 5a, without the numerical compensation
bðzÞ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
b2 ðzÞ aGðzÞ of Γ > 1, the linear-superposition rate (Γ ¼ 1) makes an obvious un­
x1 ðzÞ ¼ (42) derestimation of V–I performance in cases of high CO concentration.
2D1;eff ðk 1Þ
Fig. 5c further shows the profiles of hydrogen current fraction (ω) at
where different H2/CO compositions within the active zone in the anode,
� � which is defined as a region where 99% of electrochemical reaction
GðzÞ ¼
σt
ηact ðzÞ ηact ð0Þ
Iz
þ D1;eff x1;b þ D3;eff x3;b (43) happens similar to the concept of the boundary layer in fluid dynamics
cne F σ ion [24].
�
k ¼ Keq;WGS D1;eff D4;eff D2;eff D3;eff �
a ¼ 4kðk 1Þ D1;eff x1;b þ D2;eff x2;b (44) 4.2. Thermo-electrochemical cell-level computation
�
bðzÞ ¼ ðk 1ÞGðzÞ þ k D1;eff x1;b þ D2;eff x2;b þ D3;eff x3;b þ D4;eff x4;b
Based on the above-mentioned charge and mass transfer model in
in which Keq,WGS is the thermodynamic equilibrium constant of the electrode level, we further carried out a two-dimensional thermo-elec­
WGSR. trochemical cell-level computation by introducing the gas flow and heat
Then the profile of ω is obtained as follows. transfer along the channel. As the thickness of the positive-electrolyte-
� negative (PEN) is generally much less than the length of flow channel,
cD n F d2 x1 d2 ηact
ωðzÞ ¼ 1;eff e (45) the 2D model was reasonably simplified into a 1Dþ1D computation. The
σt dz2 dz2 assumptions of the thermal flow along the channel include: no along-
The enhancement factor can be further determined [3,28,29]. the-channel gas diffusion, negligible heat conduction of gas phase,
negligible gas pressure drop, neglecting variation of specific heat, heat

� � ��
γH γH O � ðEact;H2 Eact;CO Þ 1 1 � �γH þγH O γCO γCO
j0;ref;H2 xH22;b xH22O;b ϑH2 1 ωA=E R T Tref p 2 2 2
​Γ ¼ γCO γCO2
e (46)
j0;ref;CO xCO;b ω
xCO2 ;b A=E ð1 ϑH2 Þ p0

7
D. Feng et al. Journal of Power Sources 449 (2020) 227499

Fig. 5. (a) Comparison between experimental and simulated polarization

curves at different H2/CO compositions. The experimental data come from
Jiang and Virkar [31]. The main parameters are: total gas pressure in anode and
cathode pa ¼ pc ¼ 1atm, cell temperature T ¼ 1073.15 K, H2 and CO exchange
current density j0,H2 ¼ 6.87 � 108 A m 3 and j0,CO ¼ 0.8 � 108 A m 3, anodic
and cathodic transfer coefficients in anode α ¼ β ¼ 0.5, porosity-tortuosity ratio
in anode ε/τ ¼ 0.075, and the volume fraction of the ionic conducting phase in
the anode is 0.3. The overpotential in cathode (ηc) was calculated by I ¼ 6850.8
(pcxO2,b/p0)0.25 [xO2,E/E/xO2,b.exp(Fηc/RT) exp( Fηc/RT)], where xO2,b and
xO2,E/E ¼ 1 (1 xO2,b)exp(IRTlc/4FpcDO2,N2,eff) are oxygen molar fraction at the
cathode/flow channel and cathode/electrolyte layer interface, respectively, in
which lc is the cathode thickness. The open circuit voltage was determined by
experiments rather than the Nernst equation. (b) Value of the enhancement
factor (Γ) with respect to the variation of current density and fuel compositions.
(c) Distribution of hydrogen current fraction (ω) in the active zone in anode
when the operating current density I ¼ 2 � 104 A m 2.

conductivity and heat of reactions with respect to temperature,

neglecting radiation between solid elements, uniform electrode poten­
tial along the channel due to the high electronic conductivity of the
electrode, and only steady state is considered. H2–H2O is fuel in the
anode and air in the cathode. In consideration of equimolar counter-
diffusion of H2–H2O, the along-the-channel model is governed in mole
frame as follows.
Total mass balance of gas,
dðca ua Þ
¼0 (47)
dx

dðcc uc Þ IðxÞ
¼ Sm;c ξrib (48)
dx 2ne F
Species mass balance of gas,
dxi υi IðxÞ
ca ua ¼ Sm;a ξrib ði ¼ H2 ; H2 OÞ (49)
dx ne F

dxO2 IðxÞ
cc uc ¼ Sm;c ξrib ð1 xO2 Þ (50)
dx 2ne F
Gaseous energy balance,
dTk
ck uk Cp;k ¼ Sh;k hk ðT Tk Þ ðk ¼ a; cÞ (51)
dx
Solid energy balance,
� �
d2 T X ΔHref
0 ¼ κs þ γ Sh;k hk ðTk TÞ γ a Sm;a ξrib IðxÞ þ Vcell (52)
dx2 k¼a;c k ne F

In the above equations, ck, uk, Tk are the total molar concentration,
molar-average velocity and temperature of gas in the anode (k ¼ a) and
the cathode (k ¼ c) side, T is the solid temperature which is used as the
local operating temperature in the isothermal electrode-level model
(along z coordinate), κs the heat conductivity of the solid phase, ΔHref is
the heat of reaction of H2 þ 1/2O2 → H2O at the reference temperature
(for simplicity, the temperature dependence of heat of reaction is not
considered), I(x) is the local current density along the x coordinate
(direction of channel length).
The geometry relevant parameters include: rh the hydrodynamic
radius of the channel, Sm and Sh the specific mass or heat transfer area in
(caption on next column) channel, respectively, ξrib the rib factor, γ the ratio of the gas flow sec­
tion area to the solid section area. For a planar SOFC (PSOFC)
configuration, there is,
Dch;k Wch 1 1 Wrib
rh;k ¼ ; Sm;k ¼ ; Sh;k ¼ ; ξrib ¼ 1 þ ;
2ðDch;k þ Wch Þ Dch;k rh;k Wch (53)
�
γ k ¼ Dch;k Wch ðAcon;a þ Acon;c þ lPEN Wch ζrib Þ

where Dch and Wch the depth and width of the channel, Wrib the width of
the rib, lPEN ¼ la þ le þ lc the thickness of the PEN, in which la, le, lc are

8
D. Feng et al. Journal of Power Sources 449 (2020) 227499

the thickness of the anode, electrolyte layer and cathode, respectively,

Acon the section area of the interconnector, i.e. Acon,k ¼ Wchξribδcon
WchξribDch,k, in which δcon is the interconnector thickness. The convec­
tive heat transfer coefficient, h is estimated from the corresponding
Nusselt (Nu) number which is relevant to the channel depth/width ratio,

Nu ¼ 11:1197r4 31:4600r3 þ 35:1816r2 19:3996r þ 7:5385

(54)
r ¼ minðDch =Wch ; Wch =Dch Þ
�
hk ¼ λk Nu 4rh;k (55)

where λ is the heat conductivity of gas.

Fig. 6 shows the simulation results of the along-the-channel distri­
butions of gaseous/solid temperature, gas concentration, current density
and overpotential, in case of concurrent flow of fuel and air. The
computational variables are: la ¼ 5.4 � 10 4 m, le ¼ 1 � 10 5 m, lc ¼ 5.5
� 10 4 m, channel length L ¼ 0.088 m, Dch,a ¼ Dch,c ¼ Wch ¼ Wrib ¼ 2 �
10 3 m, δcon ¼ 2.5 � 10 3 m, inlet fuel (97% H2 3% H2O) temperature
Ta,in ¼ 973.15 K, inlet air temperature Tc,in ¼ 973.15 K, equivalent fuel
and air current density (i.e. the current density when the utilization
factor of fuel or air is unity) Ieq,a ¼ 2 � 104 A m 2, Ieq,c ¼ 4 � 104 A m 2,
and σ ion ¼ 3.34 � 104.e 10300/T S m 1, σ el ¼ 9.5 � 107/T.e 1150/T S m 1,
ε/τ ¼ 0.32/8.586, α ¼ β ¼ 0.5, γH2 ¼ 0.5, γH2O ¼ 0, j0,ref ¼ 1.0 � 108 A
m 3, Eact ¼ 130 kJ mol 1 in the anode and σ el ¼ 4.2 � 107/T.e 1150/T S

Fig. 7. Along-the-channel profile of (a) fuel, air and solid temperature, (b) local
current density, (c) overpotential and open circuit voltage, and (d) hydrogen
and oxygen molar fraction, from a counterflow computation. Solid lines
represent analytical approximation solution and discrete points with markers
are fully numerical simulation results.

m 1, ε/τ ¼ 0.4/2, α ¼ β ¼ 0.5, γO2 ¼ 0.465, j0,ref ¼ 8.385 � 108 A m 3,

and Eact ¼ 125 kJ mol 1 in the cathode, reference temperature Tref ¼
1073.15 K. Fig. 7 shows the simulated results in case of countercurrent
flow of fuel and air. The flow configuration, i.e. coflow or counterflow,
makes a significantly different temperature distribution, and also the
along-the-channel profile of the current density duet to the tradeoff
between local temperature and fuel concentration. All the simulation
show good performance of the analytical algorithm, meanwhile the
introduction of the analytical electrode-level model in closed form
contributes much higher computational efficiency than the fully nu­
merical solution for the 1Dþ1D computation.

5. Conclusions

Accurate and fast prediction of overpotential and species concen­

trations in the porous electrode is valuable for system-level optimization
and controller design in fuel cell science and technology. The conven­
tional linear approximation (CLA) of the exponent items in the Butler-
Volmer equation leads to significant discrepancies in the estimation of
the overpotential profile at the moderate and high electric loads. On the
other hand, due to keeping the intrinsic system nonlinearities, the
power-law (PL) approach can make good predictions at the high electric
Fig. 6. Along-the-channel profile of (a) fuel, air and solid temperature, (b) local loads, but it is still limited to cases of thick electrodes [25].
current density, (c) overpential and open circuit voltage, and (d) hydrogen Following the advantages of CLA, i.e. simple framework and good
molar fraction, from a coflow computation. Solid lines represent analytical predictions of the overall shape of overpotential profiles, this paper has
approximation solution and discrete points with markers are fully numerical presented an algorithm of linear translation (LT), cf. Eq. (16), based on a
simulation results. 1D dimensionless charge-transfer model. The LT approach significantly

9
D. Feng et al. Journal of Power Sources 449 (2020) 227499

improved the performance of linear approximations at the moderate and interests or personal relationships that could have appeared to influence
high operating current densities for thin electrodes. Furthermore, the LT the work reported in this paper.
approach directly found the magnitude of the minimum overpotential,
cf. Eq. (23), consolidated the mathematical base of the PL approach and Acknowledgments
therefore expanded its scope of application. As a result, the hybrid so­
lution of the LT þ PL approach, cf. Eq. (28), was convincingly validated This work was supported by the Natural Science Foundation of China
in a wide range of SOFC operations, including moderate and high [grant number 51976009 and 51606008] and the Beijing Science and
operating current densities, electrolyte- and electrode-supported elec­ Technology Project [grant number Z181100004518004].
trodes, and asymmetric anodic and cathodic transfer coefficients.
Interestingly, in the LT þ PL approach, the overpotential at the elec­ References
trode/electrolyte (E/E) interface, ηE/E was used as a parameter instead of
a performance variable in conventional SOFC models, cf. Eqs. (25) and [1] A. Kulikovsky, Analytical Modelling of Fuel Cells, Elsevier, Amsterdam, 2010.
[2] T.A. Prokop, K. Berent, H. Iwai, J.S. Szmyd, G. Brus, Int. J. Hydrogen Energy 43
(26). As obtained from the implicit exact solution of the original charge- (2018) 10016–10030.
transfer differential system, the value of ηE/E is quite accurate and can be [3] C. Bao, Y. Wang, D.L. Feng, Z.Y. Jiang, X.X. Zhang, Prog. Energy Combust. Sci. 66
directly used as the overall overpotential when the overpotential profile (2018) 83–140.
[4] C. Bao, Analytical electrochemistry and transport in solid oxide fuel cells – from
is not needed. The analytical overpotential profile was then applied to simulation to modeling. https://doi.org/10.13140/RG.2.2.25394.71366, 2018.
obtain the analytical species concentration profiles and concentration [5] A.J. Bard, L.R. Faulkner, Electrochemical Methods: Fundamentals and
overpotential based on a decoupled charge and mass transfer model. Applications, second ed., Wiley, New York, 2001.
[6] R.J. Gorte, J.M. Vohs, Ann. Rev. Chem. Biomol. Eng. 2 (2) (2011) 9–30.
Everything in the electrode-level model is in closed form, when applied [7] S.B. Adler, Solid State Ion. 135 (2000) 603–612.
in thermo-electrochemical cell-level computation, it contributes much [8] C.W. Tanner, K.Z. Fung, A.V. Virkar, J. Electrochem. Soc. 144 (1997) 21–30.
higher numerical efficiency than the fully numerical computation. [9] J.D. Nicholas, S.A. Barnett, J. Electrochem. Soc. 157 (2010) B536–B541.
[10] J.D. Nicholas, L. Wang, A.V. Call, S.A. Barnett, Phys. Chem. Chem. Phys. 14 (2012)
Note that, the local concentration in the exchange current density
15379–15392.
model is a common practice in the SOFC electrode simulation, corre­ [11] S.D. Ebbesen, S.H. Jensen, A. Hauch, M.B. Mogensen, Chem. Rev. 114 (2014)
spondingly a combination of activation and concentration overpotential 10697–10734.
is generally applied in the general current-overpotential equation. [12] A. Leonide, S. Hansmann, A. Weber, E. Ivers-Tiffee, J. Power Sources 196 (2011)
7343–7346.
However, how to correlate the local concentration (e.g. the species re­ [13] A. Kromp, A. Leonide, A. Weber, E. Ivers-Tiffee, J. Electrochem. Soc. 158 (2011)
action orders) is unclear due to the complexity of the multi-step mech­ B980–B986.
anism in SOFCs. Aiming at a performance-oriented modeling, we [14] A. Kromp, S. Dierickx, A. Leonide, A. Weber, E. Ivers-Tiffee, J. Electrochem. Soc.
159 (2012) B597–B601.
calculated the activation and concentration overpotential separately, i. [15] H.Y. Zhu, A. Kromp, A. Leonide, E. Ivers-Tiffee, O. Deutschmann, R.J. Kee,
e. firstly got the profile of the activation polarization (correlating the J. Electrochem. Soc. 159 (2012) F255–F266.
exchange current density with the bulk concentrations) and then applied [16] S. Kakac, A. Pramuanjaroenkij, X.Y. Zhou, Int. J. Hydrogen Energy 32 (2007)
761–786.
it in the mass-transfer model (the diffusion effect) to get the species [17] S.A. Hajimolana, M.A. Hussain, W.M.A.W. Daud, M. Soroush, A. Shamiri, Renew.
concentration profiles and subsequently the concentration loss. The Sustain. Energy Rev. 15 (2011) 1893–1917.
applicability and accuracy of our approach was validated in predicting [18] K.N. Grew, W.K.S. Chiu, J. Power Sources 199 (2012) 1–13.
[19] M. Andersson, J.L. Yuan, B. Sunden, Appl. Energy 87 (2010) 1461–1476.
the performance of the button cell experiments. The advantages of the [20] H.Y. Zhu, R.J. Kee, V.M. Janardhanan, O. Deutschmann, D.G. Goodwin,
two-step method will be more obvious when a more complex mass J. Electrochem. Soc. 152 (2005) A2427–A2440.
transfer model is used and a closed-form analytical solution is [21] S.H. Chan, K.A. Khor, Z.T. Xia, J. Power Sources 93 (2001) 130–140.
[22] P. Costamagna, P. Costa, V. Antonucci, Electrochim. Acta 43 (1998) 375–394.
anticipated.
[23] S. Nagata, A. Momma, T. Kato, Y. Kasuga, J. Power Sources 101 (2001) 60–71.
Although presented in the context of SOFCs, in which the BV equa­ [24] C. Bao, N.S. Cai, AIChE J. 53 (2007) 2968–2979.
tion generally provides a well-practical framework for performance- [25] C. Bao, W.G. Bessler, J. Power Sources 210 (2012) 67–80.
oriented analysis although its applicability in science is still in argu­ [26] C. Bao, Y.X. Shi, E. Croiset, C. Li, N.S. Cai, J. Power Sources 195 (2010)
4871–4892.
ments for solid-oxide systems, the algorithm in this work is generally [27] C. Bao, Z.Y. Jiang, X.X. Zhang, J. Power Sources 310 (2016) 32–40.
applicable for other kinds of electrochemical systems when the Butler- [28] C. Bao, X.X. Zhang, Electrochim. Acta 134 (2014) 426–434.
Volmer kinetics holds. [29] C. Bao, Z.Y. Jiang, X.X. Zhang, J. Power Sources 294 (2015) 317–332.
[30] C. Bao, Z.Y. Jiang, X.X. Zhang, J. Power Sources 324 (2016) 261–271.
[31] Y. Jiang, A.V. Virkar, J. Electrochem. Soc. 150 (2003) A942–A951.
Declaration of competing interest [32] Y. Matsuzaki, I. Yasuda, J. Electrochem. Soc. 147 (2000) 1630–1635.
[33] J.S. O’Brien, J.B. Giorgi, J. Power Sources 200 (2012) 14–20.
The authors declare that they have no known competing financial

10