Abstract

A continuum model for yttria-stabilized zirconia (YSZ) in the framework of non-equilibrium thermodynamics is developed. Particular attention is given to (i) modeling of the YSZ-metal-gas triple phase boundary, (ii) incorporation of the lattice structure and immobile oxide ions within the free energy model and (iii) surface reactions. A finite volume discretization method based on modified Scharfetter-Gummel fluxes is derived in order to perform numerical simulations. The model is used to study the impact of yttria and immobile oxide ions on the structure of the charged boundary layer and the double layer capacitance. Cyclic voltammograms of an air-half cell are simulated to study the effect of parameter variations on surface reactions, adsorption and anion diffusion.

Keywords

Introduction

Detailed continuum models of high temperature solid oxide electrochemical cells (SOEC)1 describe the underlying chemistry with spatially distinguished phases (oxide ion conductor, electric conductor, gas) of the triple phase boundary [1, 2, 3, 4]. Surface physics processes such as tangential diffusion and surface chemical reactions of the surface species are employed. In particular, the electron-transfer reaction at the triple phase boundary is usually modelled with Butler-Volmer-type kinetics containing overpotential, the difference of the electric potential between the metal and the bulk of the YSZ, as the driving force. The ionically or electrically conductive parts of a solid oxide cell are electroneutral in the respective bulks. The overpotential, appearing at the phase interface, is caused by formation of a charged double layer of oxide ions in YSZ and electrons in the electrode. Although the overpotential correlates with the excess concentration of oxide ions available for the electron-transfer reaction in steady-state scenarios, it cannot capture the dynamics of the double layer. Therefore, if such a model is compared to the results of a dynamic current-voltage measurement, e.g., electrochemical impedance spectroscopy or linear-sweep voltammetry, the dynamics of the double layer is underrepresented.

To determine the structure and dynamics of the space-charge layer of oxide ions in the YSZ, a generalized Poisson-Nernst-Planck (PNP) system can employed. By the generalization it is possible to account for the effect of the finite density of available lattice sites for oxide ions at the continuum level.

Such an approach was already used to capture the formation and behavior of the electrochemical double layers at electrode-electrolyte interfaces [5, 6]. The PNP system was successfully applied to the solid-state electrochemical systems, e.g., lithium batteries [7, 8, 9]. In [10], the PNP equations were already applied for proton ceramic fuel cells; however, the thermodynamics of the crystalline lattice and of the surface were not taken into account.

In this work, we apply a modeling approach for charged bulk-surface interfaces based on first principles of nonequilibrium thermodynamics [11, 12]. The resulting generalized Poisson-Nernst-Planck system is used to formulate the model of dynamics of the space-charged layer at the YSZ-metal-air triple interface. The main advantage of this approach is its consistency between the free energy (equilibrium) and fluxes (dynamics).

The paper is organized as follows. The free energy model of the bulk YSZ, capturing the crystalline structure, immobile oxide ions and elastic deformation, is developed in the “Bulk YSZ” section. The resulting chemical potentials are introduced into the gPNP model [11] after its modification for the description of the lattice velocity. In the following section, bulk metal and bulk gas phases are treated under the assumption of diffusional equilibrium. The free energy of the surface and the surface dynamics are described and developed in the the “Surface—triple phase boundary” section. The modeling approach results in a coupled system of evolution equation describing the transport of oxide ions in the bulk of electrolyte, adsorption of oxide ions from bulk to the surface and electron-transfer reaction alongside with the Poisson equation.

Using a finite volume based discretization, double layer capacitance and linear-sweep voltammetry simulations are performed in the “Simulation of a SOC half-cell” section. The performed simulations study the effects of the newly introduced concept of immobile oxide ions, the free energy parameters and the kinetic rates on the current response.

The novelty of the approach lies in the synthesis of the crystalline lattice bulk-surface free energy description and the coupled bulk-surface dynamics in non-equilibrium thermodynamics framework. Owing to this, it is possible to simulate the equilibrium behavior, e.g., the double layer capacitance, and dynamic behavior, e.g., the cyclic voltammetry, using a single model. Notable contribution to the state of the art models of YSZ is the thermodynamic treatment of the surface dynamics.

Bulk YSZ

We consider the charge transport exclusively in the isothermal electrostatic setting, therefore the temperature T is assumed to be constant and the electric field is given as E = −∇φ. Moreover, a simple material model for polarization based on a constant susceptibility χ is chosen.

General mixture and crystalline structure

Mixture quantities

We model YSZ as mixture of four constituents: zirconium and yttrium cations denoted by Zr and Y, respectively, and oxide anions. We assume that only a part of the oxide anions is freely mobile and refer to these as Om, whereas the remaining immobile oxide anions Oi are fixed to the underlying crystal structure. For referencing the different constituents of the mixture we use the index set \(\mathcal {I}_{\text {YSZ}}=\{\text {Zr},\text {Y},\text {Oi},\text {Om}\}\). Each constituent is characterized by the atomic mass mα and its atomic charges zαe0, where \(\alpha \in \mathcal {I}_{\text {YSZ}}\). The constant e0 is the elementary charge and zα is the charge number of the constituent. Multiplication of the number densities nα by mα gives the partial mass densities,

Crystalline structure

The crystalline structure of pure ZrO2 is well known (see e.g. [13]) and might be described conveniently in terms of unit crystal cells. Unit crystal cells of yttria-doped zirconia are, due to the yttria doping, difficult to be described systematically [14].

To overcome this, we introduce cation and anion spatial lattices, so that they coincide with the respective lattices in pure cubic ZrO2, i.e., locations of Zr4+ or O2− (see Fig. 1). Contrary to the pure ZrO2, the cation lattice of YSZ is occupied also by Y3+ and some of the anion lattice sites may be empty. The cation lattice unit cell is assumed to be face-centered cubic and contains 8 cations in its vertices and 4 in the centers of the faces. Each vertex site is shared by seven other unit cells and each face-center site by one additional unit cell. Hence, there are \(M^{\#}_{\text {C}}=4\) cation lattice sites belonging to one unit cell. There are \(M_{\text {A}}^{\#}=8\) anion lattice sites contained in the cation lattice unit cell, these are located inside the unit cell and not being shared by the neighboring unit cells. In general, the ratio \(m = M_{\text {A}}^{\#} / M_{\text {C}}^{\#}\) is a fixed constant that results from the given combination of materials. In the case of YSZ, we have m = 2. The spacing of the lattice can be described by a number density n# of unit crystal cells, that may be non-homogeneous in space due to the non-uniformity of the lattice. The densities of the cation lattice sites are then given as \(n_{\text {C}}^{\#} = M^{\#}_{\text {C}}n^{\#}\) while for the anion lattice sites is \(m M^{\#}_{\text {C}} n^{\#}\). We assume that all cation lattice sites are actually occupied by either zirconium or yttrium cations whereas some of the anion sites may be left unoccupied. We thus have

Illustration of the yttria-stabilized zirconia structure. The cation lattice is occupied by Zr4+ (blue) and Y3+ (orange). The anion lattice is occupied by mobile and immobile O2− (gray) or vacant (dashed)

To simplify the model, we assume the zirconium, yttrium and immobile oxide ions are bound to the lattice and thus all are transported with identical lattice velocity

Free energy and chemical potentials

The free energy density2ρψ of YSZ is assumed to be a function of temperature T, partial mass densities ρα and the electric field E. We suppose that the free energy density ρψ(T,ρα,E) can be split into four additive parts: reference energy, entropy of mixing, elastic energy and polarization energy,

Polarization energy

On top of the free charge density nF according to 2.2right, an excess charge density nP may arise in the material due to the presence of the electric field, mechanical strain, etc. (see for example [15, Chapter 2]). This excess charge is usually described by a polarization vector P so that

$$ -\text{div} \boldsymbol{P} = n^{\mathrm{P}} . $$

(2.12)

We refrain from a comprehensive discussion of constitutive modeling of polarization like, e.g., in [11] and assume that in bulk YSZ, the relaxation time of the polarization is small and the polarization vector P is proportional to the electric field E, i.e.,

The number χ is the electric susceptibility of YSZ, which for simplicity is assumed spatially homogeneous here. Integrating (2.13) such that ρψpolar vanishes for E → 0 yields the free energy density due to polarization

where \(\rho \tilde \psi = \rho \psi ^{\text {ref}}+\rho \psi ^{\text {mix}}+\rho \psi ^{\text {mech}}\). The elastic contribution to the free energy is based on a simple linear constitutive relation between the material pressure p and the number densities nα of YSZ,

Here K is the bulk modulus of YSZ and \(v_{\alpha }^{{\text {ref}}}\) are the specific volumes of the YSZ species under the reference pressure pref. In general, the specific volumes are functions of temperature and pressure, but for simplicity we assume \(v_{\alpha }^{{\text {ref}}}\) are constant.

By use of an alternative set of variables for the free energy density \(\rho \tilde \psi \) the Gibbs-Duhem relation (2.15) can be written as (cf. [6, equation A.6]),

Entropy of mixing

The entropy of mixing depends on the microscopic configuration of the mobile oxide ions in the anion lattice. We therefore consider a YSZ specimen that is homogeneous, so that nα = Nα/V, where Nα is the total number of a species in a volume V. Let W represent the number of possible realizations to arrange the mobile oxide ions on the anion lattice. Then, the mixing entropy density, according to Boltzmann’s formula, reads

Every immobile oxide ion is assumed to be fixed at a certain anion lattice site. The number of anion lattice sites available for the mobile oxide ions is therefore \((m N_{\mathrm {C}}^{{\#}} - N_{\text {Oi}})\). Thus, there are

Bulk governing equations and constitutive modeling

The electro-thermodynamic state of YSZ, occupying an interval \({\Omega }_{\text {YSZ}}\subset \mathbb {R}\) at any time t, is described by the number densities nα (\(\alpha \in \mathcal {I}_{\text {YSZ}}\)), the barycentric velocity υ and the electrostatic potential φ, which all are functions of time and position. In the isothermal electrostatic setting with a constant susceptibility, the evolution equations for the electro-thermodynamic state variables in the bulk are given by the Poisson equation, partial mass balances and the quasi-static momentum balance [11, 16],

To satisfy the second law of thermodynamics, i.e., to guarantee that the entropy production is non-negative, we insert the relations (2.25) into the entropy production (2.26) and then choose a linear relation between the diffusion flux JOm and the resulting term depending on the driving forces. We obtain

Here, mobility coefficient M may be a function of the thermodynamic variables and their derivatives, as long as it is guaranteed to be non-negative.

Incompressibility

A useful simplification of the YSZ bulk model is possible when taking the large bulk modulus K of YSZ into account. Hayashi et al. in [17] reported a bulk modulus of YSZ of K = 205 GPa at 25 ∘C and we assume it to be in a comparable order of magnitude at the operating temperature of YSZ at 600 ∘C. This motivates to study the incompressible limit \(\frac {K}{p^{\text {ref}}}\to \infty \). Under the assumption that the pressure p is bounded, we obtain from the constitutive relation (2.16) the constraint

Thus, the pressure p becomes an independent variable of the system and the sum of all number densities is independent of the pressure. For simplicity we assumed that the crystal lattice does not deform over time and that all species except species Om move with the lattice velocity. To be consistent with the incompressibility constraint (2.29), we thus have to require that the specific volume of the mobile oxide ions vanishes, i.e.,

Vanishing lattice velocity

For further simplification of the YSZ model, we assume that the lattice does not deform over time such that an appropriate reference frame can be chosen where the lattice velocity υ# vanishes,

$$ \begin{array}{@{}rcl@{}} {\upsilon}^{{\#}} = 0 . \end{array} $$

(2.32)

Then, the mass balance equations imply constant number densities for the immobile species, i.e., ∂tnα = 0 for α = Zr,Y,Oi, and the barycentric velocity is given by ρυ = ρOmυOm which can be expressed in terms of the diffusion flux of the mobile oxide ions as

The assumptions of incompressibility and vanishing lattice velocity may be also viewed alternatively as a description of the charge transport in the reference frame of the cation lattice which does not undergo any deformation.

Summary of the bulk YSZ model

The constitutive modeling above motivates to change the set of variables from the number densities \((n_{\alpha })_{\alpha \in \mathcal {I}}\) to \(\{ n_{\text {C}}^{\#},\nu ^{\#},x^{{\#}},y \}\). Due to the vanishing lattice velocity, the quantities \(n_{\text {C}}^{\#}\), x# and ν# are constant in time and are further considered as model parameters. Therefore, the thermodynamic state of the bulk YSZ is described by three quantities: filling ratio y, electrostatic potential φ and pressure p. In addition, we define the lattice volume V#, lattice mass m# and lattice charge number z# as

Let us assume that M is linearly dependent [18] on ρOm as \(M = D\frac {m_{\text {Om}}}{k_{\text {B}}}\rho _{\text {Om}}\). Eventually, the free charge and the diffusion flux of mobile oxide ions are given as

where (2.35c) was used in place of the pressure gradient term. The final form of the diffusion flux of the mobile oxide ion is proportional to the gradient of the electrochemcal potential.

The parameter x# has usually values in the range of [0,0.2] and we have \(\nu ^{\#} \in [0, \frac {1}{m}\frac {2+x^{{\#}}}{1+x^{{\#}}}]\). The remaining parameters of the YSZ model are given in Table 1.

Table 1

Characteristic values. Per-particle masses mα are used in the calculations

Temperature

T

800 ∘C

YSZ dielectric susceptibility

χ

27

Zr cation charge number

zZr

+ 4

Y cation charge number

zY

+ 3

Oxide ion charge number

zOm,zOi

− 2

Zr molar mass

MZr

91.22 g mol− 1

Y molar mass

MY

88.91 g mol− 1

O molar mass

MO

16 g mol− 1

Ratio of C/A lattices

m

2

YSZ molar fraction

x#

0.08

Ratio of immobile O2−

ν#

\([0, \frac {1}{m}\frac {2+x^{{\#}}}{1+x^{{\#}}}]\)

Specific lattice volume of YSZ

V#

3.35 × 29− 29m3

Lattice cation number density

\(n_{\text {C}}^{\#}\)

(V#)− 1

Diffusion coefficient

D

1× 10− 11 m2/s

Bulk metal and gas phase

In order to act as an electrolyte in a SOEC, the YSZ has to be connected to two different materials: a gas phase and some electric conductor. In this paper, we do no consider the internal structure of these parts of the SOEC. Therefore, we assume the gas to be equilibrated such that boundary conditions at the gas-YSZ surface can be determined easily. Although not appropriate for the use in real SOEC, we will treat the conductor as a pure metal, since this way the conductor can be almost completely removed from the model.

Bulk gas

The gas in the bulk is assumed to behave as an ideal mixture of ideal gases. We introduce the index set \(\mathcal {I}_{\text {gas}}\) of the constituents of the gas phase. For each constituent, the partial pressure is pα = cαRT. The chemical potential of a gaseous species reads

where the reference pressure is given by the standard atmospheric pressure pref = 100 kPa and \(\mu _{\alpha }^{\text {ref}}\) is the chemical potential of the pure substance.

In the bulk domain \({\mathrm {\Omega }}_{\text {gas}}\subset \mathbb {R}^{3}\), we assume that the diffusion is fast such that the chemical potentials are homogeneous in space, i.e., ∇μα = 0 for \(\alpha \in \mathcal {I}_{\text {gas}}\). Since there are no charge carriers in the gas, we assume that the electric potential φ is also homogeneous in the gas phase.

Bulk metal

For the description of the conductor, we apply the Sommerfeld model of metals (cf. [19]). The metal is considered as a mixture of positively charged metal ions M+ and free electrons e− with negligible volume and high mobility. Thus, we use the index set \(\mathcal {I}_{\text {metal}} = \{\text {M}^{+},\text {e}^{-}\}\) for the constituents. We assume the metal ions to be incompressible and thus the density of metal ions to be homogeneous in the whole metal domain \({\Omega }_{\text {metal}}\subset \mathbb {R}^{3}\) (cf. Landstorfer et al. [6]). Sufficiently far away from the metal boundary, i.e., outside of double layers, the metal is electroneutral and therefore the bulk number density \(n_{\text {e}^{-}}\) of the electrons and the corresponding bulk chemical potential \(\mu _{\text {e}^{-}}\) are material dependent constants. Neglecting electric resistance, the electric potential φ is homogeneous in the metal bulk. Moreover, we assume quasi-equilibrium in the metal such that in particular the electrochemical potential of the electrons is constant not only in the bulk but also inside double layers, i.e.,

Surface—triple phase boundary

The electrodes in solid oxide cells are combined of YSZ, metal and the gas phase. Thus, an interface model should, in principle, treat three thermodynamically distinct surfaces and one triple phase line present in the electrode. For a start, in this work, we aim at a strongly simplified 1D model of the electrodes. To incorporate the triple phase boundary into such a 1D model, we assume that the only contribution of the metal as an electric conductor is to provide free electrons for the charge transport. We make the following assumptions:

i) The YSZ surface is endowed with a thin layer of metal ions and their corresponding free electrons.

ii) The tangential transport of electrons along the surface is assumed to be fast compared to all the other treated kinetic processes.

iii) Apart from the reference free energy, the metal ions and electrons do not further contribute to the free energy and entropy of the surface.

Due to the first assumption, the electrons are transported only tangentially to the surface. The second assumption implies spatially homogeneous surface electrochemical potentials which only may change in time. The last assumption restricts the interaction between the metal and YSZ surface species. It allows to approximate the triple phase boundary by a simple surface model, which can be reduced to a 1D model in a straightforward way. There could be further contributions to the free energy, for example due to the surface tension of the metal ions, and/or due to an interaction of the metal species and the YSZ species. Since this would presumably add further parameters requiring fitting, we decided to treat the metal surface species, so that they do not introduce further interaction with the YSZ species.

A more detailed derivation of this reduction of a triple phase line into a 1D model can be found in the context of intercalation electrodes in [20].

The following derivation of the YSZ surface model is based on the general approach developed in [11, 16].

Surface constituents and basic quantities

As in the bulk, we describe the YSZ surface as a mixture of different surface constituents and apply for the surface quantities analogous notation with an underset “s” added. In the isothermal case, the surface temperature \(\underset {\text {s}}{T}\) is identical to the constant bulk temperature T and appears in the equations only as a parameter. In addition to the constituents from the metal and the bulk phases of the gas and YSZ bulk, surface reaction products may be present on the surface. Thus, the index set of all surface constituents is of the form \(\mathcal {I}_{S}=\mathcal {I}_{\text {YSZ}}\cup \mathcal {I}_{\text {gas}}\cup \mathcal {I}_{\text {metal}}\cup \mathcal {I}_{\text {react}}\), where \(\mathcal {I}_{\text {react}}\) is the index set of surface reaction products.

Each surface constituent is characterized by its surface number density \(\underset {\mathrm {s}}{n}{}_{\alpha }\), atomic mass mα and electric charge number zα. The partial mass densities \(\underset {\text {s}}{\rho }{}_{\alpha }\), the total mass density \(\underset {\text {s}}{\rho }\) and the free electric charge density for the surface are defined by

We assume that proper preparation and cutting of the bulk YSZ crystal results in the formation of a planar face which can be represented by our surface model. Therefore, as in the bulk YSZ case, the surface lattice density of cations is in certain relation to the surface density of anion lattice, i.e., the surface anion density is \(\underset {\text {s}}{m} \underset {\text {s}}{n}{{}_{\text {C}}^{\#}}\). The surface cation lattice is assumed to be fully occupied by zirconium and yttrium cations, whereas the anion lattice is partially occupied by mobile and immobile oxide ions (see Fig. 2).

The surface model needs to reflect the structure of the bulk YSZ model. The YSZ surface is defined by the cation crystal lattice. The deformation of the cation lattice therefore prescribes the surface velocity. In order to maintain the compatibility of the bulk model and the surface model, we have

On the YSZ surface gaseous species may adsorb and some reaction products may be formed. The admissible adsorption sites for gaseous species and reaction products in general depend on the lattice sites of the YSZ crystal. We assume that the density of the adsorption sites is proportional to the density of the anion surface lattice sites of YSZ. Several chemical reactions may occur. Denoting the constituents by Aα for \(\alpha \in \mathcal {I}_{S}\), the reactions can be written in the form

Surface free energy

The surface free energy can in general be assumed to be independent of the electric field. Here, we also assume that there is no elastic energy contribution and we distinguish two different entropic contributions to the free energy density. One takes into account the entropy of mixing of the mobile oxide ions on the anion lattice and the other is due to for the mixing of adsorbed gas species and reaction products on the adsorption sites. The metal ions and electrons only contribute to the reference energy. The free energy density for the surface is of the form

In general an elastic energy contribution has to be taken into account. The derivation of the energy is quite similar to the bulk. In [6] an example for a metal-electrolyte interface can be found. It turns out that if the constitutive equation of the surface tension depends only on the immobile YSZ species, and the lattice velocity υ# is equal to the surface velocity, then the remaining equations for the adsorption and surface reaction are independent of the elastic contribution. Therefore, for simplicity, we ignore the surface elasticity.

In analogous way, traces for functions in the gas bulk domain can be defined. Due to the choice of pairwise disjoint index sets for the bulk domains, most of the quantities are only defined in one of the subdomains. Therefore, we assume the simplification \(u|_{S} = \lim _{x\to S} u\). By convention, we let ν denote the outer normal of the YSZ domain.

In the planar one-dimensional approximation of the general surface mass balance equation (cf. [11, 16]), the tangential transport and the curvature related terms vanish. Only the surface chemical reactions (4.4) and mass transport normal to the surface can change the surface densities of the constituents. The surface mass balances and the remaining surface equation for the electric field in the electrostatic approximation read

We assume that Zr,YandOi are not involved in any surface reaction. Since \({\upsilon }_{\alpha } = {\upsilon }^{{\#}} = \underset {\mathrm {s}}{{\upsilon }}{}\) for α ∈{Zr,Y,Oi} according to (2.8) and (4.17), the surface mass balance (4.16a) implies that the corresponding surface number densities are constant, i.e.,

Maxwell’s surface equations in the electrostatic setting imply that the electrostatic potential is continuous at the gas-YSZ interface (see, e.g. [15] for further details). This allows us to introduce the surface electrostatic potential,

The entropy production can be split into the three contributions: \(\underset {\mathrm {s}}{\xi }{}_{\text {react}}\), \(\underset {\mathrm {s}}{\xi }{}_{\text {YSZ}}\) and \(\underset {\mathrm {s}}{\xi }{}_{\text {gas}}\), stemming from surface the reactions (4.4), adsorption from the bulk YSZ and adsorption from the gas phase, respectively. In analogous way like in the bulk, the structure of the entropy production (4.19) allows to derive constitutive equations such that the second law of thermodynamics is satisfied, i.e., the entropy production is non-negative.

Adsorption from YSZ bulk

Let us define the adsorption of oxide ions from the bulk to the surface as

where the second bracket on the right-hand side is equal to affinity of (4.20). By using a linear relation between the differences of chemical potentials and the mass flux, the entropy production \(\underset {s}{\xi }{}_{\text {YSZ}}\) is guaranteed to be non-negative,

Adsorption from gas phase

In the bulk gas phase, the fluxes are restricted by the constraint \({\sum }_{\alpha \in \mathcal {I}_{\text {gas}}} J_{\alpha } = 0\) and on the surface, (4.3) has to be satisfied. Therefore, we reformulate the entropy production due to the gas adsorption, as

with \({R_{0}^{i}} \ge 0\). The constants βi ∈ (0,1) are called symmetry factors.

Note that the arguments of the exponentials in (4.25) only depend on the surface-related quantities like surf. temperature and the chemical potentials of the surface species. Equation (4.25) is a surface generalization of the mass action kinetics. The adopted modeling approach implies the electric potential to be continuous, the overpotential5 appearing in the Butler-Volmer equation can be seen as an accounting for the potential drop due to the (unresolved) charge double layer. The presented model resolves the structure of the charged layer in the bulk YSZ and, in this sense, represents a more detailed description than the Butler-Volmer equation. In an asymptotic limit of vanishing double layer width the constitutive equation (4.25) allows to derive generalized Butler-Volmer equations for the surface reactions (see [21]).

Summary of the surface model

On the surface, we consider a single surface net reaction (5.2)right with β = 1/2. From the YSZ phase only the mobile oxide ions and from the conductor only the surface electrons are allowed to participate in this reaction. We assume fast adsorption from the gas phase, i.e., \(\mu _{\alpha }|_{S} = \underset {s}{\mu }{}_{\alpha }\) for \(\alpha \in \mathcal {I}_{\text {gas}}\).

Simulation of a SOC half-cell

We consider an YSZ-air electrode that contains the YSZ and gas bulk domains and the YSZ-gas surface S located at xS. We chose a point xB > xS outside of the double layer, located in the bulk YSZ sufficiently far away from S. Thus, the YSZ can be assumed to be electroneutral and consequently also isobaric in the YSZ bulk including xB. We assume that the pressure at xB corresponds to the outer pressure6 and the filling ratio of the anion lattice sites y at xB is determined by the crystal lattice, i.e.,

The adsorption of gaseous species is assumed to be considerably faster than the reaction and diffusion processes. Hence, the phenomenological coefficients in (4.24) for gaseous adsorbates are large, implying that the surface chemical potential and bulk chemical potential of the gas species are equal. Moreover, we assume fast dissociation, i.e., the reaction rate for the dissociation reaction (5.2) is large, and we obtain from (4.25)

Cell potential

The thin metal layer on the YSZ surface is assumed to be connected to a metal current collector, e.g., a wire. Therefore, there is an electric contact at the YSZ surface to an external circuit. Let \(\mu _{\text {e}^{-}}^{\text {ext}}\) and φext denote the (spatially homogeneous) chemical and the electrostatic potential in the current collector bulk, respectively. Assuming, that the electrochemical potential of the electrons is continuous at the surface, we can determine the contact potential \(U_{0}^{\text {ref}} = \varphi _{\text {ext}} -\underset {s}{\varphi }\) as

Due to the incompressibility of the metal bulk and the constitutive equation (4.14e) on the surface, the contact potential is a material dependent constant, i.e., \(\partial _{t} U_{0}^{\text {ref}} = 0\).

In principle, we are capable to measure the electrostatic potential φB at xB, e.g., with a suitable reference electrode. We define the half cell potential U of the solid oxide half-cell as

Thus, boundary condition for the electrostatic potential in the YSZ domain is given by the half cell potential U, and a normalization condition for φB, e.g., φB = 0.

Electric current

We are interested in the electric current I flowing through the electric wire contacted to the SOC electrode. The global mass balance equations allows us to relate the electric current I, flowing through the wire, to the quantities of the SOC electrode model as follows,

where A is the area of the cross section of the gas-YSZ interface. The derivation of formula (5.6) is summarized in the “Appendix:” section.

Double layer capacitance of blocking electrode

First, we want to investigate the equilibrium properties of the model derived above and therefore assume that no electron transfer reaction take place on the surface. This situation can be met if the contact of gas phase and YSZ is prevented by, e.g., metal layer. When an applied voltage is sustained so that the system is allowed to relax to an equilibrium state, and mobile oxide ions adsorb or desorb between the bulk and the surface and a charged layer in the bulk of YSZ is formed. We introduce the boundary layer charge QBL and the surface charge QS of the gas-YSZ interface as

respectively. Due to the 1D approximation, we are able to derive explicit representations of the bulk and surface capacitance as functions of the potential difference \(U-U_{0}^{\text {ref}}\). The homogeneity of the electrochemical potential in equilibrium, i.e.,

The impact of the mobility ratio ν# and of dielectric constant χ on the bulk layer capacitance is shown in Fig. 3. To screen a positive surface potential, a negatively charged layer in the YSZ has to be formed by occupying available anion lattice sites. Clearly, this number of available anion lattice sites is independent of the mobility ratio ν# and thus, the charge layer profile and the double layer capacitance CBL have to be independent of ν# for positive applied potentials. To the contrary, when the surface potential is more negative than the bulk, a small ν#, i.e., a large portion of the oxide anions is mobile, allows to vacate many anion lattice sides near the surface, leading to effective screening of the surface potential by a high negative charge density in the boundary layer and resulting in high capacity. The growth of the double layer capacitance for increasing χ, can be attributed to a spreading of the boundary layer due to the greater amount of the polarized charge.

Figure 4 shows the influence of ΔGA on the double layer capacitance of a blocking electrode. Negatively charged oxide ions tend to move into higher electric potential. If the adsorption energy, \({\Delta } G_{\mathrm {A}} = m_{\text {Om}}\underset {\mathrm {s}}{\mu }{}^{\text {ref}}_{\text {Om}} - m_{\text {Om}}\mu ^{\text {ref}}_{\text {Om}}\), is positive, then energy is required to for oxide ion to pass from the bulk to the surface. Stronger negative values of ΔGA foster the adsorption of oxide anions to the surface and thereby shift the surface capacitance maximum to more negative potentials. This can be seen most clearly in Fig. 4left where only the surface contribution is shown. Comparison of the Figs. 3 and 4 suggests that the bulk contribution remains undisturbed. The maxima of surface capacitance in Fig. 4left are due to the saturation of the surface for growing potential difference. The position of the maxima occurs for

Comparison to experiment

Figure 5 compares simulations with fitted data to experimentally measured capacitance curves for different temperatures [22]. For the simulation we used fitted parameters according to Table 3. The experimentally observed capacitances of the solid interface (see also [23] and [24]) are higher than for a metal electrode in aqueous solution. This is consistent with the absence of solvation shell formation in YSZ solid electrolyte.

We do not attempt to systematically adjust the model parameters to the data due to the polycrystalline nature of the YSZ studied in the experiment, instead, we try to illustrate the possible temperature dependence and the effect of the fitted parameters. As the temperature dependencies would need additional modeling efforts, as a first step, we performed the fit separately for each temperature.

It is difficult to assert that a particular oxide ion is mobile or immobile in the microscopic picture. It is suitable to consider the parameters ν# and \(\underset {\mathrm {s}}{\nu }{}^{{\#}}\) determining certain (dynamic) equilibrium between the admissible and occupied vacancies in state with vanishing macroscopic free charge density. As this is usually an effect of thermal excitations, the values of ν# and \(\underset {s}{\nu }^{{\#}}\) should depend on temperature. Also ΔGA presumably depends on the temperature.

To this end also \(\underset {\mathrm {s}}{m}{}\) was treated as a fitting parameter shared for the three cases.

Capacitive currents

In the case of time-dependent applied voltages, the current representation (5.6) simplifies in the absence of reactions to the case of no electron transfer:

Thus, the current is composed of two contributions describing the change of the surface charge and the boundary layer charge, respectively. However, unlike in the equilibrium case, QS and QBL are not uniquely determined by the applied voltage. Consider a small time depending perturbation around the half cell equilibrium potential \(\bar U\), i.e., the applied voltage is \(U(t) = \bar U + {\Delta } U(t)\). For a time scale of the perturbation considerably slower than the diffusion and adsorption, the system can be assumed to behave quasi-static and the current I can be linearized at \(\bar U\) such that

Thus, the double layer capacitance can be measured at low frequencies using impedance spectroscopy, or with cyclic voltammetry (CV) at low sweep rate. Here, sweep rate refers to the slope of voltage change during one linear cycle.

Kinetic coefficients

The blocking electrode model contains two kinetic parameters: diffusion coefficient D and adsorption rate \(\underset {\mathrm {s}}{A}{}_{0}\). If one of those parameters is small w.r.t. sweep rate, a limitation of the total current occurs. To this end the sweep rate is fixed to 1 mV s− 1 in this paragraph. Small values of adsorption coefficient limit charging and discharging of the surface oxide ions as it is shown in Fig. 6. The current due to charging of the bulk double layer is not affected by this.

Similarly, small values of D lead to limitation of the rate of charging the bulk double as documented in Fig. 7. In this case, the charging of the surface is affected, because the bulk diffusion limits the supply of the oxide ions.

Length of domain and sweep rate

Faster sweep rates affect the current response of the blocking electrode similarly as small values of the kinetic coefficients. Fast-changing voltage unveils limited rates of oxide ion transportation that can be attributed to concrete mechanisms. Figure 8 illustrates this for the oxide ion adsorption. Figure 8right in particular shows that the rate of the surface charging is limited due to the adsorption. For even greater sweep rates, the decreasing rates of current to the bulk diffusion limitation are displayed in Fig. 9left.

Cyclic voltammetry with realistic sweep rate rvolt = 1mV s− 1 is fixed in further demonstration of the basic features of the investigated system with the reaction.

Free energy parameters

Gibbs energy of reaction ΔGR is treated as an additional free energy parameter entering the model with the surface chemical reaction. The different ΔGR values (see Fig. 10left) do not alter the charging of the double layer but lead to the shift of the onset of the reaction current. Gibbs energy of adsorption ΔGA (see Fig. 10right) shifts, consistently with the blocking electrode case (cf. Fig. 4), charging of bulk a surface layer. The shift of the reaction onset occurs because ΔGA shifts the chemical potential of the surface oxide ions (cf. (5.21)). The reaction current is for either non-zero ΔGA or ΔGR in the depicted range much greater then the bulk and surface contributions.

Reaction rate

According to (5.6) surface reaction rate \(\underset {\mathrm {s}}{R}{}_{0}\) changes the relative magnitude of the reaction current; hence, it also changes the relative onset of the reaction current w.r.t bulk and surface contributions as it shown in Fig. 11. The limiting case of a small \(\underset {s}{R}{{}_{0}}\) is the blocking electrode.

The effects of D and \(\underset {\mathrm {s}}{A}{}_{0}\) are for the open system similar as for the blocking electrode case. Small values would lead to surface charging and consequently to bulk charging limitations thus hindering the reaction.

Discussion

The representation of the interface was chosen as uncomplicated as possible so that the behavior of oxide ions double layer dynamics remains unobscured. This was achieved, however, let us discuss the drawbacks of the treatment. First, in a real electrode two distinguished surfaces (YSZ, metal) are present and the electron-transfer reaction occurs near their intersection. Hence, tangential diffusion of the surface species comes into play together with the particular geometrical realization. To this end a two or three dimensional model would be required including the in-plane transport of the species. A question that naturally follows is: where exactly does the electron-transfer reaction occur, at the contact line or on one of the surfaces? Second, behavior of the metal electrons may in the close vicinity of the contact line start to display quantum effects that may result in richer behavior of the electron-transfer reaction. Third, the adsorption of gaseous species may under some circumstances limit the supply of gaseous species to the surface. Fourth, the appearing surface species depend on the particular electrode material. In particular, the nature and amount of the surface species will be different for Pt, Au or LSM electrodes. Also an additional phase of surface oxide ions with different adsorption energy might be present. Finally, one might consider production of surface oxygen O(s) for the blocking electrode (although no desorption to the gas phase is possible) and investigate the mechanical strain to the interfaces due to this.

Summary and Conclusions

A generalized Poisson-Nernst-Planck system describing YSZ|gas|metal-interface has been derived from first principles of nonequilibrium thermodynamics and numerically solved for simulating double layer capacitance and cyclic voltammetry measurements.

The core of the gPNP system is due to carefully derived free energy densities for the bulk YSZ and the YSZ|metal|gas surface capturing the main features of the YSZ crystalline nature. It is assumed that the described species, except for mobile oxide ions, are bound to the crystalline lattice. These assumptions result, using the entropy principle, in a novel form of the mobile oxide ion flux, which is a certain combination of the electrochemical potentials of all species. The charged layer in the metal is assumed to be in a diffusional equilibrium, since no transport limitations of the electrons is assumed. Finally, the formula for the electric current measured in the apparatus is derived.

A numerical model for the system has been derived and implemented in one spatial dimension using a finite volume method, specifically a variant of the Scharfetter-Gummel scheme, in the Julia programming language [25] for the details see Appendix A.

Although the model is strictly developed as isothermal, most of its parameters may depend on the temperature. Therefore, the parametric study is also aimed to demonstrate the scenarios where some of the parameters become limiting to the charge transfer of the system. Finally, the capacitance of blocking YSZ electrode taken from literature [22] is fitted with the model, the quality of the fit relies heavily on the newly introduced ratios of immobile oxide ions ν# and \(\underset {s}{\nu }^{{\#}}\). For each temperature these can be fitted alongside with ΔGA to the measured data. While the derivation of the model assumed a single crystal, the measurements had been obtained for polycrystalline YSZ. Therefore, the presented fitting results can be seen only as a first step towards a model for polycrystalline YSZ which ideally should be derived from the presented model using homogenization techniques. Moreover, the presented model can serve as a starting point for further extensions containing more sophisticated surface chemistry capable of describing the anodic and cathodic within one kinetic model.

We chose nitrogen as the reference species for the gas phase, i.e., A0 = N2.

Notes

Acknowledgments

The publication of this article was funded by the Open Access fund of the Weierstrass Institute.

Funding information

This work was supported by the German Research Foundation, DFG project no. FU 316/14-1, and by the Czech Science Foundation, GAČR project no. 19-14244J. VM received funding from the Fuel Cells and Hydrogen 2 Joint Undertaking under grant agreement no. 671481. This Joint Undertaking receives support from the European Union’s Horizon 2020 research and innovation programme, Hydrogen Europe and Hydrogen Europe research. PV and VM received partial support by grant SVV-2017-260455.

Publisher’s note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Appendix

Appendix A: Electric current

Let I be the electric current flowing through an electric wire to the gas-YSZ surface, which can be measured by an amperemeter. The current is related to the temporal change of the surface electron density and electron production on the gas-YSZ surface. For spatially homogeneous fields on a gas-YSZ surface with the area A, we have

The derivation is based on surface and bulk balance equations (see [11]), the metal model proposed in [6] and the assumption that the atomic mass of electrons are much smaller than the atomic mass of metal atoms, i.e., \(m_{\text {e}^{-}}/m_{M}\approx 0\).

To express the electron number density \(\underset {s}{n}{}_{\text {e}^{-}}\) in (A.1), we use the identity

Appendix C: The finite volume method

In order to perform the spatial discretization, we introduce collocation points \(x_{1}=x_{S}, x_{2}, \dots , x_{n-1},\)xn= xB in the simulation domain Ω = (xS,xB). The density of these points is increased in a geometric fashion towards the electrode surface at xS. Around the collocation points, we define the control volumes \(\omega _{1}=[x_{1},\frac {x_{2}+x_{1}}{2}]\), \(\omega _{i}=[\frac {x_{i}+x_{i-1}}{2},\frac {x_{i}+x_{i+1}}{2}]\) (\(i=2{\dots } n-1\)), \(\omega _{n}=[\frac {x_{n}+x_{n-1}}{2},x_{n}]\).

The finite volume discretization method used to perform the numerical simulation is based on the classical Scharfetter-Gummel scheme from semiconductor device simulation [26] which assumes constant species fluxes between neighboring control volumes, The fluxes are expressed via the unknowns in the corresponding collocation points based on an analytical solution of the flux equation. This approach automatically introduces an upwind stabilization of the discretization scheme which is necessary to handle the possibly steep electric potential gradients in the polarization boundary layer.

In order to handle the non-idealities occurring in generalized PNP models, the scheme needs to be adapted in a thermodynamically consistent manner. For an introductory discussion of the general ideas in the context of semiconductors (see 27), in the context of electrolyte simulation, a reformulation based on species activities as primary variables can be a starting point for a corresponding modification [28].

Here, we use an approach which starts from the reformulation of the species flux based on the introduction of a drift potential g(y,φ) combined of the excess chemical potential describing the non-ideality and the electrostatic potential, an idea which goes back at least to [29],

where yk,yl,φk and φl are values in computational nodes and \(B(x) := \frac {x}{\exp (x) - 1}\) is Bernoulli function. Under the assumption of jOm and \(\partial _{x} g(y, \varphi ) = g^{\prime }\) being constant, as in [26], the direct calculation of the numerical flux can be done using the integration factor \(\exp (-\frac {g^{\prime }}{\tilde {D}})\) and integrating over [xk,xl].

Copyright information

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.