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

: 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-o ﬀ 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 ﬁt the measured data. The advantage of the proposed model is that it takes lower computational load during the simulation compared to the complex phase ﬁeld models and is suitable for estimation of SOFC electric performance over time.


Introduction
Solid oxide fuel cell (SOFC) systems are a promising technology for stationary applications since they convert hydrogen (H 2 ) 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 2 ). 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 4 -CO 2 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 −2 ), and the fuel utilization of 60%. The CH 4 and CO 2 flow rates were both 43 cm 3 min −1 , whereas the airflow rate was 750 cm 3 min −1 . The SOFC output voltage decreased different output current. The mass fraction of water (y H 2 O ) is calculated as 1 − y H 2 , whereas the mass fraction of nitrogen (y N 2 ) is calculated as 1 − y O 2 .

Conservation of Mass
The mass conservation equation is written in the following form: 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 −1 K −1 .
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.

Conservation of Mass
The mass conservation equation is written in the following form: The density of a gas is denoted with ρ, v is the average velocity of mass flux, Si 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 (yi). Ideal gas constant is denoted with R = 8.314 J mol −1 K −1 .
( ) 2 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.

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)(2)(3)(4)(5) in Figure 2 depicts the mechanism that is described in the equally numbered subsection.

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)(2)(3)(4)(5) in Figure 2 depicts the mechanism that is described in the equally numbered subsection.

Surface Exchange
It is assumed the oxygen (O2) that fills the pores of solid (i.e., mixed ionic-electronic conductive-MIEC)) material is adsorbed into the cathode surface. Oxygen ions (O 2− ) are formed with a double negative charge by accepting electrons (e − ) from the cathode. The flow of O 2− (JO2−) can be modeled as follows [16]: The surface exchange coefficient is denoted with ksurf. It depends on the temperature and the partial pressure of O2. The equilibrium oxygen ion concentration is denoted with CO2−, eq and depends on the partial pressure of O2 within pores. The concentration of oxygen ions (CO2−) decreases during the increasing electric current load (CO2− < CO2−, eq).

Bulk Diffusion
Since the electronic conductivity of the cathode MIEC material is high and should not limit the electric current, the diffusion of O 2− is assumed to limit the current. The diffusion of O 2− is modeled by using Fick's first law: Dchem is chemical diffusion coefficient that depends on temperature and partial pressure of O2 [17].

Charge Transfer
At the interface between the solid MIEC and electrolyte material the O 2− 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: The anode/cathode exchange current density is denoted with i0,a/c, forward/backward exponential coefficient of anode/cathode with αf/b a/c (for simplicity, αf a = 0.5, αf c = 0.5, αb a = −0.5, αb c =

Surface Exchange
It is assumed the oxygen (O 2 ) that fills the pores of solid (i.e., mixed ionic-electronic conductive-MIEC)) material is adsorbed into the cathode surface. Oxygen ions (O 2− ) are formed with a double negative charge by accepting electrons (e − ) from the cathode. The flow of O 2− (J O 2− ) can be modeled as follows [16]: The surface exchange coefficient is denoted with k surf . It depends on the temperature and the partial pressure of O 2 . The equilibrium oxygen ion concentration is denoted with C O 2− , eq and depends on the partial pressure of O 2 within pores. The concentration of oxygen ions (C O 2− ) decreases during the increasing electric current load (C O 2− < C O 2− , eq ).

Bulk Diffusion
Since the electronic conductivity of the cathode MIEC material is high and should not limit the electric current, the diffusion of O 2− is assumed to limit the current. The diffusion of O 2− is modeled by using Fick's first law: D chem is chemical diffusion coefficient that depends on temperature and partial pressure of O 2 [17].

Charge Transfer
At the interface between the solid MIEC and electrolyte material the O 2− 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: ; at TPBs in the anode.
; at TPBs in the cathode.
The anode/cathode exchange current density is denoted with i 0,a/c , forward/backward exponential coefficient of anode/cathode with α a/c f/b (for simplicity, α a f = 0.5, α c f = 0.5, α a b = −0.5, α c b = −0.5 are Appl. Sci. 2020, 10, 1906 6 of 18 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.
Partial pressures of oxygen, hydrogen and water are denoted with p O 2 , p H 2 , and p H 2 O , respectively, and the ∆G H 2 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.

Electronic Current
Electronic current is modeled in solid MIEC material of the cathode by the continuity equation [17]: The electronic conductivity of cathode MIEC material is denoted with σ s,c , δ ads denotes delta function: δ ads ( → r ) = 1; at MIEC/gas interface in cathode 0; elsewhere (19) Similarly, the electronic current is modeled in solid MIEC material of anode: The electronic conductivity of anode MIEC material is denoted with σ s,a , δ a denotes delta function.

Ionic Current
The ionic current in anode/cathode/electrolyte is modeled with the continuity equation [17]: The ionic conductivity of anode/cathode/electrolyte is denoted with σ ion,a/c/el , δ a/c denotes delta function.

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.

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

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 Appl. Sci. 2020, 10,1906 8 of 18 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.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 18 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.

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 O2 and H2 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 timeconsuming [23]. This is why the detection in our case was restricted to 2-D grey-scale images only.

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 2 and H 2 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.

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 nondegraded 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]: Time-dependent diameter of Ni grains is denoted with d, initial diameter with d0, 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 d0 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.

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]: Time-dependent diameter of Ni grains is denoted with d, initial diameter with d 0 , 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 0 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: 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].

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 shows the list of input parameters.

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.

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.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 18 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.

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 Since the open circuit voltage (OCV) of SOFC is about 0.97 V, the measured Vcell 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 Vcell, 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]. 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 2 at non-degraded and 1.4 Ω cm 2 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 −2 , 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).

DC Analysis
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].

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 −2 , whereas 3 vol% humidified H 2 was supplied to the anode and air was supplied on the cathode side. The flow rates of hydrogen and air were 350 cm 3 min −1 and 1500 cm 3 min −1 , 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 −2 , 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 −21 m 3 h −1 ) 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 −2 vs. 2.9 µm −2 ), but another study revealed [30] that the density of active TPBs may be lower (0.74 µm −2 ). 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, H 2 ), referring to Equation (17), it can be speculated that y 0, H 2 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, H 2 ) was adjusted to 0.955 in order to simulate ∆OCV at t = 200 h.
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 = V0,c − V0,a depends on mass fraction of hydrogen (y0, H2), referring to Equation (17), it can be speculated that y0, H2 was slightly decreased due to increased humidity in the fuel or hydrogen gas leak. To achieve this, the input molar fraction of hydrogen (x0, 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 −21 m 3 h −1 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 −1 , whereas the relative degradation rate was about 25% kh −1 . 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: The P el is electric power density, P in is input power density, Φ in,H 2 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 H 2 ) only if it is assumed that each oxidized H 2 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 −1 , whereas the relative degradation rate was about 25% kh −1 . 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].  The voltage (Vda+ΔOCV) of SOFC, with degraded anode and ΔOCV, followed the reference voltage (VRef. [10]) from 0 to 200 h, but the difference between Vda+ΔOCV and VRef. [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 VRef. [10] the grain coarsening was also modeled in the SOFC cathode (growth coefficient kc = 1.4 × 10 −21 m 3 h −1 and cathode exponential coefficient κc = 3) The values were manually fitted to achieve a minimum standard deviation between the voltage (Vda/c+ΔOCV) of SOFC, with degraded anode/cathode and ΔOCV, and VRef. [10]. It should be noted that using the model for online estimation of SOFC degradation would require similar fitting procedure. The parameters ka/c and κa/c should be fitted automatically on the measured voltage profile within the first few 100 h intervals.
In Figure 7b, Vda/c+ΔOCV closely followed VRef. [10] if the degradation due to the grain coarsening was modeled in both electrodes. The maximum absolute difference between the Vda/c+ΔOCV and VRef. [10] was about 7 mV, whereas standard deviation was about 3 mV. Both values were negligibly small. The absolute degradation rate of VRef. [10] from 0 to 1000 h was estimated to be 0.14 V kh −1 , whereas the relative degradation rate was about 25% kh −1 . The Vda/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 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).
The energy yield of non-degraded cell (E ref ) equals 25.1 MWh m −2 , whereas the energy yield of degraded cell (E da/c+∆OCV ) equals 15.8 MWh m −2 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.

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 −1 and 25% kh −1 , 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. Acknowledgments: The Slovenian Research Agency (ARRS) through Research Program P2-0001 and the project PR-07383 is gratefully acknowledged. The Author would also like to kindly thank Đani Juričić, for fruitful discussions and suggestions that improved the content of this paper.

Conflicts of Interest:
The Author declare that there is no conflict of interest regarding the content of this article. Mass fraction of the i-th gas specie η [/]

Abbreviation
Conversion efficiency