1

Simulation study about geometrical shaping of electrode-electrolyte contact

F. Delloro*, M. Viviani

*Corresponding author. E-mail: f.delloro@ge.ieni.cnr.it Phone: +39 3477862093
Institute for Energetics and Interphases, National Research Council,
via De Marini 4, 16149 Genova, Italy

Abstract
The main idea investigated in the paper is that a patterned contact
between composite electrode and electrolyte can improve
performances of the assembly. Shaped contact has been explored
by means of 2D finite elements simulations. Results show the
conditions in which the idea can be profitably utilized, bearing
performance increment up to 45 % in terms of current density
drawn.

Keywords: solid oxide fuel cell, composite electrode, patterned
electrode-electrolyte contact, modeling, simulation

1. Introduction

2

Current research in Solid Oxide Fuel Cells (SOFC) systems is mainly focused on the
optimization of materials and fabrication techniques. In fact lowering cost and improving
durability and performance are required to reach commercial standards and bring SOFC
systems to industrial production scale. Particular attention is given to electrodes, which
have to exploit several mansions, i.e. reaction and transport of electrons, ions and gaseous
reactants and products. Therefore materials used for electrodes have to fulfill a number of
requirements: catalyze the reaction, optimize the volumetric density of active sites, assure
a sufficient gas phase transport through the pores, minimize ohmic losses due to ionic and
electronic transport. The difficulties in the optimization of such a material are clearly
understandable looking at number and diversity of these requirements.
Mixed ionic and electronic conductors (LSC [1], LSCF [2], BSCF [3] for the cathode,
Fe-doped Ceria [4], LSCM [5] for the anode), as well as mixtures of electronic and ionic
conductor materials have been widely studied and used for this purpose [6]. The main
advantage of using a composite electrode (mixed ionic and electronic conductor, porous
materials) is to increase greatly active triple phase boundary (TPB) length. In fact, for an
electrode made with electronic conductor only, TPB is located at the interface with the
electrolyte. For a composite electrode, instead, TPB is spread all over the volume. We
have to further distinguish between active and non-active contacts, the latter occurring
when at least one of the three phases belongs to a cluster not connected with the
respective sink or source. The importance of TPB length is readily explained looking at
the typical kinetics of an electrodic process, as the Buttler-Volmer expression (see for
example [7] or eq. 14 below). Exchange current density can be written as a function of
the active TPB length:

(1)
3

where
aTPB
[m/m
3
] is the active TPB length density and i
0
[A/m] is the exchange current
density per unit TPB length, depending on temperature, gas composition and material. A
problem arising in the development of mixed and porous materials is the low
conductivity, compared to the pure dense conductor. Especially when reaction is spread
in a consistent part of the volume, ohmic drop due to the electrode can heavily affect the
performance of the cell.
Recently also the use of functional layers at the interface with current collector and
electrolyte has been considered. To obtain a good current distribution at the
interconnect/electrode interface, a porous layer containing electrocatalyst material is
deposited and works as current collector/distributor [8]. On the same side, in order to
prevent electrode poisoning by elements normally contained in the interconnects (for
instance Cr), a very thin compact layer with possibly high electronic conductivity is
grown on the interconnect surface [9]. An additional thin functional layer can be included
in the electrode/electrolyte interface [10, 11] to improve electrode electrochemical
activity. This is often similar to the bulk electrode, but with a different grain size
distribution and morphology assuring a greater number of reaction sites.
In this paper we present an alternative strategy for enhancing the performance of the
cell. We investigated, by means of 2D finite-element simulations, the idea of designing a
structured contact between electrode and electrolyte, instead of the usual flat one. In fig.
1 some simulation domains are showed, where it is clearly visible the stem (this is how
we decided to call the electrolyte rectangle in the electrode area). It has to be considered
as a repeating unit, spreading all over the electrode-electrolyte contact area. In the
simulations we compared the performance of conventional electrode-electrolyte contact
with the structured one, in a wide range of conditions.
4

Having to fix some of the simulation parameters (e.g. conductivities, exchange
current density), we chose to test our idea in the cathode side of a hydrogen fed SOFC,
which is a main source of irreversibility in the cell. The cathode is composed of a
Lanthanum Strontium Manganite (LSM) and Yttria-stabilized Zirconia (YSZ) mixture,
the most conventional and stable couple of materials for application in SOFCs [1214].
Results obtained are of general interest for all kinds of reacting layer. For example, the
same technique can be applied to the central membrane of the IDEAL-cell (see [15]
and companion papers).

2. Mathematical model

2.1 Introduction
The overall fuel cell reaction (H
2
+ O
2
H
2
O) is separated into 3 parts, one
(oxygen reduction) taking place at the cathode and the others at the anode.

Cathode O
2
+ 2 e

O
2
(2)

Anode H
2
2 H
+
+ 2 e

2 H
+
+ O
2

H
2
O
We will focus on the cathodic electrochemical reaction only (2), having left the anode out
of our simulation domain. The model presented in this paper is developed to describe the
polarization behaviour of the cathode-electrolyte assembly of a conventional SOFC, with
the main aim of providing a tool for the comparison of the standard flat electrode-
electrolyte interface with the patterned one. In particular the model has to predict the
5

polarization losses as a function of operating conditions, cell design (stem height and
distance between adjacent stems), material properties (e.g. conductivities of ionic
species) and morphological properties (e.g. porosity, active TPB length, volume fraction
and tortuosity of percolative paths).
The electrolyte is a dense ion-conducting layer (YSZ). The cathode is a composite porous
layer, made of a mixture of ion-conducting particles (YSZ) and electron-conducting
particles (LSM). Here, electrochemical reduction of oxygen can only take place in the
proximity of a triple phase contact, where molecular oxygen, coming from porosity, and
electrons, coming from electronic conducting phase, can meet to form oxygen ions close
to the ionic conductor material, which constantly carries oxygen away from the reaction
sites towards the electrolyte.

2.2 Model assumptions
The model is based on the following assumptions:
1) Steady state conditions.
2) Temperature is uniform throughout the assembly, i.e. heat transport is neglected.
3) The porous electrodes are treated as continuum of the three involved phases. The
transport within each phase (gas species, electrons and oxygen ions, respectively) is
described using effective transport coefficients.
4) Electronic resistivity of LSM is neglected, being 2 orders of magnitude lower that
ionic resistivity of YSZ (see for example [16]).
5) No mixed electronic-ionic conduction in either the electron-conducting phase or the
ion-conducting phase.
6

6) No electronic current at the cathode-electrolyte interface. This assumption is
supported by the low number of active reaction sites at the interface compared to
those in the bulk volume of the cathode.
7) Electrochemical reaction kinetically limited by a single step, e.g. the transfer of a
single electron, with no other limitation (e.g. adsorption, dissociation, migration),
thus resulting in a standard B-V kinetic.
Notice that there are 3 variables involved in our PDE system: |
YSZ
(potential in ionic
conductive phase), c
O
and c
N
(oxygen and nitrogen partial pressures in the gas phase).

2.3 Transport equations
Recalling assumption 4), we neglect ohmic drop in the electron-conducting phase, so we
can consider the potential of this phase fixed to a reference value. This simplify our
calculations, since also the potential difference between LSM and YSZ,
LSM YSZ
| | | = A ,
depends only by |
YSZ
.
Continuity equation for oxygen ions can be written:


F i N
V O
= V


(3)
where
V
M
i is the current volumetric density corresponding to the volumetric reaction rate
of oxygen reduction, through eq. (2). Under stationary conditions, the molar flux of
oxygen ions is associated to current density in the conducting phase, which is related to
the potential difference through the Ohm law.

YSZ
eff
YSZ
O
2 2
|
o
V = =
F F
i
N

(4)
The effective conductivity is related to the dense material property as:
7



YSZ
eff
YSZ
o o q =
(5)
where we introduce the effectiveness factor q. Its effect was studied in a set of dedicated
simulations and results will be shown later.
Regarding the gas phase transport, we assume the cathode compartment fed by air and
consider cathodic atmosphere composed only by oxygen and nitrogen. We decided to
apply the Dusty Gas Model (DGM in the following) equations for 2 species to describe
mass transport in the gas phase. The model is presented in [17], but the revision of the
paper by the authors revealed an error, corrected in eq. 12 in the present paper. The DGM
visualizes the porous medium as composed by a fix structure of solid particles (dust). Gas
molecules move in the space between dust particles following the kinetic theory of gases.
For each species the total molar flux is given by the sum of diffusion and of permeation
fluxes. Total flux
d f
i i i
N N N + = appears in the continuity equation for species i:

i i
i
N S
t
c

V = +
c
c
c
(6)
where c is the porosity and S
i
are source/sink terms. Steady state diffusion is described by
the following set of equations:


) (
T
1
m
d d
k
d
i
n
i j
j ij
j i i j
i
i
y c
D
N y N y
D
N
V =

=
=

(7)
where the effective Knudsen diffusion coefficient of component i, D
i
k
, and the effective
bulk diffusion coefficients of pairs i-j, D
ij
m
, are defined as:

8


m m
2 / 1
k
8
3
2
ij ij
i
g
i
D
M
T R
r D D =
|
|
.
|

\
|
=
t


(8)
is the ratio tortuosity over porosity, r the mean pore radius and M
i
the molar mass of
species i,
m
ij
D is the binary diffusion coefficient, given in [18] as

1 1 1
2 3 / 1 3 / 1
75 . 1
m
) ( 2
) (
00143 . 0

+ =
+
=
j ij
j i ij
ij
M M M
V V pM
T
i
D

V
i
being the special Fuller et al. diffusion volume.
For the permeation flux we use the Darcy equation:

T
f
c
p By
N
i
i
V =

(9)
where the permeability B can be calculated as a function of the morphological properties
of the membrane by using empirical expressions. We used the one given by Blake-
Kozeny:

( )
2
3
2
p
1 72 c
c
t
=
d
B
(10)
is the gas mixture viscosity, evaluated as

1/2
N N
1/2
O O
1/2
N N N
1/2
O O O
M c M c
M c M c
+
+
=


Using eqs. (7) and (11) in eq. (6), we come to a linear expression:

= V
j
j ij i
N F c


(11)
Matrix elements are given in [17], with a little error: o
i
in place of o
j
in the first eq. (12).
The correct expression is:
9

+
=
+ + =
= + =
=
j
k
j j
k
i
i
k
i
i i
i j
m
ij
j
k
i
ii
k
i
j i
m
ij
i
ij
D c B
D B
D
c
D
c
c D
F
j i
D
c
D c
c
F
/ 1
/
1 1
T
T
o
o
o



(12)
Inversion of eq. (11) enables us to write continuity equations (6) for the gas components
as a function of concentrations only:

( )

V V = +
c
c

j
j ij i
i
c F S
t
c
1
c
(13)
Eqs. (13), together with eqs. (6) for the potential in the ion conducting phase, constitute a
set of PDEs which can be solved in the simulation domains, given opportune boundary
conditions. The equations are coupled by reaction source/sink terms S
i
, which in our case
can be written, due to mass and charge conservation: F i S 2 /
V O
=

,

0
N
= S .

2.4 Macro-kinetic model.
Assuming that in the oxygen reduction there is a single rate determining step and that it is
an electron transfer process, then it can be shown that reaction kinetic is of the Butler-
Volmer type:


( )
(

|
.
|

\
|

|
.
|

\
|
=
act M M
0
V V
1 exp exp q o q o
RT
F
RT
F
i i
act

(14)

where i
V
is the local volume-specific reaction rate, i
0
V
the volume-specific exchange
current density, which depends on the microstructure as shown in eq. (1). q
act
is the
activation overpotential defined such that a positive value leads to water production.
10

Symmetry factor o
M
is fixed to 0.5. In the simulations we explored a wide range of
variation of the kinetic parameter i
0
V
, from 10
5
to 10
10
A/m
3
.
The activation overpotential q
act
depends on the potential difference A| between LSM
and YSZ in the composite electrode. With the assumption that |
LSM
is a constant, q
act

becomes a function of the YSZ potential only:


( ) ( )
YSZ eq YSZ, LSM YSZ LSM YSZ act
| | | | | | | | q = = =
eq eq

(15)
Note that the equilibrium potential |
YSZ,eq
depends on local concentrations according to
the Nernst equation:


|
|
.
|

\
|
=

) (O
) O (
ln
2 2

2 2
gas 2,
eq YSZ,
a
a
F
RT
F
G
|

(16)

where a(...) means activity of the species specified in brackets.
By using eq. (16) to calculate A|
YSZ, eq
, concentration polarization is implicitly introduced
into the model.

2.5 Simulation domain and boundary conditions:

We solved the equations presented above in different domains. Imposing periodic
boundary conditions on the left/right edges, domains represent repeating units. In other
words, we obtain simulation domains cutting the cell over its symmetry axis. For the flat
contact, the basis area can be chosen randomly, because there exist infinite symmetry
planes. Boundary conditions are given in table 2.
11

The aim of our simulations was to compare current density from reference simulation, the
flat electrode-electrolyte contact (fig. 1, right), with the ones from structured contact
simulations (fig. 2, left).


Fig. 1: (left) Reference simulation domain: the flat contact. (right) stem geometry and
parameters to be varied. Horizontal size of the stem (b) was fixed to 20m in the
simulations, while d and h were varied in the range 10-40m.

Referring to fig. 1 boundary conditions are listed in table I.





Variable 1 2 3 4
|
YSZ
|
0
0 Periodic
12

p
O2
- 0.21 Periodic
p
N2
- 0.79 Periodic
Table I: Boundary conditions referring to fig. 1. |
0
is total overpotential and was fixed at
0.6 V in the simulations.

The values of parameters used in the simulations are listed in Table II, while a summary
of the constitutive equations of the model is given in Table III.

Value Units Reference
Constants
F, Faradays constant 96485.3 C mol
-1

R, gas constant 8.314472 JK
-1
mol
-1

o
M
, symmetry factor 0.5
Operating conditions
Temperature, T 1073 K
Outside pressure, P
out
1 atm
Outside oxygen partial pressure, P
O
0.21 atm
Cell design parameters
Thickness of electrode, o
ele
50 m
Thickness of electrolyte, o
ely
20 m
Stem base 20 m
Material properties
Resistivity YSZ,
YSZ
1.19 O m (7)
Porosity, c 0.5 (1)
Tortuosity, t 4 (1)
Mean pore radius, d
p
10
-6
m
Gas phase properties
Oxygen viscosity,
O
3.2572 10
-5
Kg m
-1
s
-1
(18)
Nitrogen viscosity,
N
3.8034 10
-5
Kg m
-1
s
-1
(18)
13

Special Fuller et al. diffusion volume, V
N
18.5 (18)
Special Fuller et al. diffusion volume, V
O
16.3 (18)
Table II: Model parameters.

Pysicochemical process Equation

Oxygen reduction kinetics
( )
(

|
.
|

\
|
|
.
|

\
|
=
act M M
0
V V
1 exp exp q o q o
RT
F
RT
F
i i
act

YSZ eq YSZ, act
| | q =

|
|
.
|

\
|
=

) (O
) O (
ln
2 2

2 2
gas 2,
eq YSZ,
a
a
F
RT
F
G
|

Porous electrode gas transport

oxygen
( )
n n o o o
o
c F c F
F
i
t
c
V + V V = +
c
c

1
o
1 V
2
c

nitrogen
( )
o o n n n n
n
c F c F
t
c
V + V V =
c
c

1 1
c
Potential distribution
V YSZ YSZ
2i q = V V | o

Table III: Summary of model equations.

3. Results and Discussion

Two sets of stationary simulations at fixed total overpotential have been performed. In
the first we fixed the geometry and varied the effectiveness factor q and kinetic parameter
i
0
V
. In the second we fixed q and varied geometry and i
0
V
.

3.1 Effect of the kinetic parameter
14



Fig 2: current density distribution at fixed overpotential with different i
0
V
(left 10
5
, right
10
8
A/m
3
)

The effect of the kinetic parameter is readily explained looking at fig. 2 showing
volumetric current density distributions. We remark that reaction rate differs from current
density by a factor (r
V
=i
V
/2F).
In the left picture, corresponding to a simulation with i
0
V
10
5
A/m
3
, reaction is spread all
over the volume available. Electrochemical reaction kinetic is the limiting factor, so we
will refer to similar situations as kinetic regimes. This configuration clearly represents
a non-optimized cathode, since a thicker one will show better performances, having
larger volume available for reaction.
In the right picture, corresponding to a simulation with i
0
V
10
8
A/m
3
, we notice that
reaction is limited within a narrow region, close to the electrolyte. We will refer to it as
the active region, the upper part of the electrode not hosting any reaction. Increasing
15

further the kinetic parameter, the active zone results narrower. Current density profiles at
different i
0
V
and fixed overpotential can be interpreted introducing a distance a, measured
from the electrode-electrolyte interface, at which current is almost zero. We will refer to
this distance as the active range. Simulation results of fig. 2 show that for low i
0
V
(left)
the active range a is bigger than electrode thickness o
ele
, while for higher values of i
0
V

(right) a is smaller than o
ele
. When i
0
V
approaches infinity, a tends to zero and the active
zone is ideally reduced to the surface of contact between electrode and electrolyte. In this
case activation overpotential will be negligible with respect to ohmic potential drop of the
electrolyte. If also concentration overpotential is small, well refer to similar situations as
ohmic regimes. We remark that the analysis made here on the active range is valid only
when all conditions and parameters (total overpotential, conductivities, gas transport
parameters, etc.) are kept fixed.

3.2 Variation of effectiveness factor
We decided to explore in our simulations a whole range of conditions, in particular
focusing on material properties. In this first set of simulations we studied the effect of
effectiveness factor and kinetic parameter. To this aim, we fixed the geometry choosing
the following parameters: d = 10m and h = 40m and ran the simulations changing i
0
V

on a logarithmic scale in the range 10
6
10
9
A/m
3
and q between 2 20. The interval
chosen for i
0
V
is representative of a variety of intermediate situations between activation
regime and ohmic regime. The values chosen for q start from 2, an optimistic improbable
case, and 20, value calculated by means of percolation theory when the two phases have
particles of the same size and porosity is 50%. Probably an optimized material can show
a q close to 10.
16

Results are presented in fig. 3 where the percentage of current variation with respect to
the reference case (100*(i - i
ref
)/i
ref
) is plotted as a function of q. Each curve corresponds
to a different i
0
V
. Starting from the lowest, 10
6
A/m
3
, performance of structured contact is
worse than the flat one, except when q > 18. Again here cathode thickness is not
optimized, because with increased reaction volume (i.e. removing the stem) current is
higher. The best performances for this geometry are obtained for i
0
V
in the range 10
7
-10
8

A/m
3
. When it reaches 10
9
A/m
3
, current is only slightly higher than the reference case
and it will fall below it for even higher i
0
V
(this simulation is not shown here). This trend
is clearly due to the transition to ohmic regimes: in the structured contact design, the
electrolyte has a mean thickness higher than in the reference case, so that, when
electrolyte loss becomes dominant, structured contacts show worse performances.
Resuming, in activation and ohmic regimes the geometry tested is not satisfying, while,
for intermediate situations, performance increase can reach notable values, up to 45%.
Looking at current profiles, the explanation of this fact can be given in terms of the active
range introduced before. In fact stem size (10x40m) is comparable with the active range
when i
0
V
is in the interval 10
7
-10
8
A/m
3
, which shows the best performances. In other
words, when stem size is of the same order as the active range, structured contact
enhances the performance. On the contrary, a flat contact surface is preferable with
respect to a structured one when active range and stem size are not in the same order of
magnitude.
The geometry chosen stems are tall and closely packed does not fit for very low
values of q (2-4), in the whole range of i
0
V
.

17


Fig. 3: Comparison of flat and structured contacts, changing q (x-axis) and i
0
V

(different curves). On y-axis the relative variation of total current drawn 100*(i - i
ref
)/i
ref

is plotted.

3.3 Variation of geometry
The second set of simulations performed was aimed to study the effect of different
contact geometries, fixing q and varying i
0
V
. To achieve this, we decided to vary
independently the two parameters h and d, representing respectively the height of a stem
and half distance between two adjacent stems, as shown in fig. 1. The range of variation
chosen was 10-40m for both h and d and 10
5
-10
8
A/m
3
for i
0
V
.
Results are presented in fig. 4 where the percentage of current variation with respect to
the reference case (100*(i - i
ref
)/i
ref
) is plotted as a function of h. Each curve corresponds
to a different value of d, as specified in the legends. Lets come to the analysis of results
now. In fig 4 A results for i
0
V
= 10
5
A/m
3
are shown, close to the kinetic regime. Again
18

we notice that performances are always worse than the reference, and that configurations
with taller (high h) or more compact (low d) stems show lower current densities. This
confirms the interpretation given before, i.e. active volume is the limiting factor and is
reduced by the presence of stems.



Fig. 4: Comparison of flat and structured contacts, changing h (x-axis) and d
(different curves). On y-axis the relative variation of total current drawn 100*(i - i
ref
)/i
ref

is plotted. Each frame shows results at fixed i
0
V
.

Fig. 4 B shows results for i
0
V
= 10
6
A/m
3
. It is clear the existence of an optimal choice of
the geometrical parameters d and h. This is certainly due to a balancing between
i
0
V
= 10
5
A/m
3

i
0
V
= 10
6
A/m
3

i
0
V
= 10
7
A/m
3
i
0
V
= 10
8
A/m
3

A
D
B
C
19

contrasting factors, i.e. the positive effect of enhanced ion conduction and the negative
one of active volume reduction. A more detailed set of simulations should be made here
to catch the true maximum in the range d = 15-20m and h = 30-40m.
Fig. 4 C and D show results for i
0
V
= 10
7
and 10
8
A/m
3
, respectively. In both cases the
best performances were obtained for high stems (h = 40m), close each other (d = 10m),
for the range of geometries explored. The next step here would be to study the
optimization of stems size for this range of i
0
V
, outside the studied range for d and h.
Why the presence of the stem, in certain situations, can increase the assembly
performance? The answer is that the stem acts as an highway for ions, being a dense
conductor leaning out into the porous electrode. It reduces the ohmic loss due to ionic
transport in the electrode. This interpretation explains also why stems are useful when
their size is comparable with the active range (see sec. 3.1 above). We already discussed
the case of kinetic regimes, where active volume reduction due to the stem reduces
performances, and of ohmic regimes, where stems simply increase the mean electrolyte
thickness. In intermediate situations instead, losses due to ionic transport in the electrode
are comparable with activation and other ohmic losses. Their reduction due the patterned
contact significantly reflects into a performance increase.
A topic which will undergo future research efforts is the shape of stems. Intuition
suggests that triangular sections should work better than the rectangular ones studied here.
Moreover, the possibility of pyramidal stems, even more promising, should be tested in
3D simulations.


4. Conclusions
20


Comparison between reference flat electrolyte-electrode contact and structured one,
by means of steady state finite elements simulations, showed that the structured contact
can be advantageous in certain situations. In particular the selected stem geometry can
improve the current density in a wide range of kinetic parameter i
0
V
and of effectiveness
factor. Further studies are needed to reach the optimization of the stem geometry with
different i
0
V
.
List of symbols
a Activity
B permeability m
-2

c
i
volume concentration of species i mol m
3

c
T
total volume concentration in gas phase mol m
-3

E
cell
cell voltage V
D
K
i
Knudsen diffusion coefficient for species i m
2
s
-1

D
m
ij
Binary diffusion coefficient for species i and j m
2
s
-1

d diameter M
F Faraday constant C mol
-1

G Gibbs enthalpy J mol
1

i macroscopic current density A m
-2

V
i
microscopic current density A m
-3

0
V
i
exchange current density A m
-3

aTPB

active TPB length density m
-2

M
i

molar mass of species i Kg mol
-1

21

i
N

molar flux of species i mol m
2
s
1

p
i
partial pressure of species i Pa
q effectiveness factor
<r> mean pore radius M
R
g
gas constant J mol
-1
K
-1

T temperature K
S
i
source term for species i mol m
-3
s
-1

y
i
molar fraction of species i
Greek letters
o
transfer coefficient
o
thickness M
c
porosity
|
electric potential V
q
overpotential V

viscosity Pa s

resistivity O m
o
conductivity S m
-1

t
tortuosity
Subscripts
eq equilibrium
eff effective
i chemical species i
ref reference simulation
22



Acknowledgments

Many thanks to Cristiano Nicolella (Dep. of Chemical Engineering, Universit di Pisa,
via Diotisalvi 2, 56126 Pisa, Italy) for giving the possibility to perform the simulations in
his institute.

References

[1] Tao S., Irvine J.T.S., Kilner J.A. (2005) Adv Mater 17:1734-1737
[2] Pea-Martinez J., Marrero-Lpez D., Prez-Coll D., Ruiz-Morales J.C., Nuez P.
(2007) Electrochim Acta 52:2950-2958
[3] Shao Z., Haile S.M., (2004) Nature 431:170-173
[4] Lv H., Tu H-y, Zhao B-y, Wu Y-j, Hu K-a (2007) Solid State Ionics 177:3467-3472
[5] Tao S., Irvine J.T.S. (2003) Nat Mater 2:320-323
[6] Goodenough J.B., Huang Y-H (2007) Journal of Power Sources 173:1-10
[7] H. Zhu, R.J. Kee, V.M. Janardhanan, O. Deutschmann, D.G. Goodwin, J.Electrochem.
Soc. (2005) 152:A2427
[8] Jung H.Y., Kim W-S, Choi S-H, Kim H-C, Kim J, Lee H-W, Lee J-H (2006) Journal
of Power Sources 155:145-151
[9] Piccardo P., Chevalier S., Molins R., Viviani M., Caboche G., Barbucci A., Sennour
M., Amendola R. (2006) Surf Coat Technol 201(7):4471-4475
[10] Ucida H., Arisaka S., Watanabe M. (2000) Solid State Ionics 135:347-351
23

[11] Liu M., Dong D., Peng R., Gao J., Diwu J., Liu X., Meng G. (2008) Journal of
Power Sources 180:215-220
[12] Hammouche A., Schouler E.J.L., Henault M. (1988) Solid State Ionics 28-
30:1205-1207
[13] Mizusaki J., Tagawa H., Naraya K., Sasamoto T. (1991) Solid State Ionics
49:111-118
[14] Takeda Y., Sakaki Y., Ichikawa T., Imanishi N., Yamamoto O., Mari M., Mori N.,
Abe T. (1994) Solid State Ionics 72:257-264
[15] A. Bertei, C. Nicolella, F. Delloro, W. G. Bessler, N. Bundschuh, A. Thorel
(2011) ECS Trans. 35(1):883
[16] Y. Ji, J.A. Kilner, M.F. Carolan (2005) Solid State Ionics 176:937943
[17] B. Todd, J.B. Young (2002) Journal of Power Sources 110:186-200
[18] D. Arnost, P. Schneider (1995) The Chemical Engineering Journal 57:91-99