energies

Article
Thin Solid Film Electrolyte and Its Impact
on Electrode Polarization in Solid Oxide Fuel Cells
Studied by Three-Dimensional Microstructure-Scale
Numerical Simulation
Tomasz A. Prokop 1, * , Grzegorz Brus 1,† , Shinji Kimijima 2,† and Janusz S. Szmyd 1,†
1 Department of Fundamental Research in Energy Engineering, AGH University of Science and Technology,
30-059 Krakow, Poland; [email protected] (G.B.); [email protected] (J.S.S.)
2 Department of Machinery and Control Systems, Shibaura Institute of Technology, 135-8548 Tokyo, Japan;
[email protected]
* Correspondence: [email protected]; Tel.: +48-12-617-5053
† These authors contributed equally to this work.




Received: 31 August 2020; Accepted: 25 September 2020; Published: 1 October 2020


Abstract: In this work, a three-dimensional microstructure-scale model of a Solid Oxide Fuel Cell’s
Positive-Electrolyte-Negative assembly is applied for the purpose of investigating the impact of
decreasing the electrolyte thickness on the magnitude, and the composition of electrochemical losses
generated within the cell. Focused-Ion-Beam Scanning Electron Microscopy reconstructions are used
to construct a computational domain, in which charge transport equations are solved. Butler–Volmer
model is used to compute local reaction rates, and empirical relationships are used to obtain local
conductivities. The results point towards three-dimensional nature of transport phenomena in thin
electrolytes, and electrode-electrolyte interfaces.

Keywords: solid oxide fuel cell; electrolyte; focused ion beam scanning electron microscopy;
simulation; microstructure

1. Introduction
As highly efficient energy conversion devices, capable of converting the chemical energy of fuels,
directly into electrical energy, Solid Oxide Fuel Cells (SOFCs) have received much interest from the
contemporary scientific community [1–3]. One major source of thermodynamic losses generated
in an SOFC device is the resistance to the conduction of oxide ions in the ceramic phase of the
cell. Thus, it is usually beneficial from a thermodynamic standpoint to decrease the thickness of the
electrolyte. Fabrication methods, such as thermal spray, chemical vapor deposition, physical !vapor
deposition or pulsed laser deposition allow introducing electrolytic films thinner than 10 µm (as thin
as 1 µm) [4,5]. Another interesting technique, known as Plasma-Enhanced Atomic Layer Deposition
has allowed fabrication of electrolytes thinner than 0.1 µm. A number of research teams have
managed to manufacture working anode-supported SOFC cells with electrolytes as thin as 10 nm,
close to the theoretical limit related to dielectric breakdown strength [6]. Baek et al. [6] achieved
such low thickness by applying freestanding electrolyte membranes, manufactured using silicon
micromachining techniques. Nonetheless, cells with such a thin electrolyte face a number of problems,
including durability issues, such as cracks and pinholes, as well as local irregularities stemming from
the grain shape [6].
The understanding of the impact of local microstructure inhomogeneities may be improved
through the application of a three-dimensional, pore-scale model of charge transport phenomena.

Energies 2020, 13, 5127; doi:10.3390/en13195127 www.mdpi.com/journal/energies

Energies 2020, 13, 5127 2 of 14

In a non-continuous model of a porous electrode, each node within the computational domain is
assigned a specific transport function equivalent to the transport function of the corresponding phase
within the porous ceramic-metal electrode. The required digital distribution of phases within the
microstructure may be obtained using statistical modeling, or three-dimensional reconstructions of
a real electrode, obtained using nanotomography. The latter approach has been made possible by
the application of Focused Ion Beam Scanning Electron Microscopy (FIB-SEM) to analyze SOFC
electrodes, starting with the ground-breaking works by Wilson et al. [7]. In FIB-SEM nanotomography,
layers of sample are removed using a Focused Ion Beam. Subsequently, the uncovered cross-sections
are imaged using a scanning electron microscope. A non-continuous electrode model was applied
in works by Suzue et al. [8] (for a stochastic reconstruction), Shikazono et al. [9], Kanno et al. [10],
Kishimoto et al. [11] Carraro et al. [12] and others.
So far, modeling efforts regarding the impact of electrolyte thickness on electrochemical
performance have been relatively sparse, as most works focus on the thermo-mechanical properties of
the device [13,14]. A study by Park et al. [15] involved a continuous-electrode, cell-scale model of an
SOFC unit cell, including cases with electrolyte thickness ranging from 80 µm to 100 µm. Despite the
sensitivity to local microstructure irregularities, as well as the potential to introduce hierarchical
microstructures, and mesoscale electrode-electrolyte interface structures [16], the efforts to model
electrochemical transport phenomena at microstructure scale have been rare. Iwai et al. [17] used
a pore-scale model with non-continuous phase distribution to compare loss generation for several
electrolyte thicknesses ranging from 1 µm to 10 µm, demonstrating increase in standard deviation of
current densities within electrolytes thinner than 5 µm. Additionally, some influence of cathodic grain
size was predicted by the model. Kishimoto et al. [18] furthered this analysis to discuss the impact of
fabricating a grooved electrolyte.
The literature survey unravels the knowledge gap and the necessity to investigate the impact
of thin electrolyte on microstructure-scale potential distributions in an SOFC. The purpose of the
research described in this work is to apply a microstructure-scale, non-continuous model to analyze the
impact of cross-electrolyte effects on the performance of an SOFC positive-electrolyte-negative (PEN)
assembly more extensively. The cross-electrode interactions are present if there exists a difference
between the performance predicted by a model which assumes uniform current flow through the
electrolyte, and the model in which no such assumption is made. Previously we have constructed and
validated a model of this type, using it for the purpose of analyzing heterogeneity [19,20], and aging
phenomena [21]. However, in these studies, the electrodes were investigated individually. In our
current paper, we present a model, where the computational domain includes both the anode, and the
cathode. The decomposition of thermodynamic losses is presented alongside with distribution of
charge transfer and electric potential within a cell. A number of Positive-Electrolyte-Negative assembly
simulations with different electrolyte thicknesses are studied, and effects of cross-electrode phenomena
are assessed in particular.

2. Mathematical and Numerical Model

The modeled chemical reaction includes the half reactions of the combustion of hydrogen:

O2 + 2e− → 2O2−
(
on the anode,
2H2 + O2 → 2H2 O 2− −
(1)
2H2 + 2O → 2H2 O + 2e on the cathode.

The charge transfer model is based on the Poisson conservation with the source term being the
local, volumetric charge transfer rate:

(O2− )
(

−i , j = ion
∇ · σj ∇φj = , (2)
i , j = el (e− )
Energies 2020, 13, 5127 3 of 14

where φel (V) and φion (V) are the electron- and the ion- conducting phase potentials. σel and
σion (Ω−1 m−1 ) are the conductivities in electron-conducting and oxide phases, while i A m−3

is the volumetric charge transfer rate, computed using the Butler–Volmer model:

αano Fηact
 ano
αfrw Fηact
    
ano ano

 i0,tpb ltpb exp − exp − bcw anode
RT RT


i= "
cath
!
cath
!# , (3)
cath cath cath cath 2α Fη act 2α Fη act
frw
− exp − bcw

 (i0,dpb Adpb + i0,tpb ltpb ) exp cathode


RT RT

i0,tpb (A m−1 ) is the specific Triple Phase Boundary (TPB) exchange current density, i0,dpb (A m−2 )
is the specific Double Phase Boundary (DPB) density, ltpb (m m−3 ) is the local Triple Phase Boundary
(TPB) length density, and Adpb (m2 m−3 ) is the local Double Phase Boundary (DPB) length density.
Superscripts ‘ano’ and ‘cath’ indicate values specific to the anodic and the cathodic reaction sites
respectively. αfwd and αbcw are, respectively, forward and backwards charge transfer coefficient.
The specific exchange current densities and charge transfer coefficients for the cathode are taken from
studies by Matsuzaki et al. [22], Miyoshi et al. [23], and Kim et al. [24]. For the anode, the data by
de Boer [9,25] and Holtappels et al. [26] is used. The model assumes non-continuous distribution
of phases within the computational domain, which is based directly on the digital reconstructions
from FIB-SEM nanotomography. The TPB and the DPB densities are non-zero only in proximity to
the reaction sites. For the sake of simplicity, gas composition gradient in the active layer is neglected.
The cell analyzed for the purpose of this study is composed of a Ni-YSZ anode, a YSZ electrolyte,
and an LSCF-GDC cathode. The microstructure reconstructions discussed in this paper are taken from
our previous study [21,27]. The parameters of each electrode are discussed in Table 1. In our previous
research, the transport phenomena on each of the electrodes were solved individually, while the
potential drop on the electrolyte was computed separately. Such an approach was sufficient when
thicker electrolytes were analyzed. To perform simulations involving thin electrolytes, the model
was expanded and tailored to analyze cathodic, and anodic transport phenomena within a single
three-dimensional computational domain. Fourier boundary condition is implemented. For each
data point, a potential difference is designated arbitrarily. Open Circuit Voltage (OCV) is set as
the reference potential value at the electron-conducting boundary of the anode. The scheme of the
boundary conditions is visualized in Table 2. The set of equation is discretized using the Finite Volume
Method (FVM), with each finite volume in the computational domain corresponding to a voxel within
the source FIB-SEM data. The solution is achieved using Successive Over-Relaxation (SOR) method
with local linearization of the source term. The methodology is visualized in Figure 1. The resulting
model equations are solved using in-house numerical code written in C++.

Figure 1. A microstructure-scale non-continuous model of an SOFC electrode.

Energies 2020, 13, 5127 4 of 14

Table 1. The electrode parameters used in the cross-electrode phenomena study.

Phase

ltpb (m m−3 ) Adpb (m2 m−3 ) ψf τ ψf τ ψf τ

Cathode
LSCF GDC Pore
2.11 × 1012 1.63 × 106 0.34 3.95 0.27 10.32 0.39 2.42
Anode
Ni YSZ Pore
4.44 × 1012 0.349 4.95 0.462 3.009 0.159 26.159

Table 2. Scheme of Fourier Boundary Conditions—Positive-Electrolyte-Negative assembly of a cell

composed of a Nickel (Ni)—Yttrium Stabilized Zirconia (YSZ) anode, YSZ electrolyte, and Lanthanium
Strontium Cobaltite Ferrite (LSCF)—Gadolinum Doped Ceria (GDC) cathode.

zb,an

jel= F[ϕel,b,an-ϕel(z1)]
jion= F[ϕion,b,an-ϕion(z1)]
Ji= F[pi,an-p(z1)]

z1
Pore
x
Ni
y
z YSZ
GDC

Pore

LSCF
z2
jel= F[ϕel,b,an-ϕel(z2)]
jion= F[ϕion,b,an-ϕion(z2)]
Ji= F[pi,an-p(z2)]
zb,cat

z = zb,an z = zb,cat x ∈ {0, xb } y ∈ {0, yb }

∂pH2 ∂pH2
pH2 = pH2 ,b ∂x = 0 ∂y = 0
∂pH2 O ∂pH2 O
pH2 O = pH2 O,b ∂x = 0 ∂y = 0
∂pO2 ∂pO2
pO2 = pO2 ,b ∂x = 0 ∂y = 0
∂pN2 ∂pN2
pN2 = pN2 ,b ∂x = 0 ∂y = 0
∂φel ∂φel
φel = φel,an ∂x = 0 ∂y = 0
∂φel ∂φel
φel = φel,cat ∂x = 0 ∂y = 0

More details regarding the mathematical, and the numerical model have been provided in our
previous papers [19,21].

3. Results
Validation of the model was carried out by comparison of simulated cell potential to experimental
data. The experimental data—microstructure reconstructions and current-voltage relationships—were
obtained from a SolidPower S.p.A stack, consisting of 9 anode-supported cells connected in
series [27,28]. Each cell consisted of a 240 µm Ni-YSZ anode, 10 µm Yttria-stabilized Zirconia (YSZ)
electrolyte and a 50 µm Gadolinium Doped Ceria (GDC) Lanthanum Strontium Cobalt Ferrite (LSCF)
cathode. The microstructure data is presented in Table 1. The results of the comparison are presented
in Figure 2. For each data point in the chart, the horizontal coordinate corresponds to the voltage
Energies 2020, 13, 5127 5 of 14

measured experimentally, while the vertical coordinate refers to the voltage predicted by the model.
A point located on the line represents perfect fit. While some discrepancy between the experiment
and the simulation is observed, the agreement is satisfactory, given that reaction and conduction
parameters are taken from open literature, rather than being fitted to the experimental data set.

Comparison to experimental data Comparison to experimental data

1.2 1.2

1.1 1.1
Predicted voltage (V)

Predicted voltage (V)

1.0 1.0

0.9 0.9

0.8 0.8
T= 1073 K pH2 = 70 kPa
T= 1023 K pH2 = 60 kPa
0.7 T= 973 K 0.7 pH2 = 50 kPa
T= 923 K pH2 = 40 kPa
exp = exp =

0.6 0.7 0.8 0.9 1.0 1.1 1.2 0.6 0.7 0.8 0.9 1.0 1.1 1.2
Measured voltage exp (V) Measured voltage exp (V)

(a) (b)
Figure 2. Comparison of the experimental cell voltage to simulated cell voltage. Simulation parameters:
p = 100,000 Pa, pH2 = 60, 000 Pa, T = 1023 K, unless specified otherwise. Parametric study involving
various: (a) Temperatures (b) Pressures.

The simulation was performed for virtual electrodes with different electrode thicknesses, ranging
from 0.10 µm to 10 µm. Since a thin portion of the anode is considered, concentration losses are
neglected. The simulated polarization curves are presented in Figure 3. Each line in Figure 3
represents combined thermodynamic losses on positive and negative electrodes for the given current.
The distribution of curves demonstrates that as the electrolyte becomes thinner, the losses decrease,
although not at linear rate. This indicates the need to explore this phenomenon to a greater extent.
Three-dimensional current and potential distributions are presented in Figure 4. It can be seen that
most of the voltage loss due to reaction and conduction irreversibility occurs on the side of the cathode,
which is relatively thin, considering the thickness of the active layer (indicated by potential gradient).
The irregular distribution of ion-conducting phases has a greater impact on the cathode-side, due
to high tortuosity of the GDC phase, and lower conductivity of the LSCF phase for this particular
electrode. On the other hand, the YSZ phase has high volume ratio and relatively low tortuosity,
which—combined with high TPB density—results in lower electrochemical losses.

0.100
0.12 0.099
0.098
0.097
0.10 0.096
excluding electrolyte losses
Total overpotential (V)

0.095
2100 2150 2200 2250 2300
0.08

0.06

0.04 Electrolyte thickness: 0.1 m

Electrolyte thickness: 0.25 m
Electrolyte thickness: 0.5 m
Electrolyte thickness: 1 m
0.02 Electrolyte thickness: 2.5 m
Electrolyte thickness: 5 m
Electrolyte thickness: 10 m
0 500 1000 1500 2000 2500 3000
Mean current density j (Am 2)

Figure 3. Total cell overpotential. Simulation parameters: p = 100,000 Pa, pH2 = 20,000 Pa,
pO2 = 21,000 Pa, T = 1023 K.
Energies 2020, 13, 5127 6 of 14

(a) η = 0.12 V.

(b) η = 0.06 V.
Figure 4. The distribution of ion-conducting phase potential in cells with different electrolyte
thicknesses. Simulation parameters: p = 100,000 Pa, pH2 = 20,000 Pa, pO2 = 21,000 Pa, T = 1023 K.
Energies 2020, 13, 5127 7 of 14

The generated overpotential can be divided into the components related losses occurring due
to the resistance to ionic current conduction (ohmic overpotential), and the reaction irreversibilities
(activation overpotential) for each electrode. The results of overpotential decomposition, together with
the active layer potential and current distributions are presented in Figure 5. The subfigures presented
on the left-hand side show the charge transfer rate between ion and electron-conducting phases
(blue line), together with the electric potential of the ion-conducting phase (red dash-dot line), and the
electron-conducting phase (red dotted line) for the total overpotential η = 0.125 V. It can be seen that
the charge transfer rate distributions are jagged and uneven—this is related to the non-continuous
distribution of reaction sites within the computational domain, which is based on the FIB-SEM sample.
In the electrodes, the potentials distributions are similar to hyperbolic tangent functions, while in
the electrolyte, they appear to be linear. This holds true for all the discussed distributions. As the
electrolyte thickness is decreased, the current generated on the electrodes increases. Subfigures
presented on the right-hand side show decomposition of total overpotential for different current
densities. The maximum current density marked on the charts refers to the highest current density
obtained within the presented simulations. For the analyzed current range, the losses appear to
be divided more or less evenly between the reaction irreversibilities (activation overpotential) and
conduction irreversibilities (ohmic overpotential). Cathodic losses are roughly twice as high as the
losses on the anode. Most of the change in polarization appears to be the result of changing electrolyte
overpotential. While, as expected, the share of losses generated on the electrolyte shrinks as its
thickness decreases, the composition of losses generated on the electrodes does not change significantly.
It is clear that as the electrolyte becomes thicker, less current is generated, and more of the overpotential
is devoted to overcoming the resistivity of the electrolyte. This is the reason why a comparative analysis
of the problem requires normalization, for example by means of isolating and subtracting the losses
generated on the electrolyte.
Figure 6 displays the potential losses generated on the electrolyte for cells of different thicknesses.
As the cell becomes thinner, the losses do not decrease linearly. For every case, halving the electrolyte
thickness results in a smaller decrease of overpotential, than what would be expected from a
linear relationship. The three-dimensional model is able to predict this behavior since it captures
the horizontal component of the current flow through the electrolyte, occurring due to the local
irregularities in the microstructure at the electrode-electrolyte interface. Thus, the pathways for the
ionic current are longer than what would be expected from a one-dimensional model assuming that
the ionic current in the electrolyte flows in a straight line.

0.925 0.13
Charge transfer rate i (107 Am 3)

Max. current density

100 0.12
0.900 0.11 Cathodic
75
Electric potential (V)

0.10 activation
overpotential
overpotential (V)

0.875 50
Decomposition of

0.09
0.850 25 0.08
0.07
0 Cathodic ohmic
0.825 0.06 overpotential
25 0.05
0.800 50 0.04
ion - ion conductor potential
Electrolyte
75 0.03 overpotential
0.775 el - electron conductor potential
0.02 Anodic ohmic
100 overpotential
0.750 i - Charge transfer rate 0.01 Anodic activation
overpotential
0 10 20 30 40 50 60 0 500 1000 1500 2000 2500 3000
Distance from the anodic channel z ( m) Mean current density j (Am 2)
(a)

Figure 5. Cont.
Energies 2020, 13, 5127 8 of 14

0.925 0.13

Charge transfer rate i (107 Am 3)

Max. current density
100 0.12
0.900 0.11
Electric potential (V) 75 Cathodic
0.10 activation
overpotential

overpotential (V)
0.875 50

Decomposition of
0.09
0.850 25 0.08
0.07
0
0.825 0.06 Cathodic ohmic
25 overpotential
0.05
0.800 50 0.04
ion - ion conductor potential 0.03 Electrolyte
75 overpotential
0.775 el - electron conductor potential
0.02 Anodic ohmic
100 overpotential
0.750 i - Charge transfer rate 0.01 Anodic activation
overpotential
0 10 20 30 40 50 60 0 500 1000 1500 2000 2500 3000
Distance from the anodic channel z ( m) Mean current density j (Am 2)
(b)

0.925
Charge transfer rate i (107 Am 3) 0.13
Max. current density
100 0.12
0.900 0.11
75 Cathodic
Electric potential (V)

0.10 activation
overpotential (V)
0.875 50 overpotential
Decomposition of

0.09
0.850 25 0.08
0.07
0
0.825 0.06 Cathodic ohmic
25 0.05 overpotential
0.800 50 0.04
ion - ion conductor potential 0.03 Electrolyte
0.775 el - electron conductor potential 75 overpotential
0.02 Anodic ohmic
overpotential
100 0.01 Anodic activation
0.750 i - Charge transfer rate overpotential
0 10 20 30 40 50 60 0 500 1000 1500 2000 2500 3000
Distance from the anodic channel z ( m) Mean current density j (Am 2)
(c)

0.925 0.13
Charge transfer rate i (107 Am 3)

Max. current density

100 0.12
0.900 0.11
75 Cathodic
Electric potential (V)

0.10 activation
overpotential (V)

0.875 50 overpotential
Decomposition of

0.09
0.850 25 0.08
0.07
0
0.825 0.06
25 Cathodic ohmic
0.05 overpotential
0.800 50 0.04
ion - ion conductor potential 0.03 Electrolyte
0.775 el - electron conductor potential 75 overpotential
0.02 Anodic ohmic
overpotential
100 0.01 Anodic activation
0.750 i - Charge transfer rate overpotential
0 10 20 30 40 50 60 0 500 1000 1500 2000 2500 3000
Distance from the anodic channel z ( m) Mean current density j (Am 2)
(d)
Figure 5. Cont.
Energies 2020, 13, 5127 9 of 14

0.925 0.13

Charge transfer rate i (107 Am 3)

Max. current density
100 0.12
0.900 0.11
Electric potential (V) 75 Cathodic
0.10 activation

overpotential (V)
0.875 50 overpotential

Decomposition of
0.09
0.850 25 0.08
0.07
0
0.825 0.06
25 Cathodic ohmic
0.05 overpotential
0.800 50 0.04
ion - ion conductor potential 0.03 Electrolyte
0.775 el - electron conductor potential 75 overpotential
Anodic ohmic
0.02 overpotential
100 0.01 Anodic activation
0.750 i - Charge transfer rate overpotential
0 10 20 30 40 50 60 0 500 1000 1500 2000 2500 3000
Distance from the anodic channel z ( m) Mean current density j (Am 2)
(e)

0.925
Charge transfer rate i (107 Am 3) 0.13
Max. current density
100 0.12
0.900 0.11
75
Electric potential (V)

0.10 Cathodic
activation
overpotential (V)
0.875 50 overpotential
Decomposition of

0.09
0.850 25 0.08
0.07
0
0.825 0.06
25 Cathodic ohmic
0.05 overpotential
0.800 50 0.04
ion - ion conductor potential 0.03 Electrolyte
0.775 el - electron conductor potential 75 overpotential
Anodic ohmic
0.02 overpotential
100 0.01 Anodic activation
0.750 i - Charge transfer rate overpotential
0 10 20 30 40 50 60 0 500 1000 1500 2000 2500 3000
Distance from the anodic channel z ( m) Mean current density j (Am 2)
(f)

0.925 0.13
Charge transfer rate i (107 Am 3)

Max. current density

100 0.12
0.900 0.11
75
Electric potential (V)

0.10 Cathodic
activation
overpotential (V)

0.875 50 overpotential
Decomposition of

0.09
0.850 25 0.08
0.07
0
0.825 0.06
25 Cathodic ohmic
0.05 overpotential
0.800 50 0.04
ion - ion conductor potential 0.03 Electrolyte
0.775 el - electron conductor potential 75 overpotential
Anodic ohmic
0.02 overpotential
100 0.01 Anodic activation
0.750 i - Charge transfer rate overpotential
0 10 20 30 40 50 60 0 500 1000 1500 2000 2500 3000
Distance from the anodic channel z ( m) Mean current density j (Am 2)
(g)
Figure 5. Overpotential decomposition together with active layer potential and current distributions
for η = 0.125 V. Simulation parameters: p = 100,000 Pa, pH2 = 20,000 Pa, pO2 = 21,000 Pa, T = 1023 K.
Electrolyte thickness: (a) 10 µm; (b) 5 µm; (c) 2.5 µm; (d) 1 µm; (e) 0.5 µm; (f) 0.25 µm; (g) 0.1 µm.
Energies 2020, 13, 5127 10 of 14

0.020
Electrolyte thickness: 0.1 m
0.018 Electrolyte thickness: 0.25 m
Electrolyte thickness: 0.5 m
0.016 Electrolyte thickness: 1 m
Electrolyte thickness: 2.5 m
Electrolyte thickness: 5 m
0.014 Electrolyte thickness: 10 m
on the electrolyte

0.012
Losses (V)

0.010

0.008

0.006

0.004

0.002

0 300 600 900 1200 1500 1800 2100 2400 2700 3000
Mean current density j (Am 2)
Figure 6. The voltage losses on the electrolyte. Simulation parameters: p = 100,000 Pa, pH2 = 20,000 Pa,
pO2 = 21,000 Pa, T = 1023 K.

As suggested before, the overpotential-current relationships derived from the three-dimensional

model were corrected to remove the influence of electrolyte by subtracting the averaged potential
drop on the electrolyte itself. The corrected polarization curves are presented in Figure 7. Interestingly,
the following tendency emerges: the thinner the electrolyte, the smaller the losses on the electrodes.
The differences among the polarization curves excluding the losses on the electrolyte suggest the
presence of cross-electrode phenomena, as the electrode-only losses decrease to a greater extent than
what would be expected from just the removal of the current-resisting material. Most of the difference
is present on the cathode side, while the anodic losses remain mostly unaffected by the decrease of
electrolyte thickness, increasing slightly as the distance to the cathode shrinks.

0.12 0.100 ic
0.099 thod
0.098 odic + ecnatial
0.097 An erpot
0.096 ov
0.10 0.095
2100 2150 2200 2250 2300
excluding electrolyte losses

0.08
Total overpotential (V)

tial
0.06 oten
erp
odi c ov
cath
ic +
0.04
A nod

al
0.02 Anodic overpotenti
Electrolyte thickness: 0.1 m
Electrolyte thickness: 0.25 m
0.00 0.0215 Electrolyte thickness: 0.5 m
al Electrolyte thickness: 1 m
0.0213
0.0211 Anodic overpotenti Electrolyte thickness: 2.5 m
0.0209 Electrolyte thickness: 5 m
0.02 0.0207
0.0205 Electrolyte thickness: 10 m
2100 2110 2120 2130 2140 2150
0 500 1000 1500 2000 2500 3000
Mean current density j (Am 2)
Figure 7. The cell overpotential excluding the voltage losses on the electrolyte. Simulation parameters:
p = 100,000 Pa, pH2 = 20,000 Pa, pO2 = 21,000 Pa, T = 1023 K.
Energies 2020, 13, 5127 11 of 14

Figure 8 depicts the share of overpotential components for different electrolyte thicknesses.
Intuitively, the total share of losses generated on the electrodes is smaller for thicker electrolytes.
However, the model predicts that the share of losses on the cathode, in relation to electrode-only losses,
increases. Thus, it appears that the majority of electrode performance improvement in the presence
of a thin electrolyte occurs on the cell’s cathode. These tendencies are observed for both the ohmic
overpotential, and the activation overpotential. The results from Figure 7, and Figure 8 indicate that
the electrolyte thickness indirectly impacts the overpotential by value higher than its ohmic resistance.

Percentage of losses on electrodes

anodic activation overpotential anodic activation overpotential
70%
anodic ohmic overpotential 70%
anodic ohmic overpotential
Percentage of total losses

electrolyte overpotential cathodic ohmic overpotential

cathodic ohmic overpotential cathodic activation overpotential
60% cathodic activation overpotential 60%

50% 50%

40% 40%

30% 30%

20% 20%

1.0 2.5 5.0 10.0 1.0 2.5 5.0 10.0

Electrolyte thickness ( m) Electrolyte thickness ( m)
(a) Total losses (b) Losses on the electrodes
Figure 8. Overpotential component at different electrolyte thicknesses as percentage of the total losses,
and as percentage of the losses on electrodes. Simulation parameters: p = 100,000 Pa, pH2 = 20,000 Pa,
pO2 = 21,000 Pa, T = 1023 K, η = 0.125 V.

4. Conclusions
In this paper, a microstructure-scale transport model was used to compute potential fields and
current distributions in a positive-electrolyte-negative assembly of a Solid Oxide Fuel Cell for several
cases, in which different electrolytes of different thicknesses were considered. The model has shown
that reducing electrolyte thickness below 10 µm yields diminishing returns, as its relationship to
overpotential is non-linear. This is likely an effect of a horizontal component to current pathways,
resulting from microstructure irregularities. Although, the results show little effect of electrolyte
thickness on the composition of electrode overpotential, some differences were predicted nonetheless.
The current-overpotential relationships indicated that the electrode-only losses would decrease as
the electrolyte became thinner. As the current density increases, so do the differences, although they
do not seem to exceed 5%. This behavior could be explained by the local phase distribution affecting
the pathways of ions penetrating from the cathode to the anode. Additionally, it was predicted
that the bulk of the increase happens on the cathode, as the anodic losses did decrease by a small
margin. These tendencies were observed for both the activation overpotential, and the ionic conduction
overpotential. The relatively low magnitude of the electrode-only losses may be related to the high
degree of homogeneity in the microstructures used for the purpose of this study. If the electrodes
were more anisotropic, for example featuring a periodic inhomogeneity, the cross-electrode effects
would likely be more prominent. The results suggest that transport on a thin electrolyte is not
one-dimensional, pointing towards the importance of using a three-dimensional microstructure scale
model to tackle this research problem. When electrodes are analyzed separately, which is common in
the literature, the obtained results may differ compared to the simulation conducted on the entire PEN
structure, particularly when the electrolyte is thinner than 10 µm.
Energies 2020, 13, 5127 12 of 14

Author Contributions: Conceptualization, S.K.; methodology, S.K.; software, T.A.P.; validation, T.A.P.;
formal analysis, T.A.P.; investigation, G.B.; resources, and G.B.; data curation, T.A.P.; writing—original draft
preparation, T.A.P.; writing—review and editing, G.B., S.K. and J.S.S.; visualization, T.A.P.; supervision, J.S.S. and
S.K.; project administration, G.B.; funding acquisition, G.B. All authors have read and agreed to the published
version of the manuscript.
Funding: This work was supported by the National Science Centre of Poland (SONATA-10, Grant No.
UMO-2015/19/D/ST8/00839). The authors are grateful for the support.
Acknowledgments: This research was supported in part by PL-Grid Infrastructure. Matplotlib Python library
was used for data visualization [29].
Conflicts of Interest: The authors declare no conflict of interest.

Nomenclature
Abbreviations
DPB Double Phase Boundary
FIB Focused Ion Beam
GDC Gadolinum Doped Ceria
LSCF Lanthanum Strontium Cobalt Ferrite
OCV Open Circuit Voltage
PEN Positive Electrolyte Negative
SEM Scanning Electron Microscopy
SOFC Solid Oxide Fuel Cell
SOR Successive Over-Relaxation
TPB Triple Phase Boundary
YSZ Yttrium-Stabilized Zirconia
Roman symbols
F Faraday constant A s mol−1
i Charge transfer rate A m−3
i0,tpb Equilibrium exchange current density at TPB A m−1
i0,dpb Equilibrium exchange current density at DPB A m−2
j Mean charge transfer rate A m−2
ltpb TPB density m m−3
Adpb DPB density m2 m − 3
R Universal gas constant J mol−1 K−1
T Temperature K
x,y Planar coordinates m
z Depth (distance from anodic channel) m
Greek symbols
α Charge transfer coefficient -
η Overpotential V
φ Electrical potential V
σ Electrical conductivity Ω −1 m−1
τ Tortuosity -
ψ Phase volume fraction -
Subscripts
act Activation
b Boundary (bulk)
dpb Double Phase Boundary
H2 Hydrogen
H2 O Water vapor
O2 Oxygen
i A substance
ion Oxide ion conducting phase
tpb Triple phase boundary
dpb Double phase boundary
0 Equilibrium
Energies 2020, 13, 5127 13 of 14

Superscripts
ano anodic
cath cathodic

References
1. Pianko-Oprych, P.; Hosseini, S.M. Dynamic Analysis of Load Operations of Two-Stage SOFC Stacks Power
Generation System. Energies 2017, 10, 2103. [CrossRef]
2. Fang, X.; Zhu, J.; Lin, Z. Effects of Electrode Composition and Thickness on the Mechanical Performance of a
Solid Oxide Fuel Cell. Energies 2018, 11, 1735. [CrossRef]
3. Gandiglio, M.; De Sario, F.; Lanzini, A.; Bobba, S.; Santarelli, M.; Blengini, G.A. Life Cycle Assessment of a
Biogas-Fed Solid Oxide Fuel Cell (SOFC) Integrated in a Wastewater Treatment Plant. Energies 2019, 12, 1611.
[CrossRef]
4. Coddet, P.; Liao, H.L.; Coddet, C. A review on high power SOFC electrolyte layer manufacturing using
thermal spray and physical vapour deposition technologies. Adv. Manuf. 2014. [CrossRef]
5. Noh, H.S.; Lee, H.; Kim, B.K.; Lee, H.W.; Lee, J.H.; Son, J.W. Microstructural factors of electrodes affecting
the performance of anode-supported thin film yttria-stabilized zirconia electrolyte (1 µm) solid oxide fuel
cells. J. Power Sources 2011. [CrossRef]
6. Baek, J.D.; Liu, K.Y.; Su, P.C. A functional micro-solid oxide fuel cell with a 10 nm-thick freestanding
electrolyte. J. Mater. Chem. A 2017, 5, 18414–18419. [CrossRef]
7. Wilson, J.R.; Kobsiriphat, W.; Mendoza, R.; Chen, H.Y.; Hiller, J.M.; Miller, D.J.; Thornton, K.; Voorhees, P.W.;
Adler, S.B.; Barnett, S.A. Three-dimensional reconstruction of a solid-oxide fuel-cell anode. Nat. Mater. 2006,
5, 541–544. [CrossRef]
8. Suzue, Y.; Shikazono, N.; Kasagi, N. Micro modeling of solid oxide fuel cell anode based on stochastic
reconstruction. J. Power Sources 2008, 184, 52–59. [CrossRef]
9. Shikazono, N.; Kanno, D.; Matsuzaki, K.; Teshima, H.; Sumino, S.; Kasagi, N. Numerical Assessment of
SOFC Anode Polarization Based on Three-Dimensional Model Microstructure Reconstructed from FIB-SEM
Images. J. Electrochem. Soc. 2010, 157, B665–B672. [CrossRef]
10. Kanno, D.; Shikazono, N.; Takagi, N.; Matsuzaki, K.; Kasagi, N. Evaluation of SOFC anode polarization
simulation using three-dimensional microstructures reconstructed by FIB tomography. Electrochim. Acta
2011, 56, 4015–4021. [CrossRef]
11. Kishimoto, M.; Iwai, H.; Saito, M.; Yoshida, H. Three-Dimensional Simulation of SOFC Anode Polarization
Characteristics Based on Sub-Grid Scale Modeling of Microstructure. J. Electrochem. Soc. 2012, 159, B315–B323.
[CrossRef]
12. Carraro, T.; Joos, J.; Rüger, B.; Weber, A.; Ivers-Tiffée, E. 3D finite element model for reconstructed
mixed-conducting cathodes: I. Performance quantification. Electrochim. Acta 2012, 77, 315–323. [CrossRef]
13. Baek, J.D.; Yoon, Y.J.; Lee, W.; Su, P.C. A circular membrane for nano thin film micro solid oxide fuel cells
with enhanced mechanical stability. Energy Environ. Sci. 2015. [CrossRef]
14. Chen, Z.; Wang, X.; Brandon, N.; Atkinson, A. Numerical Study of Solid Oxide Fuel Cell Contacting
Mechanics. Fuel Cells 2018, 18, 42–50. [CrossRef]
15. Park, J.; Kim, D.; Baek, J.; Yoon, Y.J.; Su, P.C.; Lee, S. Effect of Electrolyte Thickness on Electrochemical
Reactions and Thermo-Fluidic Characteristics inside a SOFC Unit Cell. Energies 2018, 11, 473. [CrossRef]
16. Chen, Y.; Zhang, Y.; Baker, J.; Majumdar, P.; Yang, Z.; Han, M.; Chen, F. Hierarchically oriented macroporous
anode-supported solid oxide fuel cell with thin ceria electrolyte film. ACS Appl. Mater. Interfaces 2014,
6, 5130–5136. [CrossRef] [PubMed]
17. Iwai, H.; Kadomiya, R.; Kishimoto, M.; Saito, M.; Yoshida, H. Numerical analysis of cross-electrode
interaction in SOFCs with thin electrolyte. In Proceedings of the 13th European SOFC & SOE Forum 2018,
Luzern, Switzerland, 3–6 July 2018.
18. Kishimoto, M.; Sasaki, M.; Iwai, H.; Yoshida, H. Numerical assessment of mesoscale modification of thin
electrolyte in anode-supported solid oxide fuel cells. In Proceedings of the 13th European SOFC & SOE
Forum 2018, Luzern, Switzerland, 3–6 July 2018; p. A1305.
Energies 2020, 13, 5127 14 of 14

19. Prokop, T.; Berent, K.; Iwai, H.; Szmyd, J.S.; Brus, G. A three-dimensional heterogeneity analysis of
electrochemical energy conversion in SOFC anodes using electron nanotomography and mathematical
modeling. Int. J. Hydrogen Energy 2018, 43, 10016–10030. [CrossRef]
20. Prokop, T.; Berent, K.; Szmyd, J.S.; Brus, G. A three-dimensional numerical assessment of heterogeneity
impact on a solid oxide fuel cell’s anode performance. Catalyst 2018, 8, 503. [CrossRef]
21. Prokop, T.A.; Berent, K.; Mozdzierz, M.; Szmyd, J.S.; Brus, G. A Three-Dimensional Microstructure-Scale
Simulation of a Solid Oxide Fuel Cell Anode—The Analysis of Stack Performance Enhancement After a
Long-Term Operation. Energies 2019, 12, 4784. [CrossRef]
22. Matsuzaki, K.; Shikazono, N.; Kasagi, N. Three-dimensional numerical analysis of mixed ionic and electronic
conducting cathode reconstructed by focused ion beam scanning electron microscope. J. Power Sources 2011,
196, 3073–3082. [CrossRef]
23. Miyoshi, K.; Miyamae, T.; Iwai, H.; Saito, M.; Kishimoto, M.; Yoshida, H. Exchange current model for
(La0.8 Sr0.2 )0.95 MnO3 (LSM) porous cathode for solid oxide fuel cells. J. Power Sources 2016, 315, 63–69.
[CrossRef]
24. Kim, Y.T.; Jiao, Z.; Shikazono, N. Evaluation of La0.6 Sr0.4 Co0.2 Fe0.8 O3−δ− Gd0.1 Ce0.9 O1.95 composite cathode
with three dimensional microstructure reconstruction. J. Power Sources 2017, 342, 787–795. [CrossRef]
25. de Boer, B. Hydrogen Oxidation at Porous Nickel and Nickel/Yttria Stabilised Zirconia Cermet Electrodes.
Ph.D. Thesis, Universiteit Twente, Enschede, The Netherlands, 1998.
26. Holtappels, P.; de Haart, L.G.J.; Stimming, U. Reaction of Hydrogen/Water Mixtures on Nickel-Zirconia
Cermet Electrodes: I. DC Polarization Characteristics. J. Electrochem. Soc. 1999, 146, 1620–1625. [CrossRef]
27. Brus, G.; Iwai, H.; Szmyd, J.S. An Anisotropic Microstructure Evolution in a Solid Oxide Fuel Cell Anode.
Nanoscale Res. Lett. 2020, 15, 3. [CrossRef] [PubMed]
28. Mozdzierz, M.; Berent, K.; Kimijima, S.; Szmyd, J.S.; Brus, G. A Multiscale Approach to the Numerical
Simulation of the Solid Oxide Fuel Cell. Catalysts 2019, 9, 253. [CrossRef]
29. Hunter, J.D. Matplotlib: A 2D graphics environment. Comput. Sci. Eng. 2007, 9, 90–95. [CrossRef]

c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access
article distributed under the terms and conditions of the Creative Commons Attribution
(CC BY) license (http://creativecommons.org/licenses/by/4.0/).