Multiscale Modeling of Degradation of Full Solid Oxide Fuel Cell Stacks

(cid:1) A multiscale degradation model for a full solid oxide fuel cell stack is developed. (cid:1) The model simulates 38 thousand hours of the stack lifetime in 1 h and 15 min. (cid:1) The model is validated with experimental data and predicts common long-term trends. (cid:1) Nickel coarsening, chromium poisoning, and interconnect oxidation are considered. (cid:1) Degradations under potentiostatic and galvanostatic operation modes are compared.


Introduction
Solid oxide fuel cell (SOFC) is a promising candidate for stationary power plants and mobile applications in constant operation (e.g. larger ships) due to its high efficiency and fuel flexibility. However, its market uptake is currently hindered by high costs stemming from its current low volume production, mechanical failures due to e.g. crack formation and creep, and degradation limiting its lifetime. The degradation could be induced mechanically [1] or chemically such as interconnect oxidation [2], carbon deposition [3], nickel particle coarsening [4], and the poisoning of the cathode and anode by chromium [5] and sulfur [6], respectively.
While degradation has been widely studied for single cells, where the conditions of the cell are well monitored, then the degradation inside a SOFC stack is more uncertain due to the various gradients of temperature, flows, and species concentrations [7]. Comprehensive analysis is essential to understand and mitigate degradation to improve the long-term performance of the SOFC. However, the interplay of the degradation phenomena, i.e. their multiphysics nature, is a hurdle against developing a precise prediction for each degradation phenomenon [8].
Many studies have been done on the degradation phenomena and their effects on SOFC performance. For instance, the following studies have been devoted to SOFC degradation due to the chromium (Cr) poisoning of the cathode [9]: discusses Cr poisoning mechanism and kinetics [10], investigates the effects of the electrodes and electrolytes compositions on the Cr degradation rate [11], studies Cr poisoning in porous perovskite cathodes [12], performs thermodynamic analysis of Cr poisoning [13], investigates the effects of the current density and cathode thickness on the Cr poisoning [14], studies Cr poisoning in different cathode electrode types [15], shows the effect of the humidity of the air supplied to the cathode side on the Cr poisoning, and [16] investigates the effect of the Cr concentration on the Cr poisoning rate.
The interconnect Cr, which is used to protect it against breakaway corrosion, is volatilized by the water vapor in the airflow and moves to the cathode electrode. Experiments conducted by Konysheva et al. [11,13] and Horita et al. [16] showed that Cr is distributed close to the triple-phase boundary (TPB) of the LSM-YSZ cathode electrode under a load current, while it is distributed uniformly under opencircuit voltage (OCV) operation. Moreover, the authors found that Cr poisoning of the cathode electrode should be through an electrochemical reaction of Cr at the TPB as they observed higher Cr deposition for higher activation overpotentials. Nakajo et al. [17] developed a model for the Cr poisoning that captures the cathode electrode TPB reduction based on the activation overpotential and reaction rate of the Cr oxide deposition. Miyoshi et al. [18] applied the Cr poisoning model to a 3D model consists of the electrolyte and anode and cathode electrodes. It was demonstrated that the model predicts common trends for Cr poisoning reported in the experimental studies. Moreover, it was shown that Cr poisoning increases the ohmic overpotential in addition to the activation overpotential, which is due to the shifting of active reaction zones from the electrode/electrolyte interface to the outer surface of the electrode. Nickel (Ni) particle coarsening is the primary degradation phenomenon for the anode side of the SOFC [19,20]. The Ni coarsening is due to the tendency of the particles to reduce free surface energy, Ostwald ripening mechanism, which reduces the TPB length and electronic conductivity of the anode electrode. Ni coarsening and its effects on SOFC performance have been demonstrated in experimental studies, e.g. Refs.
[21e23]. Vassen et al. [24] developed the first Ni coarsening model that predicts Ni particle size increase based on their surface diffusion. Nakajo et al. [17] used a simplified form of the model developed in Ref. [24] with the effect of the water vapor mole fraction. It was shown that the model predicts common trends for the Ni coarsening, except for very low fuel to water vapor ratios. Other Ni coarsening models could be found in the literature, e.g. Refs. [25,26]. However, they are claimed to be insufficient as the fitted Ni surface diffusion coefficient is two to three orders of magnitude smaller than the experimental ones [27]. Gao et al. [27] developed an analytical model with two fitting parameters, which could be reduced to one, for the Ni coarsening. The model was validated against experimental data sets with very high accuracy.
Another important degradation phenomenon is oxidation (corrosion) of the interconnect (IC). Metallic ICs are used for the SOFC due to their low cost and high thermal conductivity [28]. The typical metallic IC steels, e.g. Crofer22APU, contain chromium [29,30], which forms a chromium layer at the IC surface that retards its further corrosion. However, the chromium oxide scale increases the ASR, as it has a lower electronic conductivity. Oxide scale growth over the IC is typically modeled with Wagner's law [31,32].
While the different phenomena have been studied at electrode and electrolyte scale [18], cell scale [19], or at best in a repeating unit of a stack [17], then transient simulations of degradation at different locations in a full SOFC stack have not been computationally feasible. A single steady-state simulation with conventional stack models takes in the order of 40 h [7,33]. Thus, solving a transient model with such approaches would take in the order of 600 h, depending on the operational conditions and number of iterations for each time step, etc. This would thus hinder any practical usage of such a stack model. However, recent developments of a multiscale model, using so-called homogenization in the solid oxide cell (SOC) stack modeling, have significantly reduced computational resources and a consequent increase of speed [34].
This paper concerns the development of a multiscale degradation model for a full SOFC stack, which can simulate the whole lifetime of a stack, e.g. 38 thousand hours within 1 h and 15 min on a workstation. This is done by taking advantage of a recently developed multiscale modeling approach for the multiphysics stack model. Indeed, this study aims to show a route for investigating the long-term performance of the stack that can be simulated on a desktop workstation in a relatively short time while including the known degradation phenomena reported in the literature. The focus of this work is developing the computational method, and we have thus used the developed and verified degradation models from the literature [17e19] with their parameters for the electrodes and IC characteristics. The method allows for exchanging any i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( x x x x ) x x x other model into the framework. Using the individual degradation models to predict a convoluted response of the stack is reasonable as similar degradation trends are expected at cell and stack levels [35,36]. However, other stack-scale induced degradation phenomena might be involved over the longterm performance of the stacks, which could be added to the model whenever they are understood in the future.
It will be shown that the model predicts common longterm trends reported in the literature and captures the evolutions and interactions of different degradation phenomena over the whole stack for the first time. Moreover, it is demonstrated how such a model can be used to investigate the effects of different operation modes and conditions on the long-term performance of the stack, which could be a useful guide for the SOFC stack developers.
The paper is organized as follows: In Section Methodology, the approach for integrating the models of the degradation phenomena into the multiscale model of the SOFC stack is explained. The multiscale model includes the transport equations of mass, momentum, species, charge, and heat. In addition to the previous work [34] on the homogenized stack model, this model also considers laminar and turbulent flows in the fuel and air manifolds, respectively, and their coupling to Darcy's Law in the header and active domains. In Section Validation, model validation is performed and discussed. In Section Results, the evolutions of the stack performance and degradation phenomena are investigated for different operating parameters of galvanostatic and potentiostatic operation modes, operation conditions (inflow temperature and load current), and flow configuration (co-and counter-flow).

Methodology
Here, the multiscale model of the SOFC stack is described. Section Multiscale modeling concept discusses the concept of the multiscale model used in this study, including an overview of the previous model development. Detail of the studied stack is given in Section Details on the studied stack. Section Making the homogenized stack model transient presents the development of the transient model for the stack. The free flow in the manifolds and the coupling to the flow in the headers and active domain are described in Section Flow in the manifolds and coupling to the flow in the headers. The cell model is described in Section Cell model. Finally, the integration of the models for the degradation phenomena of Ni coarsening, Cr poisoning, and IC oxidation into the multiscale model is described in Section Integration of the degradation models.

Multiscale modeling concept
The stack model is based on a multiscale modeling approach for the transport equations of mass, momentum, species, charges, and heat over the active domain, headers, and sealing domains of the stack, shown in Fig. 1c. To understand the multiscale modeling approach, this can best be explained by considering how electrodes in cell models are modeled today. Rather than modeling all the phases of the electrode materials including the pores, an effective medium approach is taken. Hence, effective diffusion coefficients and conductivities for the porous microstructure are used, as shown in Fig. 1aʹ and [33,34,37,38], rather than actual porous geometry, shown in Fig. 1a. In some cases, the actual geometry, Fig. 1a, is modeled to obtain these effective parameters, e.g. Refs. [39,40]. However, the full cell is not modeled with all the geometrical details, as this would be an unnecessary computational complication. The multiscale stack model used in the current work applies the same principle e rather than modeling the details of the repeating unit (Fig. 1b), effective parameters are prescribed for an effective (homogeneous) medium (Fig. 1bʹ). Eventually, these effective parameters will be used to represent the interior of the stack (Fig. 1c), as explained in detail in Ref. [34].
An effective medium for each of the stack domains must thus be provided (active domain, headers, and sealing domains). A detailed discussion of the governing equations for the active domain is given in Ref. [34]. The same homogenization approach is used for the extension of the model to the headers and sealing domains. The variable definitions, e.g. permeabilities and thermal conductivities, are updated with the dimensions of the headers and sealing domains, which are given in Fig. 2. More details can be found in our recent work [37]. The same partial differential equations (PDE) are thus used for the mass and heat transport models in the manifolds as in the active domain. The different geometries do, however, result in different effective properties (permeability, etc.).
Degradation depends on local parameters inside the cell, typically at interfaces or junctions or in other components such as the interconnect. First-principles models should be used to deeply understand these degradation mechanisms, but this is far too circumstantial for practical use when designing stacks. Therefore, a degradation model is typically obtained from an experiment with descriptive parameters correlating to the observed behavior, as in the works described in the introduction. For instance, the activation overpotential, h act , in the air electrode heavily influences the Cr poisoning [17]. In the experiments where such correlations have been extracted, the overpotential is the effective overpotential determined by electrochemical impedance spectroscopy (EIS), which describes the effective reaction over several sites through the thickness of the air electrode (Fig. 1b).
In the framework here laid out, this overpotential is determined using a cell model embedded in the homogeneous medium for the repeating volume, Fig. 1bʹ. Thus, the overpotential is calculated based on the local conditions at the particular point in the stack (temperature, gas compositions, and current). But, as these conditions are available at every point in the stack, the local overpotentials will also be available at every point. If we zoom sufficiently into the repeating unit, we know that this is not the case. On the other hand, we also know that the variation of e.g. temperature does not vary significantly inside the stack over the mmsized repeating unit. Thus, this approach will provide us with a first overall estimate of the parameters needed for a stack degradation model, as a degradation simulation with a fully geometrically resolved stack model is not computationally feasible today.  . 2 shows the schematic of the repeating unit and the full stack, including the manifolds, active domain, headers, and sealing domains. The right half of the stack in the x-direction is modeled due to the symmetry of the stack configuration in the xdirection. The air enters and exits the stack through the front and back manifolds, respectively (blue lines in Figs. 1 and 2). Fuel enters and exits the stack through the side manifolds at the front and back, respectively (red lines in Figs. 1 and 2). Thus, we have a co-flow configuration for the base case. It should be mentioned that inlets and outlets are located at the bottom of the manifolds.

Details on the studied stack
The cell in the stack is from Forschungszentrum Jü lich and is represented through a model developed by Leonide et al. [41,42]. The base operating conditions and geometry parameters are listed in Table 1. We considered a stack of 100 cells, which is in the range of the modern stack designs. While some manufacturers tend to work with shorter stacks to minimize scrap in the production due to mechanical failures, then other works with various stack sizes. For instance, the company FuelCell Energy has developed stacks with 45, 150, and 350 cells [43]. Concerning model development, the number of cells only affects the flow regime in the manifolds as a higher number of cells results in higher flows of reactants that need to be supplied to the cells. The flow in the air manifolds is affected more, as the air flow is also used for cooling and blown in excess to the stoichiometrically needed, as opposed to that of the fuel side, where fuel utilization is pushed as high as possible. The 100-cell stack model thus has a turbulent flow in the air manifolds. However, for the shorter   Making the homogenized stack model transient As the applied framework for establishing the model, COMSOL Multiphysics, has transient versions of all the PDEs (called "physics" in COMSOL). Hence, once the multiscale approach has been established and effective parameters determined, then setting the transient version of the model and those representing all homogenized representations of the phenomena in the stack is straightforward. However, several considerations should still be taken for the specific application of simulating the degradation of the stack transiently. These are described here.
The transient governing equations for the active domain, headers, and sealing (but not the manifold) domains are list in Table 2. In these domains, the flow is friction-dominated, as explained in Ref. [34], and can thus be represented by Darcy's Law, Eq. (1). This has a single form relating the average velocity to the pressure drop over the domain. Moreover, it has to be mentioned that the transient model is developed based on the transient terms, which can be found on the left-hand side of the governing equations, as shown in Table 2.
No special treatment is applied to the governing equations to capture fast transients such as fast electrical impedance, which must be included for the sub-second transients (high-frequency changes). For instance, a constant relative permittivity of one is used, which removes the fast transients through the remanent electric displacement, term D in Eq. (4). This is reasonable as we are only considering slow changes occurring over thousands of hours. If fast electrical responses were to be modeled, an equivalent circuit with more terms needs to be included.
For the thermal response, the transient model here also includes the thermal mass term and thus the heat capacity of the stack. As seen later, we are operating the stack with constant load current or constant stack potential, not studying e.g. the effect of shutdowns. Therefore, with the slow changes over time from degradation, it would be a good approximation not to exclude the thermal mass as well. With the inclusion, the model can capture changes over minutes and hours, but not sub-seconds. Since this study aims to model the degradation phenomena over the lifetime of the stack (i.e. thousands of hours without fast transients), the governing equations in their transient forms indicated in Table 2 are thus sufficient.

Flow in the manifolds and coupling to the flow in the headers
The homogenization approach is not used for the manifolds as they are free domains without sectioning, and the flow is not friction-dominated as in the interior of the stack, so the full Navier-Stokes equations must be solved. Thus, the Species Heat i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( x x x x ) x x x governing equations for the free flow in the manifolds (laminar or turbulent flow) differ from the header/active domain, where Darcy's Law is used to effectively describe the laminar flow in the channels inside the stack [34]. Therefore, different flow models are used for the manifolds and their adjacent header domains, and continuity of the pressures and mass conservation on their common interfaces are done as described in the following.
To ensure continuity of the flows over the shared interfaces of the free and porous mediums, boundary conditions are defined on these boundaries for each flow regime based on information (velocity or pressure) from the adjacent flow regime. For the porous mediums, the pressures on their common boundaries with the free flows are set equal to the pressures of the adjacent free flows on these boundaries. Moreover, velocities of the free flows on the common boundaries with the porous mediums are calculated based on the velocities of the porous mediums on these boundaries through mass conservation of the reactants crossing these boundaries.
In the header and active domains, the different gas species only occupy a part of the volume of the repeating element. To represent this in the homogenized model, where the gas species occupies the whole volume, the density of the different species is adjusted accordingly. Thus, at the boundary between manifolds and headers, the mass flux is made equal by adjusting the velocity. This corresponds to the actual change in velocity going from a free domain (manifold) to more confined flow channels in the header. Therefore, at the interface between the manifolds and headers, the velocity of the free flow is coupled through the following relationship: Navier-Stokes (NS) equations are used to model flow distribution in the fuel manifolds as the flow regime within them is laminar with a Reynolds number (Re) of about 57. NS equations are: where u is the velocity vector, and r and m are the density and dynamic viscosity of the fluid, respectively. The airflow in the air manifolds is turbulent as its Re is about 5400. The k-u turbulence model with Wilcox revised formulations [44] is used to model the turbulent flow in the air manifolds. The k-u model is a Reynolds-averaged Navier-Stokes (RANS) model that solves the NS equations coupled to two extra partial differential equations for the turbulence kinetic energy, k, and specific rate of dissipation of the kinetic energy, u, to capture the flow turbulence. The k-u model equations are as follow: Here, m T is the eddy viscosity, which is defined as m T ¼ r k u ; and a, b 0 , b * 0 , s u , and s * k are constants, which are set to 0.52, 0.072, 0.09, 0.5, 0.5, respectively. It should be noted that these are the default values considered in COMSOL Multiphysics. Moreover, default COMSOL wall function formulations are used, which ignore the flow in the buffer sublayer and consider a non-zero velocity at the walls based on an analytic solution for the flow in the viscose sublayer. Table 3 summarizes the governing equations for the cell model. The cell model base is the same as the one presented in our previous study [34]. However, in this study, definitions of the activation and ohmic overpotentials are changed to capture the effects of the degradation phenomena on cell performance. Activation overpotentials are calculated from the Butler-Volmer (BV) equations (15) and (16), which capture the effects of the TPBs reductions given in equations (28) and (31) due to Ni coarsening and Cr poisoning, respectively. Moreover, the ohmic overpotential is defined as a function of the conductivities of the electrodes and electrolyte materials, Eq. (21), to capture the effect of the Ni coarsening on the reduction of the electrical conductivity of the anode electrode, Eq. (24).

Cell model
It should be mentioned that we do not consider the extra resistances against the transport of the species to the reaction sites under the rib regions and the transport of the charges from these sites. In other words, we suppose that all the reactions occur over an electrode active area placed under the channel. This is a reasonable assumption for the trapezoidshaped channel geometry with the high channel to rib ratio used in this study, as shown in Fig. 2, in addition to the modern wavy channels used in SOC stacks. Nonetheless, the overpotentials will be improved in our future studies to make the model more accurate, especially for the stacks with rectangular channels.

Integration of the degradation models
Degradation phenomena considered in this development of the model are Ni particle coarsening in the anode electrode, Cr poisoning of the cathode electrode, and oxidation of the IC. The analysis is done for an intermediate-temperature SOFC stack with Ni-YSZ, YSZ, and LSM-YSZ materials for the anode electrode, electrolyte, and cathode electrode, respectively. The LSM-YSZ electrode is not used in modern cell designs, but a very well-developed model of its degradation is available in the literature [17,18], as opposed to that of the LSCF.
Governing equations for the degradation phenomena of Ni coarsening, Cr poisoning, and IC oxidation are summarized in Table 4. The values of the parameters used for the BV for the anode reaction rate [45] i an ¼ i 0;an l an;eff exp 2F BV for the cathode reaction rate [46] i Anode exchange current density [47] i 0;an ¼ 31:4 p À0:03 Cathode exchange current density [48] i 0;ca ¼ 2:14 Â 10 5 p 0:376 Anode concentration overpotential [41] Cathode concentration overpotential [41] Ohmic overpotential [49] h Electronic conductivity of the anode s e;an ¼ s 0 Average particle coordination number TPB density of the anode electrode l an ¼ 2pminðr Ni ; r YSZ Þðsin qÞ 3ð1 À 4Þj Ni 4pr 3 Effective TPB density of the anode electrode Cr poisoning [18] Cr poisoning reaction TPB density of the cathode electrode 1 l ca Reaction rate of the Cr oxidation IC oxidation [17] Wagner's law for the oxide scale growth vh 2 Conductivity of the oxide scale ASR of the oxide scale ASR IC;ox ¼ h IC;ox s IC;ox (35) i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( x x x x ) x x x degradation models are listed in Table 5. It should be noted that the degradation models are local models not averaged over the modeling domain, so the parameters given in Table  5 are for the electrodes and IC, i.e. these are not effective parameters. Therefore, the homogenization approach is just applied to the multiphysics model to resolve the effective modeling variables, such as temperature, for the degradation models.

Nickel particle coarsening
Ni coarsening is described by an analytical model, Eq. (22), based on surface diffusion between adjacent Ni particles [19,27]. This model predicts experimental data very well for various operating conditions via fitting parameters b and l. A percolation theory, described in Ref. [50], is used to determine effects of the Ni coarsening on the TPB density, Eq. (28), and the electronic conductivity of the anode electrode, Eq. (24).

Chromium poisoning
An electrochemical deposition due to oxidation of volatile Chromium species at the cathode electrode TPB is believed to be the main reason for the Cr poisoning degradation [18]. The reduction of the TPB of the cathode electrode due to the Cr oxide deposition is modeled via an empirical relation, Eq. (31), based on the reaction rate of the Cr oxidation obtained from the BV equation, Eq. (32).

Corrosion of the interconnect steel
Wagner's law, Eq. (33), is used to calculate the growth of the oxide scale on the IC. Area-specific-resistance (ASR) of the oxide scale, Eq. (35), is calculated by the ratio of its thickness and conductivity, Eq. (34). Parameters used for the IC degradation model are taken from Ref. [17] for Crofer22APU. Overall, the end-of-life (EOL) of the stack is dominated by the accelerated degradation rate of the Cr poisoning for the operation condition considered in this study, as shown and discussed in section Overall evolutions of the degradation phenomena and related model variables. However, other EOL metrics are necessary for the operating conditions with alleviated Cr poisoning, e.g. higher operation temperature presented in section Effects of the gas inflow temperatures. A proper metric could be the Cr content of the IC as it forms a chromium layer over the IC to protect it against corrosion. The protective chromium layer losses Cr through the volatilization of CrO 2 (OH) 2 , which is the source for the Cr poisoning of the cathode electrode. It has been observed that the drop of the Cr content of the IC below 16% (breakaway limit) leads to oxidation of the iron content, which is an accelerating process of IC oxidation [51,52]. Conservation of the Cr is used to calculate its remaining content in the IC: where n is the mole number, x the mole fraction, r the density, h the height of the respective layers, A IC the surface of the IC, M the molar mass, and subscripts rem, i, and o denote remaining, initial, and out, respectively. Two assumptions are considered for deriving the above equations: 1) height of the IC does not change due to loss of its Cr content, and 2) a uniform oxide scale with a height of h Cr 2 O 3 forms over the IC. These assumptions are not quite accurate but sufficient for a reasonable estimation of the remaining Cr content in the IC. Substituting Eqs. 37e39 in Eq. (36) gives the remaining Cr content in the IC: It should be noted that height of the oxide scale, h Cr 2 O 3 , is calculated by Wagner's law, Eq. (33).

Numerical approach
COMSOL Multiphysics software is used to couple the governing equations for the electrochemical and degradation models and solve them numerically via the FEM approach. The right half of the stack in the x-direction is modeled, as shown in Fig. 3, to take advantage of symmetry to reduce the computational cost of the model.
A non-uniform quadrilateral mesh is used for the top surfaces of the active domain, headers, and sealing domains, and a triangular mesh is used for the top surfaces of the manifolds and sealing domains adjacent to the fuel manifolds, Fig. 3. This mesh is swept through the modeling domain non-uniformly. A relatively coarse mesh density in the z-direction is used as it does not influence the important quantities for the degradation (mass, current, and temperature distribution in the active domain). It should be sufficient for the manifolds as the flow is through a straight channel with a fixed cross-section and under fully developed conditions. Moreover, a benefit of using RANS models, including komega, is that they work fine with common meshes. A denser mesh is considered close to the edges of the modeling Fig. 3 e Mesh distribution.   Mesh-independency of the results is ensured as running the model for a finer mesh with 730 thousands of DOF, i.e. more than double the number of DOF, results in very low relative changes of 0.04, 0.34, and 0.06% for the stack voltage, maximum current density, and maximum temperature, respectively.
The Fully Coupled solver (in Comsol), i.e. coupling all the governing equations of the electrochemical and degradation models and solving all unknowns simultaneously, leads to convergence issues. Therefore, the Segregated solver (in Comsol), which subdivides the problem into several steps and solves them sequentially, is used to solve the model. The steps Fig. 5 e Model validation against the experimental data from Nishida et al. [33] for the 18-cell FZJ Mark-F SOFC stack: (a) comparisons of the polarization, temperature, and power curves for the power curve operation given in Table 6; and (b) comparisons of the temperature distributions along the section A-Ashown in Fig. 4 and for the steady-state operation given in Table 6. The abbreviations exp and sim denote experiment and simulation data, respectively. Voltage subscripts represent the cellsʹ number, 1 for the bottom cell and 18 for the top cell. Locations of the temperatures T 9 and T 10 are shown in Fig. 4. (c) Fuel stream-tubes colored by hydrogen mole fraction for the steady-state operation given in Table 6. i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y x x x ( x x x x ) x x x considered for the Segregated solver are: 1) flow distributions for the fuel and air; 2) turbulent variables for the airflow in the air manifolds; 3) mass transport of the species; 4) current, voltage, overpotentials, and degradation variables; and 5) heat transfer in the solid and fluids. The system of the equations at each step is solved by Newton iterations with direct solver PARDISO (parallel sparse direct solver).
A stationary run of the model without degradation submodels is used as an initial condition for the transient run of the model. The stationary model could not be run directly for the operation condition considered in this study, i.e. 800 Pa (8 mbar) airflow and 70% fuel utilization. Hence, initialization steps are considered for the stationary run. First, the model is run for 10 Pa (0.1 mbar) airflow and a load current of 1 A. Afterward, they are increased gradually to their final values of 800 Pa airflow and load current of 25 A.
The final run of the stationary model is used as an initial condition for further changes in the operating conditions in the stationary model and transient runs of the model. Stationary runs of the model with small changes in the operation conditions take less than 5 min. Moreover, a transient run of the model for the base operation condition, listed in Table 1, takes 1 h and 15 min for 38 kh of transient operation of the SOFC stack. The time 38 kh is right before the rise of convergence issues due to the accelerated degradation of the stack, mainly the intense increase of the cathode activation overpotential, details in section Overall evolutions of the degradation phenomena and related model variables and Fig. 6.
It takes a long time, e.g. up to a day, to resolve the severe voltage drop/high resistance regime, i.e. the failure of the stack over a short time, e.g. 2 kh? Resolving this critical regime of the stack operation for a long time also has little practical relevance, as the divergence is a sign of fast degradation and thus also the approach towards the end of life for the stack. Even if the model does not crash, it is not essential to describe the stack behavior in this final failure evolution in detail. Instead, it provides an overall estimate of the lifetime, which a certain design and operating condition lead to. This study aims to show that our stack-scale model can predict the common degradation trend of initial decay, linear regime, and predict the initiation of the final accelerated regime [53]. Therefore, the transient runs of the model are stopped after a day if they have not crashed. The EOL is thus very precisely defined here, as this would be up to the end-user of the technology to specify their requirements. As some of the fast degradations can be simulated, it is however Fig. 6 e Evolutions of (a) stack voltage and volume averages of (b) activation overpotential of the cathode electrode, (c) activation overpotential of the anode electrode, (d) Ni particle radius, (e) normalized TPB density of the cathode electrode, (f) normalized TPB density of the anode electrode, (g) electronic conductivity of the anode electrode, (h) ASR of the oxide scale over the interconnect, and (i) comparisons of ASRs of the activation overpotentials of the cathode and anode electrodes and oxide scale over the interconnect.
certain that the EOL used here is not just because of some other numerical instability.
The model run times provided here refers to simulation on a high-end workstation with an Intel Core i9 3.5 GHz 12-core processor.

Validation
Validation of the stack model In the experiment, a simulated reformate gas with 10% pre-reforming of liquefied natural gas (LNG), corresponding to 10% humidified hydrogen, is considered for the anode side, and the air is supplied to the cathode side. In this section, the operating conditions and dimensions of the modeling domain components (active area, headers, sealing, and manifolds) are taken from Ref. [33]. Thus, a model of the stack presented in Ref. [33] has been generated with dimensions for the stack components corresponding to this stack, as shown in Fig. 4.
It should be noted that dimensions of the cell components (channels, ribs, supporting layers, electrodes, and electrolyte) are not given in Ref. [33]. Therefore, the FZJ cell geometry and characteristics are taken from Refs. [41,42], which is also used in the degradation study in Section Results. In this model, the flow inlet boundary conditions are modeled as average inlet velocities based on the flow rates of the species at the inlets [33].
Moreover, since Re is less than 2300 in the air manifolds for both operating conditions considered in the experimental   7 e Evolutions of spatial distributions of (a, aʹ, and aʺ) current density, (b, bʹ, and bʺ) activation overpotential of the cathode electrode, (c, cʹ, and cʺ) TPB density of the cathode electrode, (d, dʹ, and dʺ) activation overpotential of the anode electrode, (e, eʹ, and eʺ) TPB density of the anode electrode, and (f, fʹ, and fʺ) stack temperature. Columns from left to right are for time instances of 0, 20, and 37.6 kh, respectively. data, the flow in the air manifolds is modeled as laminar flow. Table 6 lists the main operating parameters for the experiments in Ref. [33], which are used for the model validation. Fig. 5a shows the comparisons of the polarization, temperature, and power curves from the model and experimental data for the power curve operation given in Table 6. A very good agreement is seen between the experimental and simulation results for lower current densities with a slight deviation for higher current densities. The maximum relative errors for the temperatures and power data are about 4.5 and 5.8%, respectively. Fig. 5b shows the comparison of the temperature distribution along section A-Ashown in Fig. 4 and the steadystate operation given in Table 6. It is seen that the model predicts the temperature trend along the stack but underestimates it with the highest deviation of 2.9% close to the fuel inlet. These results are reasonable as the modeling was performed without attempting to calibrate the results to the specific cells used in the experiment and without detailed knowledge about the stack design. The differences would, however, be interesting to study further in a consequent study.

Validity of the degradation models
Each of the degradation models presented in Section Integration of the degradation models is validated as far as possible in the respective sources. The fidelity of the current stack model relies on these sources. As mentioned in the introduction, it is not the purpose of this work to propose alternatives to the existing component degradation models but rather demonstrate that simulation of a full stack in 3D over ten thousand hours is possible. It should be mentioned that the degradation models are local models not averaged over the modeling domain, i.e. the degradation models are added Fig. 9 e Evolutions of spatial distributions of (a, aʹ, and aʺ) current density, (b, bʹ, and bʺ) activation overpotential of the cathode electrode, and (c, cʹ, and cʺ) stack temperature under potentiostatic operation mode. Columns from left to right are for time instances of 0, 20, and 37.6 kh, respectively. to our multiphysics model as the form developed in their corresponding references with the characteristic parameters for the electrodes and IC, given in Table 5.
Stack experiments are expensive and complex to interpret. No literature on the response of the interior of a stack tested, as in Ref. [54], for degradation is currently available in the literature. This type of data is not even available to studies limited to cells or sub-cell layers, e.g. Refs. [18,19]. For an absolute validation, the tremendous task of producing well-characterized and consistent cell and stack degradation experiments must thus be conducted. However, in the meantime, models as this can be used to assess the effects of operating a SOFC stack for various conditions by combining existing models.
Here, it is assumed that the components experience similar degradation in a single component test as in the stack. This is consistent with the experimental data reported for the SOFC stack degradation, e.g. Haart et al. [35] and Comminges et al. [36], showing similar degradation regimes of initial decay, linear regime, and final accelerated performance drop. Special degradation mechanisms, e.g. pronounced corrosion at the contact between the cell and interconnect, may exist but are poorly understood and documented. These can be included once adequately understood. Indeed, as discussed further in Section Results, the stack model in this work captures all the common degradation mechanisms reported in the literature [17e19].

Results
Overall evolutions of the degradation phenomena and related model variables Fig. 6 shows the evolutions of the stack voltage and volume averages of the degradation phenomena and model variables affected by them. The stack voltage experiences degradation in three steps: initial decay, linear regime, and final accelerated regime, as shown in Fig. 6a. The initial decay is due to an increase of ASRs of the activation overpotential of the anode electrode, Fig. 6c, and oxide scale over the IC, Fig. 6h. It should be noted that the reduction of the electronic conductivity of the anode electrode (Fig. 6g) contributes to this initial decay. However, its contribution is negligible compared to others due to the high conductivity level of the anode electrode even after its degradation through the Ni coarsening. Fig. 10 e Comparisons of the evolutions of (a) stack voltage and volume averages of (b) activation overpotential of the cathode electrode, (c) activation overpotential of the anode electrode, (d) Ni particle radius, (e) normalized TPB density of the cathode electrode, (f) normalized TPB density of the anode electrode, (g) electronic conductivity of the anode electrode, (h) ASR of the oxide scale over the interconnect, and (i) minimum of the remaining chromium in the interconnect for the inflow temperatures of 700, 730, and 760 C.
Anode activation overpotential increases rapidly at the beginning and saturates at a steady level, which follows an opposite trend for the anode TPB (Fig. 6f) due to the Ni particle growth shown in Fig. 6d. ASR of the oxide scale over the IC shows an almost linear behavior after a rapid increase initially, as indicated in Fig. 6h and [19,32], which corresponds to the linear regime of the stack voltage drop.
The final accelerated drop of the stack voltage is due to an exponential behavior of the cathode activation overpotential, Fig. 6b. The reason for such a behavior lies in the fact that the cathode activation overpotential is a direct function of the current density and reciprocal of the cathode TPB density, Eq. (16). Cr poisoning in the cathode electrode reduces the cathode TPB density, Eq. (31) and Fig. 6e, while the average current density is fixed. Fig. 6i shows comparisons of the volume averages of ASRs for the activation overpotentials of the anode and cathode electrodes and oxide scale over the IC. It is seen that the dominant ASR is the one from the cathode activation overpotential. Therefore, Cr poisoning is the primary degradation phenomena for the operating conditions considered here, which is consistent with the experimental data reported in Ref. [55] for a low operating temperature of 750 C. Percentage of the ASRs corresponding to the degradation phenomena of Cr poisoning, IC oxidation, and Ni coarsening at the EOL of the stack are about 87, 12, and 1%, respectively.
All in all, it can be concluded that the model predicts evolution trends for voltage, degradation phenomena, and related model variables consistent with common behaviors reported in the numerical studies [18,19,53] and experimental studies [17,35]. Fig. 7 shows the evolutions of spatial distributions of the current density, stack temperature, activation overpotentials, and TPB density of the anode and cathode electrodes at time instances of 0, 20, and 37.6 kh. At the beginning, the current density is higher at the inlet, i.e. front face of the active domain where fuel is supplied to the stack, and reduces towards the outlet, as shown in Fig. 7a. Higher current density results in higher activation overpotential at the inlet. This can be seen clearly for the anode activation overpotential, Fig. 7d, but not for the cathode activation overpotential, Fig. 7b, since the scale is fixed for the highest level at the EOL of the stack, which is about 30 times higher than the initial level for the cathode activation overpotential.

Spatial distribution of the degradation phenomena and related model variables
Higher cathode activation overpotential leads to higher Cr poisoning close to the inlet, Eq. (32), which reduces the cathode TPB density at the inlet, Eq. (31) and Fig. 7cecʺ. The reduction of the cathode TPB density at the inlet due to Cr poisoning results in a reduction of the current density at the inlet, as shown in Fig. 7aʺ. The vertical variation of Cr poisoning is primarily due to the vertical temperature variations seen in Fig. 7fefʺ.
To deliver a constant load current, the current density needs to be increased locally in the middle of the stack to compensate for the reduced current density at the inlet. However, the cathode activation overpotential still increases at the inlet owing to the reduction of the cathode TPB density, which is about 9 orders of magnitude higher than the reduction of the current density.  Unlike the high increment of the cathode activation overpotential, which is due to a significant reduction of the cathode TPB density, the anode activation overpotential does not increase much. Equations (15) and (16) describe how the anode and cathode activation overpotentials are direct functions of the current density and reciprocal of the TPB density. For the anode side, the activation overpotential decreases again after some time at the inlet and follows the current density evolution, Fig. 7dʺ. This is due to the saturation of the anode TPB density at a level half of its initial level, as illustrated in Fig. 7eeeʺ, which makes the anode activation overpotential only a function of the current density. Fig. 7fefʺ show the evolution of the stack temperature. It is seen that degradation does not significantly change the temperature distribution pattern and increases the maximum temperature by about 5 C, which is due to the increase of the current density. The temperature variation does, however, influence different degradation rates.

Galvanostatic vs potentiostatic
To the best of the authors' knowledge, most of the SOFC degradation tests are done under galvanostatic operation mode, which extracts a constant load current from the stack. However, applying a constant voltage to the stack, i.e. potentiostatic operation mode, can enhance the stack's lifetime, as shown in Fig. 8. Running the stack under the potentiostatic operation mode alleviates the increase of the cathode activation overpotential. Fig. 9 shows the evolutions of spatial distributions of the current density, cathode activation overpotential, and stack temperature under potentiostatic operation mode at time instances of 0, 20, and 37.6 kh. Comparison of the Fig. 7bʺ and Fig  9bʺ indicate that the maximum level of the cathode activation overpotential is reduced about 50 mV for the potentiostatic operation mode, while the evolution trend for the cathode activation overpotential is the same for both potentiostatic and galvanostatic modes ( Fig. 7bʺ and Fig 9bʺ). The reason for this similar behavior is the evolution of the current density, Fig. 9aeaʺ.
The current density follows a trend similar to the one for the galvanostatic mode, i.e. higher drop at the inlet, but it reduces everywhere in the stack for the potentiostatic mode as the stack is not supposed to deliver a constant load current. This higher current density reduction results in a higher stack power reduction for the potentiostatic operation mode, as shown in Fig. 8i. Potentiostatic operation mode improves the stack lifetime, about 10%, while lessens the stack output power more than the galvanostatic operation mode. The stack power output is 1.89 and 2.15 kW at EOL under potentiostatic and galvanostatic operation modes, respectively. However, the total energy generated, integral of power over the stack lifetime, under potentiostatic and galvanostatic operation modes are 86.7 and 84.3 MWh, respectively. This indicates that the gain in the stack lifetime under the potentiostatic operation mode is high enough to compensate for the lower power generation of the stack and even improve its overall energy generation. Of course, the application of the stack might not allow for a drop in potential over time.
It should be noted that we considered the EOL at the step right before the intense reduction of the stack voltage and load current density, which is induced by the high gradient increase of the overpotentials. We suppose that the end behavior should not affect the comparisons mentioned above on the output powers; however, in practice, the stack would not be operated in such a regime with fast accelerating degradation.
Lower current densities under the potentiostatic operation mode lead to the reduction of the stack temperature, Fig. 9cecʺ, so that the maximum temperature is reduced about 10 C, which is opposite to the trend for the temperature evolution under the galvanostatic mode, Fig. 7fefʺ. All the slight differences shown in Fig. 8, other than those for the cathode activation overpotential and TPB density, resulting from the differences in the temperature evolutions under the potentiostatic and galvanostatic operation modes. These behaviors are discussed in detail for the effect of the inflow temperature in section Effects of the gas inflow temperatures. Fig. 10 shows the effect of the gas inflow temperatures on the evolutions of the stack performance, degradation phenomena, and related model variables. Inflow temperatures of 700, 730, and 760 C are considered for both anode and cathode sides. One can see that increasing the temperature reduces the cathode activation overpotential significantly, Fig. 10b, which means that Cr poisoning can be impeded by increasing the temperature, Fig. 10e; similar behavior is reported in the literature [18,53,56].

Effects of the gas inflow temperatures
Temperature affects the anode activation overpotential the same way, i.e. higher the temperature, lower the anode activation overpotential and TPB density, Fig. 10c and f. The reason for this reciprocal relationship between the temperature and overpotentials is that higher temperature increases the exchange current densities, Eqs. (17) and (18), which have Fig. 13 e Comparisons of the evolutions of (a) stack voltage and volume averages of (b) activation overpotential of the cathode electrode, (c) activation overpotential of the anode electrode, (d) Ni particle radius, (e) normalized TPB density of the cathode electrode, (f) normalized TPB density of the anode electrode, (g) electronic conductivity of the anode electrode, and (h) ASR of the oxide scale over the interconnect for the co-and counter-flow configurations. a reciprocal relationship with the activation overpotentials, Eqs. (15) and (16).
Higher temperatures alleviate the Cr poisoning but accelerate the rate of oxidation of the IC. Figs. 10h and i show that higher temperature increases ASR of the oxide scale over the IC and decreases the remaining Cr content of the IC, respectively. Therefore, it is concluded that the dominant degradation phenomenon shifts from Cr poisoning to IC oxidation for higher temperatures. Fig. 11 demonstrates this shift by comparing the ASRs of the cathode activation overpotential and oxide scale over the IC for different inflow temperatures. This could be seen in Fig. 10a, as the linear regime of the stack voltage drop is dominant for higher temperatures. Fig. 10d shows the effect of the inflow temperature on the Ni coarsening. It is seen that higher temperatures accelerate the growth of Ni particles, which is consistent with the trends reported in Refs. [17,20]. It should be mentioned that the value of b, which is a fitting parameter for the Ni coarsening model and used in Eq. (23), is interpolated for each temperature based on the values reported in Ref. [20] for 700, 800, and 900 C. The effect of the temperature increase for higher load current can be seen for all variables in Fig. 12 except for the cathode activation overpotential. As explained in section Spatial distribution of the degradation phenomena and related model variables, current density directly affects the cathode activation overpotential, Eq. (16). Therefore, it could be concluded that the effect of the current density counteracts the effect of the temperature on the cathode activation overpotential, as indicated in Fig. 12b. All in all, higher load currents lead to higher degradation rates for all degradation regimes, as shown in Fig. 12a.

Co-flow vs. counter-flow
Here, the long-term performance of the stack is investigated for the co-and counter-flow configurations. Fig. 13 shows that the counter-flow configuration enhances the stack lifetime, which is due to the coincidence of the maximum temperature and current density [53]. Higher temperatures and lower current densities reduce the degradation rate of the stack, as demonstrated in the previous sections. Overall, the maximums of the current density and temperature are located at the fuel inlet and air outlet, respectively. Therefore, they are coincident for the counter-flow configuration.
The maximum temperature for the counter-flow configuration is about 3 C higher than the co-flow configuration. However, a wider portion of the stack is covered with higher temperatures for the co-flow configuration so that the average temperature for the co-flow configuration is about 10 C higher than the counter-flow configuration. Fig. 13h confirms this behavior by showing a lower ASR for the oxide scale over the IC for the counter-flow configuration. It should be noted that there is a direct relationship between the temperature and ASR of the oxide scale over the IC, as demonstrated in section Effects of the gas inflow temperatures and Fig. 11.

Conclusion
In this work, a three-dimensional multiscale model of a solid oxide fuel cell full stack that can describe the degradation over thousands of hours of operation is developed. The degradation phenomena considered here are: 1) nickel particle coarsening in the anode electrode, 2) chromium poisoning of the cathode electrode, and 3) oxidation (corrosion) of the interconnect.
The electrochemical model is based on a multiscale approach, where the transport equations of mass, momentum, species, charges, and heat are all represented through effective material behaviors in a continuum (homogenization). Therefore, the operation of 38 thousand hours can be computed for the entire stack within 1 h and 15 min on a modern workstation computer. Besides the addition of the transient degradation models to the homogenization modeling approach, this work also considers how to couple the turbulent free flow in the manifolds to the laminar friction-dominated flows of the header and active domains of the stack interior.
The model is validated against the experimental data for an 18-cell Jü lich Mark-F stack for steady-state operation. Moreover, results show that the model predicts common trends reported in the literature for evolutions of stack performance, degradation phenomena, and the model variables affected by the degradation phenomena.
The developed model can thus be used to investigate the impact of various changes to the operating conditions on the different degradation phenomena and long-term performance of a full stack. The long-term performance of the stack is investigated under the galvanostatic and potentiostatic operation modes, different operating conditions (inflow temperature and load current), and flow configurations (co-and counter-flow). For the specific choice of the stack and cell characteristic, it is shown that potentiostatic operation mode, moderate temperature, lower load current, and counter-flow configuration alleviate the degradation rates and so improve the long-term performance of the stack.
The developed model can therefore be used to optimize the operating conditions for a range of operating parameters for long-term operation analysis of solid oxide cells due to its high run speed.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.