A Model of Solid Oxide Fuel Cell Degradation on a Microstructural Level

Abstract

The growth of nickel (Ni) particles in the porous anode is one of the most critical issues in solid oxide fuel cells (SOFC). It reduces the density of triple-phase boundaries (TPBs) over time and increases the polarization resistance of SOFC. Most of the three-dimensional models that are used to simulate this phenomenon in detail are numerically exhausting and as such intractable for on-line applications. This work presents a two-dimensional, microstructural model of reduced complexity as a trade-off between the numerical load and the level of detail. The model of Ni agglomeration is based on the power-law coarsening theory. The resulting model was validated by comparing the relative density of TPBs and the cell voltage to the experimentally measured values. It was shown that the calculated values closely fit the measured data. The advantage of the proposed model is that it takes lower computational load during the simulation compared to the complex phase field models and is suitable for estimation of SOFC electric performance over time.

1. Introduction

Solid oxide fuel cell (SOFC) systems are a promising technology for stationary applications since they convert hydrogen (H₂) directly into electricity at high conversion efficiency (η) [1]. The η of chemical to electrical energy is usually higher compared to the gas turbines or engines with internal combustion, which are limited by the Carnot cycle [2]. SOFCs can be fueled with some other hydrocarbon gases (e.g., methane) that need to be reformed externally in a reformer [3] or internally within the cell [3,4]. When the cell is fueled with hydrocarbon gas, the exhaust gas also contains carbon dioxide (CO₂). The major drawbacks of using SOFC systems are relatively high manufacturing costs and a fast degradation rate. The overall degradation of SOFC is commonly a combination of electrochemical degradation, structural degradation, and mechanical failure [5]. One of the most critical issues is structural degradation of the porous anode due to nickel (Ni) agglomeration. Ni agglomeration increases the polarization resistance of SOFC since it decreases the density of electrochemical reaction sites within the anode active layer.

Numerous experimental and theoretical studies have been conducted to elucidate the degradation processes. Electric performance losses were studied when the SOFC was fed with fuel containing poisonous gas [6]. It was found that sulfur (S) adsorbs on the nickel (Ni) surface within the porous anode layer and inhibits the reaction sites at triple-phase boundaries (TPBs). Consequently, an obvious drop in SOFC output voltage was observed in the initial few hours of the test. The degradation mechanisms occurring within a planar, nickel-yttria-stabilized zirconia (Ni-YSZ), anode-supported, SOFC with direct internal reforming of dry CH₄-CO₂ mixtures were investigated during the ageing tests [7]. One of the tests was conducted for about 300 h at the temperature of 770 °C, current of 15 A (current density of 0.3 A cm⁻²), and the fuel utilization of 60%. The CH₄ and CO₂ flow rates were both 43 cm³ min⁻¹, whereas the airflow rate was 750 cm³ min⁻¹. The SOFC output voltage decreased by about 20% in less than 300 h due to the deactivation of Ni particles as a consequence of carbon formation and deposition [7].

Industrial size (10 cm × 10 cm), planar SOFCs were tested under various polarization conditions to show their degradation rate [8]. Post-mortem analysis of the microstructure revealed that the average size of Ni particles increased within the anode. Small Ni agglomeration occurred in the cell that operated under activation polarization (low current density), whereas large Ni agglomeration was observed in the cell operated under concentration polarization (high current density) [8]. This implies that the agglomeration of Ni particles has detrimental effects on electro-catalytic performance. The degradation reflects an increased polarization resistance of the anode that is attributed to the reduced density of TPBs within the anode due to the Ni grain coarsening [9].

Another paper [10] aims to predict the degradation of SOFC due to the growth of Ni particles. A theoretical approach to assess the agglomeration of Ni particles was adopted. A simple power law was employed to model the growth of average particle diameter over time. An excellent agreement with experimentally observed changes of the microstructure during the initial 1000 h of the test was accomplished. The change of TPB density was accounted for accordingly.

An improved phase field model was proposed recently [11] to simulate the microstructure evolution of the SOFC electrode. Three different sets of parameters for the YSZ/Ni anode, YSZ/lanthanum strontium manganite (LSM) cathode, and reference electrode were used to tune the interfacial energies in the model. The model was validated by simulating the grain growth and particle coarsening. The exponential growth of grain diameter with time was in good agreement with the simulation results presented in other studies. However, the electrical performances of presented SOFC microstructures were not calculated.

Many other studies deal with electric performances of SOFC, but the degradation models presented therein [12] cannot address the microstructural peculiarities of the anode and cathode layers, such as the spatial distribution of reaction sites, since these models are basically zero- or one-dimensional (0-D or 1-D). There is a paper [13] that addresses degradation of the SOFC due to different mechanisms, but it relies on very complex relations between the model parameters, such as length of TPBs and diameter of Ni grains, and microstructural properties. However, these relations are hard to prove when a limited set of measured data is available.

The majority of the approaches that rely on detailed 3-D models are numerically exhausting. In this work a trade-off between the numerical load and the level of precision of the degradation model was sought. Therefore, a reduced-complexity, 2-D model of Ni agglomeration in the anode layer was derived. The goal of this work was to study how the SOFC electric performance evolves with time. A simplified grain growth model was fitted by the density of TPBs that follows the degradation of microstructure. It was shown that the model is capable of predicting SOFC performance over a long period of time. Thus, the main contribution of this paper is a novel, two-dimensional (2-D) model that enables description of porous electrode properties on a microstructural level, but keeps numerical load at a reasonable level.

2. Methods and Modeling

The model takes into account species, mass, and charge conservation laws. The model is governed by partial differential equations that are described in Section 2.1, Section 2.2 and Section 2.3. Since the analytic solution of the system of equations could not be found, the equations were discretized with the finite difference method (FDM), adjusted according to the boundary conditions, as presented in Section 2.4, and solved numerically. A grey-scale image was used as the input of the model to define the cross-section profile of the SOFC, different phases (i.e., solid, electrolyte, and gas), and TPBs, as described in Section 2.5 and Section 2.6. The degradation model is presented in Section 2.7, whereas the input parameters are listed in Section 2.8. The nomenclature used throughout the paper can be found in Abbreviation and Appendix A.

2.1. Conservation of Species

Reduction Equation (1) and oxidation Equation (2) describe chemical reactions that occur at the TPBs within the porous cathode and anode:

$$\frac{1}{2}{\mathrm{O}}_2+{2\mathrm{e}}^{-}{=\mathrm{O}}^{2-}$$

(1)

$${\mathrm{H}}_2+{\mathrm{O}}^{2-}{=\mathrm{H}}_2{\mathrm{O}+2\mathrm{e}}^{-}$$

(2)

Since oxygen (O₂) in air and hydrogen (H₂) in fuel are consumed when SOFC operates as the electric source, their mass fractions (y_i) within the pores are different from the mass fractions in the gas channels, so the conservation equation of gas species is written in the following form:

$$\frac{\partial (\rho y_{\mathrm{i}})}{\partial t}+\overrightarrow{\nabla }\cdot (\rho y_{\mathrm{i}}\overrightarrow{v})+\overrightarrow{\nabla }\cdot \overrightarrow{j_{\mathrm{i}}}=S_{\mathrm{i}}$$

(3)

Gas density is denoted with ρ, y_i is mass fraction of the i-th gas species (index i or j denotes H₂, H₂O, O₂ or N₂), v is the average velocity of mass flux, j_i is the i-th component of the diffusive mass flux, S_i is source or sink of the i-th gas species. Diffusion of gas species is modeled by considering the effective diffusion coefficient D_i,j,eff. The j_i is obtained by the following equation:

$$\overrightarrow{j_{\mathrm{i}}}=-{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{N}-1}\rho \cdot D_{\mathrm{i},\mathrm{j},\mathrm{eff}}}\cdot \overrightarrow{\nabla }{y}_{\mathrm{j}}=-{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{N}-1}\rho \cdot {\mu }^{\frac{3}{2}}\cdot D_{\mathrm{i},\mathrm{j}}}\cdot \overrightarrow{\nabla }{y}_{\mathrm{j}}$$

(4)

N denotes the number of gas species in gas mixture (N = 2 since only two gas species are present). The diffusion coefficient D_i,j of a binary mixture of gas species i and j can be estimated as follows [14]:

$$D_{\mathrm{i},\mathrm{j}}=\frac{{10}^{-7}\cdot T^{\frac{7}{4}}}{p\cdot M_{\mathrm{i},\mathrm{j}}^{\frac{1}{2}}\cdot {\left(V_{\mathrm{i}}^{\frac{1}{3}}+V_{\mathrm{j}}^{\frac{1}{3}}\right)}^2}$$

(5)

$$M_{\mathrm{i},\mathrm{j}}={\left(\frac{1}{M_{\mathrm{i}}}+\frac{1}{M_{\mathrm{j}}}\right)}^{-1}$$

(6)

T is the temperature in Kelvin (K), p is the pressure in bar (1 bar = 10⁵ Pa), M_i/j is the molar mass (g mol⁻¹) of gas species in the binary mixture, and V_i/j is the diffusion volume in cm³. Detailed values of M_i and V_i can be found in one of our previous papers [15]. The S_i is calculated at the TPBs, where all three phases (i.e., solid, electrolyte, and gas) meet, otherwise it is zero. The fuel (H₂) is consumed and water (H₂O) is produced at the TPBs in the anode as the SOFC operates, so the sink/source term (S_H2/S_H2O) depends on anode current (i_a). Faraday constant is denoted with F = 96,485 As mol⁻¹.

$$S_{{\mathrm{H}}_2}=-\frac{M_{{\mathrm{H}}_2}\cdot i_{\mathrm{a}}}{2\cdot F}$$

(7)

$$S_{{\mathrm{H}}_2\mathrm{O}}=+\frac{M_{{\mathrm{H}}_2\mathrm{O}}\cdot i_{\mathrm{a}}}{2\cdot F}$$

(8)

The oxygen (O₂) is consumed at the TPBs in the cathode as the SOFC operates, so the sink term (S_O2) depends on cathode current (i_c):

$$S_{{\mathrm{O}}_2}=-\frac{M_{{\mathrm{O}}_2}\cdot i_{\mathrm{c}}}{2\cdot F}$$

(9)

It should be noted that nitrogen (N₂) is neither consumed nor produced, but its mass fraction depends on the mass fraction of oxygen, which changes within the pores as the SOFC operates at different output current. The mass fraction of water (y_H2O) is calculated as 1 − y_H2, whereas the mass fraction of nitrogen (y_N2) is calculated as 1 − y_O2.

2.2. Conservation of Mass

The mass conservation equation is written in the following form:

$$\frac{\partial \rho }{\partial t}+\overrightarrow{\nabla }\cdot (\rho \overrightarrow{v})={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{N}}{S}_{\mathrm{i}}}$$

(10)

The density of a gas is denoted with ρ, v is the average velocity of mass flux, S_i represents sink or source of gas species, N is the number of gas species in the air or fuel. Since an ideal gas is assumed, the density of air (ρ_air) and fuel (ρ_fuel) depends on temperature (T), total pressure (p), and mass fraction of gas species (y_i). Ideal gas constant is denoted with R = 8.314 J mol⁻¹ K⁻¹.

$${\rho }_{\text{air}}=\frac{p\cdot M_{{\mathrm{O}}_2}\cdot M_{{\mathrm{N}}_2}}{R\cdot T\cdot \left(y_{{\mathrm{O}}_2}\cdot M_{{\mathrm{N}}_2}+(1-y_{{\mathrm{O}}_2})\cdot M_{{\mathrm{O}}_2}\right)}$$

(11)

$${\rho }_{\text{fuel}}=\frac{p\cdot M_{{\mathrm{H}}_2}\cdot M_{{\mathrm{H}}_2\mathrm{O}}}{R\cdot T\cdot \left(y_{{\mathrm{H}}_2}\cdot M_{{\mathrm{H}}_2\mathrm{O}}+(1-y_{{\mathrm{H}}_2})\cdot M_{{\mathrm{H}}_2}\right)}$$

(12)

It should be noted that discretization of the divergence terms in Equations (3) and (10) should be done carefully, since ρ_air/fuel changes within pores, and thus across each discretization step. If the central differencing scheme is used, ρ_air/fuel should be evaluated in the middle of the two adjacent discretization points. A good approximation is obtained if the average mass fraction is used in Equations (11) and (12). Figure 1 shows the aforementioned situation schematically.

2.3. Conservation of Charge

The charge transfer through the modeled SOFC is described by different charge exchange mechanisms that are modeled in the following subsections. A schematic diagram is shown in Figure 2 to illustrate these mechanisms. Each number (1–5) in Figure 2 depicts the mechanism that is described in the equally numbered subsection.

2.3.1. Surface Exchange

It is assumed the oxygen (O₂) that fills the pores of solid (i.e., mixed ionic-electronic conductive—MIEC)) material is adsorbed into the cathode surface. Oxygen ions (O²⁻) are formed with a double negative charge by accepting electrons (e⁻) from the cathode. The flow of O²⁻ (J_O²⁻) can be modeled as follows [16]:

$$\overrightarrow{n}\cdot {\overrightarrow{J}}_{{\mathrm{O}}^{2-}}=-k_{\text{surf}}\cdot \left(C_{{\mathrm{O}}^{2-},\mathrm{eq}}-C_{{\mathrm{O}}^{2-}}\right)$$

(13)

The surface exchange coefficient is denoted with k_surf. It depends on the temperature and the partial pressure of O₂. The equilibrium oxygen ion concentration is denoted with C_{O²⁻, eq} and depends on the partial pressure of O₂ within pores. The concentration of oxygen ions (C_O²⁻) decreases during the increasing electric current load (C_O²⁻ < C_{O²⁻, eq}).

2.3.2. Bulk Diffusion

Since the electronic conductivity of the cathode MIEC material is high and should not limit the electric current, the diffusion of O²⁻ is assumed to limit the current. The diffusion of O²⁻ is modeled by using Fick’s first law:

$${\overrightarrow{J}}_{{\mathrm{O}}^{2-},\text{bulk}}=-D_{\text{chem}}\cdot \overrightarrow{\nabla }{C}_{{\mathrm{O}}^{2-}}$$

(14)

D_chem is chemical diffusion coefficient that depends on temperature and partial pressure of O₂ [17].

2.3.3. Charge Transfer

At the interface between the solid MIEC and electrolyte material the O²⁻ are exchanged. The local charge transfer voltage (overpotential) is the driving force for this process [17]. The current density is modeled by the Butler–Volmer equation:

$$i_{\mathrm{a}/\mathrm{c}}=\{\begin{array}{c}{i}_{0,\mathrm{a}}\cdot \left[\text{exp}\left(\frac{{\mathsf{\alpha }}_{\mathrm{f}}^{\mathrm{a}}\cdot F\cdot (V_{\mathrm{s}}-V_{\text{el}}-V_{0,\mathrm{a}})}{R\cdot T}\right)-\text{exp}\left(\frac{{\mathsf{\alpha }}_{\mathrm{b}}^{\mathrm{a}}\cdot F\cdot (V_{\mathrm{s}}-V_{\text{el}}-V_{0,\mathrm{a}})}{R\cdot T}\right)\right];\text{ at TPBs in the anode}.\\ i_{0,\mathrm{c}}\cdot \left[\text{exp}\left(\frac{{\mathsf{\alpha }}_{\mathrm{f}}^{\mathrm{c}}\cdot F\cdot (V_{\mathrm{s}}-V_{\text{el}}-V_{0,\mathrm{c}})}{R\cdot T}\right)-\text{exp}\left(\frac{{\mathsf{\alpha }}_{\mathrm{b}}^{\mathrm{c}}\cdot F\cdot (V_{\mathrm{s}}-V_{\text{el}}-V_{0,\mathrm{c}})}{R\cdot T}\right)\right];\text{ at TPBs in the cathode}.\end{array}$$

(15)

The anode/cathode exchange current density is denoted with i_0,a/c, forward/backward exponential coefficient of anode/cathode with ${\mathsf{\alpha }}_{\mathrm{f}/\mathrm{b}}^{\mathrm{a}/\mathrm{c}}$ (for simplicity, ${\mathsf{\alpha }}_{\mathrm{f}}^{\mathrm{a}}=0.5$, ${\mathsf{\alpha }}_{\mathrm{f}}^{\mathrm{c}}=0.5$, ${\mathsf{\alpha }}_{\mathrm{b}}^{\mathrm{a}}=-0.5$, ${\mathsf{\alpha }}_{\mathrm{b}}^{\mathrm{c}}=-0.5$ are used), V_s is solid potential, V_el is electrolyte potential, and V_0,a/c is anode/cathode potential that can be obtained from Nernst equation, respectively.

$$V_{0,\mathrm{a}}=0$$

(16)

$$\begin{array}{l}{V}_{0,\mathrm{c}}=\frac{-\Delta G_{{\mathrm{H}}_2}-k_{\text{th}}\cdot T-R\cdot T\cdot \text{ln}\left(\frac{p_{{\mathrm{H}}_2}\cdot \sqrt{p_{{\mathrm{O}}_2}}}{p_{{\mathrm{H}}_2\mathrm{O}}}\right)}{2\cdot F}=\\ =\frac{-\Delta G_{{\mathrm{H}}_2}-k_{\text{th}}\cdot T-R\cdot T\cdot \text{ln}\left(\frac{y_{{\mathrm{H}}_2}\cdot M_{{\mathrm{H}}_2\mathrm{O}}}{(1-y_{{\mathrm{H}}_2})\cdot M_{{\mathrm{H}}_2}}\cdot \sqrt{\frac{y_{{\mathrm{O}}_2}\cdot M_{{\mathrm{N}}_2}}{y_{{\mathrm{O}}_2}\cdot M_{{\mathrm{N}}_2}+(1-y_{{\mathrm{O}}_2})\cdot M_{{\mathrm{O}}_2}}}\right)}{2\cdot F}\end{array}$$

(17)

Partial pressures of oxygen, hydrogen and water are denoted with p_O2, p_H2, and p_H2O, respectively, and the ΔG_H2 denotes Gibbs free energy for hydrogen. The temperature coefficient k_th is calculated from Equation (17), considering open circuit voltage of measured SOFC (OCV = V_0,c − V_0,a), operating temperature, and partial pressures (or mass fractions) of gas species.

2.3.4. Electronic Current

Electronic current is modeled in solid MIEC material of the cathode by the continuity equation [17]:

$$\overrightarrow{\nabla }\cdot \left(-{\sigma }_{\mathrm{s},\mathrm{c}}\cdot \overrightarrow{\nabla }{V}_{\mathrm{s}}\right)+2\cdot F\cdot {\delta }_{\mathrm{ads}}(\overrightarrow{r})\cdot k_{\text{surf}}\cdot \left(C_{{\mathrm{O}}^{2-},\text{eq}}-C_{{\mathrm{O}}^{2-}}\right)=0$$

(18)

The electronic conductivity of cathode MIEC material is denoted with σ_s,c, δ_ads denotes delta function:

$${\delta }_{\text{ads}}(\overrightarrow{r})=\{\begin{array}{c}1;\text{ at MIEC}/\text{gas interface in cathode}\\ 0;\text{ elsewhere}\end{array}$$

(19)

Similarly, the electronic current is modeled in solid MIEC material of anode:

$$\overrightarrow{\nabla }\cdot \left(-{\sigma }_{\mathrm{s},\mathrm{a}}\cdot \overrightarrow{\nabla }{V}_{\mathrm{s}}\right)+{\delta }_{\mathrm{a}}(\overrightarrow{r})\cdot i_{\mathrm{a}}=0,{\delta }_{\mathrm{a}}(\overrightarrow{r})=\{\begin{array}{c}1;\text{ at TPB in anode}\\ 0;\text{ elsewhere}\end{array}$$

(20)

The electronic conductivity of anode MIEC material is denoted with σ_s,a, δ_a denotes delta function.

2.3.5. Ionic Current

The ionic current in anode/cathode/electrolyte is modeled with the continuity equation [17]:

$$\overrightarrow{\nabla }\cdot \left(-{\sigma }_{\text{ion},\mathrm{a}/\mathrm{c}/\text{el}}\cdot \overrightarrow{\nabla }{V}_{\text{el}}\right)-{\delta }_{\mathrm{a}/\mathrm{c}}(\overrightarrow{r})\cdot i_{\mathrm{a}/\mathrm{c}}=0,{\delta }_{\mathrm{a}/\mathrm{c}}(\overrightarrow{r})=\{\begin{array}{c}1;\text{ at TPB in anode}/\text{cathode}\\ 0;\text{ elsewhere}\end{array}$$

(21)

The ionic conductivity of anode/cathode/electrolyte is denoted with σ_ion,a/c/el, δ_a/c denotes delta function.

2.4. Boundary Conditions

Boundary conditions (BCs) have to be specified at different interfaces of the modeled SOFC structure in order to obtain a closed solution of the system of coupled partial differential equations. Figure 3 shows the interfaces A–K within the modeled domain, whereas a description of BCs at these interfaces is given in

Table 1. The boundary conditions.

Interfaces/Variables	Vs (V)	Vel (V)	yi (/)	v (m s−1)	C (mol m−3)
A	$V_{\mathrm{s}}=V_{0,\mathrm{c}}-V_{\mathrm{b}}$	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\text{el}}=0$	/	/	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{{\mathrm{O}}^{2-}}=0$
B	/	/	$y_{\mathrm{i}}=y_{{0,\mathrm{O}}_2}$	$\partial v_y/\partial y=0$	/
C	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\mathrm{s}}=0$	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\text{el}}=0$	/	/	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{{\mathrm{O}}^{2-}}=0$
D	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\mathrm{s}}=0$	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\text{el}}=0$	$\overrightarrow{n}\cdot {\overrightarrow{j}}_{\mathrm{i}}=0$	$v_x=0$	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{{\mathrm{O}}^{2-}}=0$
E	$V_{\mathrm{s}}=V_{0,\mathrm{a}}$	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\text{el}}=0$	/	/	/
F	/	/	$y_{\mathrm{i}}=y_{{0,\mathrm{H}}_2}$	$\partial v_y/\partial y=0$	/
G	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\mathrm{s}}=0$	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\text{el}}=0$	/	/	/
H	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\mathrm{s}}=0$	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\text{el}}=0$	$\overrightarrow{n}\cdot {\overrightarrow{j}}_{\mathrm{i}}=0$	$v_x=0$	/
I	/	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\text{el}}=0$	/	/	/
J	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\mathrm{s}}=0$	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\text{el}}=i_{\mathrm{a}/\mathrm{c}}$	/	/	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{{\mathrm{O}}^{2-}}=0$
K	/	$\overrightarrow{n}\cdot {\overrightarrow{J}}_{\text{el}}=0$	$\overrightarrow{n}\cdot {\overrightarrow{j}}_{\mathrm{i}}=0$	$v_y=0$	/

.

2.5. Modeled Structure

An image of a SOFC cross-section profile was used to define the geometry of modeled domains (cathode, electrolyte, and anode) within the SOFC. The image had to be processed properly to detect the three different phases (solid, electrolyte, and gas) present within the real SOFC. Different image processing algorithms can be found in the literature [18,19,20]. Many ideas of how to choose a global threshold have been already presented [19]. Most of these methods derive the threshold value using the properties of an image histogram. A similar problem of how to detect different phases from threshold values has been already addressed in Ref. [20]. It is assumed that the grey values of the phases are normally distributed and sufficiently apart from each other. Hence their superposition results in a local minimum between the two phases when inspected in the pixel histogram [20].

A simple algorithm for detecting the three different phases was developed in our case. First, a vector with 256 locations was created. Second, the pixels in the image that have a certain grey scale value (GSV) were counted and the sum was stored at the GSV-th location in the vector. This procedure was repeated for each GSV ∈ [0,1,…,255]. Third, three local maximums of the vector were sought. Each GSV at which a local maximum was found corresponded to one of the three different phases, considering these phases had normal distribution of the pixel number with respect to this GSV. Fourth, two local minimums were sought by considering each one as located between two maximums. Fifth, two threshold values were obtained from indexes of the locations in the vector, at which two local minimums were found. An example of detection of two threshold values that differentiate the three different phases from the grey-scale image is shown in Figure 4.

2.6. Solid-Electrolyte-Gas Interfaces

Three different phases (solid, electrolyte, and gas phases) were initially set to zero values at each discretized location (w, z) of the modeled SOFC. These phases had to be unambiguously determined based on two threshold values, which were resolved as described in Section 2.5. The following criteria were considered:

• if (GSV < Threshold 1) then Gas (w, z) = 1;
• else if ((GSV ≥ Threshold 1) and (GSV < Threshold 2)) then Electrolyte (w, z) = 1;
• else Solid (w, z) = 1;

In the next step, only the connected (percolated) solid and electrolyte phases were determined in between the electric contacts, which were denoted as interfaces A and E in Figure 3. The connected gas phases to interface B and F were determined. For example, delta function δ_a/c (w, z) = 1 since TPB was active within the porous anode or cathode if all three phases ere percolated and met in a single point (w, z). Otherwise, δ_a/c (w, z) = 0 since TPB was inactive. This situation is shown schematically in Figure 5a. The upper part of Figure 5b shows active TPBs within the modeled SOFC. The lower part of Figure 5b shows the actual 2-D profile of the mass fraction (Y) of O₂ and H₂ that corresponds to the percolated gas phase within the porous cathode and anode, respectively, as seen in the upper part. Please note that TPB detection actually needs four neighboring pixels (i.e., a small image of 2 × 2 pixels) that represent different phases. It is obvious that TPB detection strongly depends on both threshold levels. Consequently, there is some uncertainty in this kind of detection, as it is impossible to determine whether a TPB at the interface really exists or not. Moreover, the reconstruction of porous material in this paper was based on the 2-D grey-scale image, which may introduce additional uncertainty. Detailed analysis and detection of TPB would require three-dimensional (3-D) tomography techniques, e.g., focused ion beam-scanning electronic microscopy (FIB-SEM) [21] and X-ray computed tomography (CT) [22]. Those provide more accurate SOFC microstructure characterization. These cutting-edge 3-D characterization techniques are very expensive and time-consuming [23]. This is why the detection in our case was restricted to 2-D grey-scale images only.

2.7. The Degradation Model

Degradation mechanisms of SOFC have been studied experimentally [24,25] and theoretically [11] to investigate how the microstructure evolves over long-term operation under different conditions, such as temperature, current density [24,25], and initial composition of the porous anode, and the average size of Ni and YSZ grains within the microstructure [11].

Since a real SOFC operates at high temperatures (at T = 800–1000 °C), the microstructure evolution, such as grain growth or particle coarsening, is significant [11]. The microstructure becomes increasingly coarser and denser with increasing testing temperature, current density, and testing time [24]. At 750 °C an increase in grain size occurs mainly in the bulk and interface compared with a non-degraded cell. At 850 °C and above, the grain coarsening in the bulk becomes more pronounced than at 750 °C. The coarsened interface region is about 1–2 μm in width at 850 °C and up to 5 μm at 950 °C. It was also seen that at 850 °C the process of grain coarsening and interface densification continues with increasing testing time up to 17,500 h [24]. From these results one can note that microstructure evolution generally progresses with increasing temperature and testing time. Ostwald-ripening is an important mechanism of Ni coarsening in SOFC anode [26]. It occurs due to the transport of volatile Ni species via the gas phase and diffusion of vacancies, driven by different grain sizes.

Despite the detailed studies, only a few papers [27,28] deal with the problem of how the microstructural degradation actually affects the electric performance of SOFC. In this paper the electric performance of a degraded SOFC was calculated by adjusting fewer parameters compared to the number of parameters that need to be set in the phase field model [11]. Moreover, the validation of such a model would require detailed analysis with three-dimensional (3-D) tomography techniques. These techniques have shown that the Ni (YSZ) grains grow exponentially with time [9]. Therefore, the degradation of SOFC due to grain growth was modeled by considering the power law [9]:

$$d^{\kappa }=d_0{}^{\kappa }+k\cdot t$$

(22)

Time-dependent diameter of Ni grains is denoted with d, initial diameter with d₀, growth coefficient with k and time with t. The exponential coefficient is denoted with κ. It equals 2 for grain growth and 3 for particle coarsening due to Ostwald-ripening [11].

However, Equation (22) has to be slightly adapted for numerical implementation. It is assumed that d₀ is determined by the minimum distance between the two adjacent points where active TPBs are located. This simplification is done due to problems with detecting grain boundaries of the Ni (YSZ) particles and their actual diameter by processing a grey-scale image. Moreover, reconstruction of the porous anode according to the phase field model would be necessary if actual grain growth was considered rigorously [11]. That is beyond the scope of the present paper, so the following algorithm was proposed.

First, the active TPBs are resolved as described in Section 2.6. Second, two vectors are defined for storing two distances (d_1,min, d_2,min) from one to the closest two TPB points. Third, d_1,min and d_2,min are found. Fourth, the SOFC model is evaluated at each iteration (q) using a fixed time step (Δt). Fifth, the following condition is checked at each iteration and at each active TPB point:

$$\text{if }\left({\left(d_{1,\text{min}}\left(\mathrm{w},\mathrm{z}\right)\left[\mathrm{q}\right]+d_{2,\text{min}}\left(\mathrm{w},\mathrm{z}\right)\left[\mathrm{q}\right]\right)}^{\kappa }<{\left(d_{1,\text{min}}\left(\mathrm{w},\mathrm{z}\right)\left[\mathrm{q}\right]\right)}^{\kappa }+\mathrm{q}\cdot k\cdot \Delta t\right)\text{ then }{\delta }_{\mathrm{a}}\left(\mathrm{w},\mathrm{z}\right)=0$$

(23)

This algorithm runs iteratively as long as it is defined by a maximum number of iterations q. It should be noted that after a certain number of iterations, the condition in Equation (23) is fulfilled at some active TPB point and this TPB is deactivated (δ_a (w, z) = 0). In this way, the particle diameter is increased. It is important to note the number of active TPBs is decreased, so this model enables fitting the parameters k and n to achieve good approximation to results of measured or calculated relative TPB density from the literature [10,11].

2.8. Model Parameters

Detailed descriptions of measured SOFC, including the lanthanum strontium cobalt ferrite (LSCF) cathode, Ni-scandia stabilized zirconia (Ni-SSZ) active (or support) anode layer, and electrolyte bilayer, can be found in Ref. [10]. The physical parameters are adopted from the literature, or fitted to match the desired characteristics of measured SOFC.

Table 2. Input parameters of the SOFC model.

Parameter	Symbol	Value	Unit	Reference
Pressure	p	1.0	bar	[10]
Anode thickness	da	20	μm	[10]
Cathode thickness	dc	20	μm	[10]
Electrolyte thickness	de	10	μm	[10]
Anode electronic conductivity	σs,a	4800	S m−1	[15]
Cathode electronic conductivity	σs,c	1600	S m−1	[15]
Anode ionic conductivity	σel,a	2.0	S m−1	[29]
Cathode ionic conductivity	σel,c	2.0	S m−1	[29]
Electrolyte ionic conductivity	σel,e	2.0	S m−1	[29]
Anode exchange current density	i0,a	1200	A m−2	[10]
Cathode exchange current density	i0,c	2200	A m−2	[10]
Surface exchange coefficient	ksurf	2.0 × 10−5	m s−1	[17]
Concentration of oxygen ions	CO2-	8.3 × 104	mol m−3	[17]
Chemical diffusion coefficient	Dchem	4.4 × 10−10	m2 s−1	[17]
Series resistance of anode support layer	Rs,a	0.5	Ω cm2	[10]
Double layer capacitance of anode	Cdl,a	0.5	F cm−2	[10]
Molar fraction of hydrogen in fuel	x0,H2	0.970	/	[10]
Molar fraction of water in fuel	x0,H2O	0.030	/	[10]
Molar fraction of oxygen in air	x0,O2	0.207	/	[10]
Molar fraction of nitrogen in air	x0,N2	0.793	/	[10]

shows the list of input parameters.

3. Results and Discussion

The SOFC electric performance is simulated by using the model according the parameters in Table 2. We start with AC analysis because it is intended to adjust the anode/cathode exchange current density (i_0,a/c) and qualitatively validate the simulation results with measured electro-impedance spectrum (EIS) in Ref. [10]. Note that the adjustment of i_0,a/c should be done initially for online application of the model in practice, by fitting the calculated and measured EIS. Exact fitting of the EIS was not of the primary concern in this study case, but rather to achieve similar trends of simulated and measured EIS curves.

3.1. AC Analysis

The modeled cell was perturbed at OCV, according to the measurement results in Ref. [10], with sinusoidal voltage with the amplitude 14 mV. The frequency (f) of voltage signal was set to each of 31 logarithmically equally spaced values of frequency from 0.1 Hz to 100 Hz to simulate the EIS. The temperature was set to 900 °C. Other input parameters of the model can be found in Table 2. Figure 6 shows EIS of (a) non-degraded (i.e., at time t = 0 h) and (b) degraded SOFC (i.e., at time t = 1000 h). The low-range frequency arc was attributed to the anode charge transfer, whereas the mid-range frequency arc was correlated with charge transfer and surface exchange process in the cathode.

The anode and cathode polarization resistances increase over time (i.e., from 0 to 1000 h) due to the reduced density of TPBs within the anode and cathode. The increased polarization resistance of the degraded SOFC indicated a correct trend regarding the time evolution of measured EIS, as can also be seen in Figure 8c in Ref. [10]. The measured series resistances were higher than in this case, which may be attributed to the wiring and contact resistances, which were not considered in the presented model. Moreover, if considering only the series resistances (about 1.05 Ω cm² at non-degraded and 1.4 Ω cm² at the degraded cell) that can be deduced from measured EIS, the voltage drop of SOFC would be about 0.47 V and 0.63 V at the current density (J) of 0.45 A cm⁻², so the direct current (DC) voltage of the SOFC (V_cell) should be lower than shown in Figure 8a in Ref. [10]. The measured V_cell of non-degraded SOFC was about 0.56 V (at time t = 0 h), whereas V_cell of degraded SOFC was about 0.42 V (at time t = 1000 h).

Since the open circuit voltage (OCV) of SOFC is about 0.97 V, the measured V_cell should have been 0.50 V and 0.34 V when the aforementioned voltage drops were considered, respectively. Due to these discrepancies, the model was fitted only to reproduce the degradation of V_cell, as presented in Section 3.2, since we focus on calculating the electric performance of SOFC over time. However, a profound inspection of the measurement equipment and physical parameters of the real SOFC device would be necessary to obtain a better fitting of EIS, whereas a comparison of impedances could be done only if detailed measurement data were available in Ref. [10].

3.2. DC Analysis

In the second set of experiments, the DC analysis was performed. A comparison between calculated and measured electric performance of SOFC was done to address more realistic degradation due to the coarsening of Ni particles within the anode. A fraction of backscattered electron (BSE) microscope pictures were taken from Figure 3a in Ref. [10] to define an approximation to the real structure of porous anode. The experimental degradation test was performed for 1000 h at the temperature of 900 °C and at current density of 0.45 A cm⁻², whereas 3 vol% humidified H₂ was supplied to the anode and air was supplied on the cathode side. The flow rates of hydrogen and air were 350 cm³ min⁻¹ and 1500 cm³ min⁻¹, respectively [10]. Due to this, the temperature and current density of the modeled SOFC was set to T_cell = 900 °C and J_cell = 0.45 A cm⁻², respectively, during the simulation from 0 to 1000 h. Time increment Δt = 100 h was used since the experimental data was defined at 100 h intervals in Ref. [10]. The anode exponential coefficient (κ_a = 3) and the anode growth coefficient (k_a = 2.1 × 10⁻²¹ m³ h⁻¹) were manually fitted to match the relative TPB density of the measured SOFC in Ref. [10].

Figure 7a shows relative TPB density in the anode (TPB_a) and in the cathode (TPB_c) as a function of time. The calculated relative TPB_a density had a similar profile as that in Figure 8d in Ref. [10]. The initial absolute value of TPB_a (t = 0 h) was approximately four times lower (0.7 μm⁻² vs. 2.9 μm⁻²), but another study revealed [30] that the density of active TPBs may be lower (0.74 μm⁻²). Since the TPBs that lie in the third dimension (i.e., in the plane of 2-D grey-scale image) cannot be resolved with our algorithm, it was obvious why the obtained density of TPBs from 2-D grey-scale image was lower. Although the TPB density was lower in our case than it should have been (or than it was in the real SOFC), it could be compensated for with higher activity of TPBs (i.e., higher exchange current density, i_0,a/c) to achieve a similar electric performance of modeled SOFC. It is important to note the presented model approximated the spatial distribution of TPBs in real SOFC. Due to this, it could be expected that the simulation results may be more reliable compared to those obtained with 0-D or 1-D models.

Figure 7b shows the calculated voltage (V_da) of SOFC, with degraded anode, which is almost the same as the reference voltage (V_Ref.[10] in Figure 8a in Ref. [10]) from 0 to 100 h, whereas a higher voltage difference (0.03–0.06 V) were noticed from 200 h to 1000 h. This difference may be partly attributed to a noticeable drop of open circuit voltage (ΔOCV = 0.02 V) that occurs in real SOFC at t = 200 h, as seen in Figure 8a in Ref. [10], and partly to some other degradation processes that possibly were going on concurrently. The ΔOCV was attributed to different agents. Since OCV = V_0,c − V_0,a depends on mass fraction of hydrogen (y_{0, H2}), referring to Equation (17), it can be speculated that y_{0, H2} was slightly decreased due to increased humidity in the fuel or hydrogen gas leak. To achieve this, the input molar fraction of hydrogen (x_{0, H2}) was adjusted to 0.955 in order to simulate ΔOCV at t = 200 h.

The voltage (V_da+ΔOCV) of SOFC, with degraded anode and ΔOCV, followed the reference voltage (V_{Ref. [10]}) from 0 to 200 h, but the difference between V_da+ΔOCV and V_{Ref. [10]} that remained was probably attributed to some other degradation processes. The grain coarsening possibly occurred in the SOFC cathode due to high temperature (900 °C) [24], although there was no data about it in Ref. [10]. In order to reproduce the higher degradation rate of V_{Ref. [10]} the grain coarsening was also modeled in the SOFC cathode (growth coefficient k_c = 1.4 × 10⁻²¹ m³ h⁻¹ and cathode exponential coefficient κ_c = 3) The values were manually fitted to achieve a minimum standard deviation between the voltage (V_da/c+ΔOCV) of SOFC, with degraded anode/cathode and ΔOCV, and V_{Ref. [10]}. It should be noted that using the model for online estimation of SOFC degradation would require similar fitting procedure. The parameters k_a/c and κ_a/c should be fitted automatically on the measured voltage profile within the first few 100 h intervals.

In Figure 7b, V_da/c+ΔOCV closely followed V_{Ref. [10]} if the degradation due to the grain coarsening was modeled in both electrodes. The maximum absolute difference between the V_da/c+ΔOCV and V_{Ref. [10]} was about 7 mV, whereas standard deviation was about 3 mV. Both values were negligibly small. The absolute degradation rate of V_{Ref. [10]} from 0 to 1000 h was estimated to be 0.14 V kh⁻¹, whereas the relative degradation rate was about 25% kh⁻¹. The V_da/c+ΔOCV had similar degradation rates. These rates were very high, most probably due to high operating temperature 900 °C, which was close to the melting point of the materials, so the agglomeration of grains progressed very quickly, as already argued in Ref. [10].

After the fitting procedure, the model was used to calculate the degradation of SOFC electric performance in the range of 10,000 h, which was a common operating period [24]. Two additional plots, as seen in Figure 7c,d, were added to show the calculated degradation of the modeled SOFC. It can be noticed that the relative TPB density of anode (cathode) dropped to about 0.18 (0.21), whereas the V_da/c+ΔOCV dropped to about 0.28 V at 10,000 h. In other words, in this period of time, the V_da/c+ΔOCV decreased by approximately 50% if considering the initial V_da/c+ΔOCV value of 0.56 V. To show how the degradation affects the electric conversion efficiency (η) of the SOFC, the η was evaluated at every Δt = 100 h from 0 to 1000 h, according to the following equation:

$$\mathsf{\eta }=\frac{P_{\text{el}}}{P_{\text{in}}}=\frac{P_{\text{el}}}{{\mathsf{\Phi }}_{\text{in},}{}_{{\mathrm{H}}_2}\cdot \Delta H_{{\mathrm{H}}_2}}=\frac{2\cdot F\cdot V_{\text{cell}}}{\Delta H_{{\mathrm{H}}_2}}$$

(24)

The P_el is electric power density, P_in is input power density, Φ_in,H2 is input molar flux of hydrogen. It can be noticed that η depends on the output voltage of SOFC (V_cell), Faraday constant (F), and standard enthalpy of formation for hydrogen (ΔH_H2) only if it is assumed that each oxidized H₂ contributes two electrons, Equation (2). Figure 8a shows the η of the SOFC as a function of time. The reference conversion efficiency (η_{Ref. [10]}) dropped rapidly from 0 to 200 h and moderately from 200 h to 1000 h (by approximately 0.07 and 0.04, respectively). A similar trend of the conversion efficiency (η_da/c+ΔOCV) of SOFC, with degraded anode/cathode and ΔOCV, was observed. The absolute degradation rate of η_{Ref. [10]} from 0 to 1000 h was estimated to be 0.11 kh⁻¹, whereas the relative degradation rate was about 25% kh⁻¹. Similar degradation rates of η_da/c+ΔOCV were noticed. These values were high, but it should be noted that the real SOFC degradation test was performed at high temperature (900 °C). The details about the experimental test can be found in Ref. [10].

It is important to note that the proposed model reproduces reliable results which are closely related to the measured ones. Figure 8b shows the electric power density (P_cell) and energy density (E_cell) of SOFC that can be calculated by Equations (25) and (26).

$$P_{\text{cell}}\left(t\right)=J_{\text{cell}}\left(t\right)\cdot V_{\text{cell}}\left(t\right)$$

(25)

$$E_{\text{cell}}(t)={\displaystyle \underset{0}{\overset{t}{\int }}{J}_{\text{cell}}\left(t\right)\cdot V_{\text{cell}}\left(t\right)}\cdot dt\prime \cong {\displaystyle \sum _{i=1}^{i=\mathrm{N}}{J}_{\text{cell}}\left[\mathrm{i}\right]\cdot V_{\text{cell}}\left[\mathrm{i}\right]\cdot \Delta t}$$

(26)

The energy yield of non-degraded cell (E_ref) equals 25.1 MWh m⁻², whereas the energy yield of degraded cell (E_da/c+ΔOCV) equals 15.8 MWh m⁻² at t = 10,000 h. The yield of electrical energy, generated by the real SOFC over a long time, would be much lower (by approximately 37%) than calculated by considering the initial electric performance (i.e., the non-degraded SOFC).

As a result, the overall cost per generated MWh of electric energy could be estimated when a high power SOFC generator was considered. This estimation is beyond the scope of the present study since the input data of the model is based on laboratory scale SOFC.

The model is very useful tool to predict the degradation rate and is suitable for online monitoring of the SOFC. However, the proposed Ni coarsening model should be improved in the future in order to consider the temperature, current density, and gas composition dependent growth of particles.

4. Conclusions

In this paper, the reduced-complexity, two-dimensional, microstructural model of solid oxide fuel cells (SOFC) was proposed since a numerically tractable model is sought for online evaluation of SOFC electric performance degradation. The SOFC degradation that occurs at high temperature due to the nickel agglomeration in the anode was modeled by using power law coarsening theory. The model was validated by comparing the electro-impedance spectra (EIS) of the modeled SOFC and measured EIS characteristics. The anode and cathode polarization resistances increased over the simulated time, which indicated the trends which are also observed in the real SOFC device.

The model was also validated based on the measurement data. The modeled structure was based on the backscattered electron (BSE) microscope cross-section picture of the measured SOFC. The input parameters were adjusted according to the known material parameters or fitted to match the relative density of TPBs profile and initial performances of the measured SOFC. The output voltage of the modeled SOFC closely follows the measured voltage from 0 to 1000 h when considering the degradation due to the grain growth in both electrodes. The absolute and relative degradation rate of electric conversion efficiency (η) equals 0.11 kh⁻¹ and 25% kh⁻¹, respectively. The electric energy yield (E_cell) of the degraded SOFC over the simulated time of 10,000 h was about 37% lower than it would be in the E_cell of non-degraded SOFC. Moreover, the overall cost per generated MWh of electric energy could be estimated when considering a commercial SOFC stack for high electric power.

Regarding the presented results, the model was found suitable for online estimation of the degraded SOFC electric performances. Future work is going to improve the proposed degradation model by adding the temperature, current density, and gas composition dependent growth of particles.