A Three-Dimensional Time-Dependent Model of the Degradation Caused by Chromium Poisoning in a Solid Oxide Fuel Cell Stack

1. Introduction

Solid oxide fuel cells (SOFCs) are a promising future power supply technology due to their high efficiency, fuel flexibility, low exhaust emissions and ability to be coupled with the heating sectors [1,2]. Currently, the mainstream material of the interconnects is steel. The metal interconnect (MIC) provides strong support and satisfying conductivity at the intermediate temperature, but it is found that the volatile chromium species CrO₂(OH)₂ and CrO₃ can be released by the oxide scale, (e.g., Cr₂O₃) formed on the surface of MIC, via the following reactions

$$2\mathrm{C}{\mathrm{r}}_2{\mathrm{O}}_3\left(\mathrm{s}\right)+3{\mathrm{O}}_2\left(\mathrm{g}\right)+4{\mathrm{H}}_2\mathrm{O}\left(\mathrm{g}\right)\rightleftarrows 4\mathrm{C}\mathrm{r}{\mathrm{O}}_2{\left(\mathrm{O}\mathrm{H}\right)}_2\left(\mathrm{g}\right)$$

(1)

$$2\mathrm{C}{\mathrm{r}}_2{\mathrm{O}}_3\left(\mathrm{s}\right)+3{\mathrm{O}}_2\left(\mathrm{g}\right)\rightleftarrows 4\mathrm{C}\mathrm{r}{\mathrm{O}}_3\left(\mathrm{g}\right)$$

(2)

The gaseous chromium species lead to a drastic degradation in the air electrode through chromium poisoning [3,4,5] that has to be mitigated [6]. Numerical modelling can be a helpful assistance in studying chromium poisoning, as it can quantitatively predict the voltage degradation under different operation conditions and stack designs. However, only limited models of chromium poisoning can be found in the current literature. A representative work comes from reference [7], which performed a one-dimensional time-dependent simulation and studied the triple-phase-boundary (TPB) and voltage evolution due to chromium poisoning in the composite air electrode (La,Sr)MnO₃/ZrO₂−Y₂O₃ (LSM/YSZ). A Butler-Volmer expression was used to predict the decrease of TPB as a result of chromium deposition. In other words, the model only considers chromium poisoning based on electrochemical reactions. The model was later extended to three dimensions [8]. The TPB was found more likely to be degraded at the interface between the air electrode and the electrolyte (AEL/EELC). A similar model was applied to study the local degradation in a 18-cells SOFC stack [9]. It was shown that chromium poisoning was the dominant degradation in the last period of operation.

Although the above models could be used to study chromium poisoning and hence simulate the voltage degradation, they were designed for LSM/YSZ air electrodes and only considered electrochemical chromium poisoning. The state-of-art air electrode consists of a mixed ionic and electronic conductor, La_0.58Sr_0.4Co_0.2Fe_0.8O_3-δ (LSCF), which exhibits a different scenario of degradation owing to chromium poisoning. For example, the formation of strontium chromate on the surface of air electrodes made by LSCF was widely observed after fuel cell operation [4,10,11]. More importantly, strontium chromate was probably generated via chemical reactions due to several findings: 1) SrO segregation on the surface of LSCF was observed at the high temperature [12,13], 2) It was found that SrO could react with gaseous chromium species to produce strontium chromate based on thermodynamic calculations [14], and 3) strontium chromate was found on the surface of LSCF under the open-circuit-voltage condition [11]. The formation of strontium chromate should not be ignored in the model as it not only lowers the conductivity, but also decreases the reactive area of the air electrode.

Another drawback in the above models is that the effect of the protective coating on the inner surface of air channels is missing. However, the protective coating exists in a real SOFC stack and an inferior protective coating surely leads to more severe chromium poisoning. For example, it has been shown that [15,16] a SOFC stack with a protective coating prepared by wet-powder-spraying (WPS) was more strongly degraded by chromium poisoning than one with a protective coating prepared by atmospheric plasma-spraying (APS), as the former had a higher porosity whereas the latter provided a fairly dense coating.

To improve the understanding of chromium poisoning in a real SOFC stack, a physical model coupling both electrochemical reactions and chemical reactions is useful. In this study, a three-dimensional computational fluid dynamics (CFD) model of a single channel of a F10 SOFC stack in Forschungszentrum Jülich GmbH (FZJ) considering chromium poisoning on the air electrode side has been developed to investigate the degradation of the F10 SOFC stack [17] performance. The voltage degradation due to decreasing reactive area during 100 kh of operation is simulated. The effects of temperature, the partial pressure of the steam in the air inlet and different protective coatings on chromium poisoning are discussed.

2. Geometry of the model

The model considers one channel of the F10 SOFC stack in FZJ. The geometry is shown in Figure 1, whose details can be found in reference [18]. On the fuel side, a Ni mesh is chosen as the gas distributor. A thick fuel electrode support layer provides the necessary mechanical strength for the cell. A thin fuel electrode is applied where the electrochemical reaction takes place. Both the fuel electrode and the support layer are made of Ni/8 mol % ZrO2−Y2O3 (8YSZ). On the air side, LSCF is chosen as the material for the air electrode. There is a contact layer (CL) between the MIC and the air electrode (AEL). It can be made of LSCF or La_0.97Mn_0.4Co_0.3Cu_0.3O_3-δ (LCC12). Unlike LSCF, LCC12 serves additionally as the chromium getter since it can absorb chromium gaseous species by forming a stable Mn-Cr spinel phase [16]. Besides, there is no strontium segregation inside LCC12.

Additionally, there should be a protective coating on the inner surface of the air channel in a real SOFC stack to prevent chromium poisoning. It is simplified as a thin surface (no volume) in the model due to the computational cost, and its influence on chromium poisoning is considered in the boundary condition for the mass transport of chromium gaseous species. Details can be found in Section 3.2.

3. Numerical model

The numerical model in this work is the extension of the previous work, where a steady-state model is developed [18,19]. By adding a degradation model of chromium poisoning, the model becomes time-dependent. Only the governing equations and parameters related to chromium poisoning are present in this work. Information about the steady-state model and parameters of the material properties, such as the conductivity and porosity of the electrodes, can be found in the reference [18].

3.1. Assumptions

There are three important assumptions for the model regarding the chromium poisoning. First, chromium poisoning is only considered on the air side. In principle, it can happen on the fuel side as well. However, it was demonstrated that chromium poisoning insignificantly affected the degradation of the fuel electrode in the fuel cell mode owing to a rather low partial pressure of chromium gaseous species [20]. Secondly, only CrO₂(OH)₂ is considered for chromium poisoning. Thermodynamic calculations pointed out that even in the case of low water partial pressure (e.g., ∼ 0.1% atm), CrO₂(OH)₂ is dominant over CrO₃ at the typical operation temperature of the SOFC stacks considered here (600 ∼ 800 °C) [3]. Thirdly, the degradation only arises from the decreasing reactive area of LSCF due to chromium deposition processes. Other factors that can also lead to degradation, such as fractures of material and lowered conductivity due to the formation of secondary phases, are currently not considered in the model.

3.2. Mass transport of CrO₂(OH)₂

The mass transport of CrO₂(OH)₂ is governed by

$$\frac{\partial \left(c{Y}_{\mathrm{CrOH}}\right)}{\partial t}+\nabla ·\left(c\overrightarrow{U}{Y}_{\mathrm{CrOH}}\right)+\nabla ·{\overrightarrow{N}}_{\mathrm{CrOH}}=Sourc{e}_{\mathrm{mole},\mathrm{CrOH}}$$

(3)

where c is the molar density of the air, $Y_{\mathrm{CrOH}}$ is the molar fraction of CrO₂(OH)₂ in the air, $\overrightarrow{U}$ is the velocity of the air that is solved from the momentum transport equation [18], ${\overrightarrow{N}}_{\mathrm{CrOH}}$ is the molar diffusion flux of CrO₂(OH)₂ and $Sourc{e}_{\mathrm{mole},\mathrm{CrOH}}$ is the source term of moles of CrO₂(OH)₂, which is related with the adsorption of CrO₂(OH)₂ on the LSCF surface and its definition can be found in Appendix A. Fick’s law [21,22] is used to solve the molar flux, given by the following equation

$${\overrightarrow{N}}_{\mathrm{CrOH}}=-D_{\mathrm{CrOH},\mathrm{air}}\nabla \left(c{Y}_{\mathrm{CrOH}}\right)$$

(4)

where $D_{\mathrm{CrOH},\mathrm{air}}$ is the diffusion coefficient of CrO₂(OH)₂ in the air.

In order to solve Equation (3), the boundary condition of $Y_{\mathrm{CrOH}}$ needs to be determined. It is assumed that CrO₂(OH)₂ only originates from the inner surface of the air channel according to the reaction (1). Followed by the reference [23,24], the reaction (1) reaches the thermal equilibrium on the inner surface of the air channel and therefore the molar fraction of CrO₂(OH)₂ at the surface of an air channel without any coatings, $Y_{\mathrm{CrOH},\mathrm{uncoated}}$, can be calculated as

$$Y_{\mathrm{CrOH},\mathrm{uncoated}}={\left(K_{\mathrm{CrOH}}·Y_{{\mathrm{H}}_2\mathrm{O}}^4·Y_{{\mathrm{O}}_2}^3\right)}^{0.25}$$

(5)

where $K_{\mathrm{CrOH}}$ is the thermal equilibrium constant of reaction (1), $Y_{{\mathrm{H}}_2\mathrm{O}}$ is the molar fraction of H₂O in the air and $Y_{{\mathrm{O}}_2}$ is the molar fraction of the O₂ in the air. According to the database summarized by Ebbinghaus [25], $K_{\mathrm{CrOH}}$ is

$$K_{\mathrm{CrOH}}={10}^{-12262/T\left[K\right]-5.2}$$

(6)

where T is the temperature. It should be noted that Equation (5) is valid for surfaces without any coatings. However, in a real FZJ SOFC stack, there is usually a protective coating on the inner surface of the air channel preventing chromium poisoning. To consider the protecting effect, the surface-averaged molar fraction, $Y_{\mathrm{CrOH}}$, is used. For each mesh grid, the surface-average molar fraction at the channel surface is calculated by including the porosity of the protective coating, ${\epsilon }_{\mathrm{coating},i}$, to Equation (5) as follows:

$$\begin{array}{cc}\hfill Y_{\mathrm{CrOH}}& =\frac{A_{\mathrm{uncoated}}{Y}_{\mathrm{CrOH},\mathrm{uncoated}}}{A}\hfill \\ \hfill & ={\epsilon }_{\mathrm{coating},i}·{\left(K_{\mathrm{CrOH}}·Y_{{\mathrm{H}}_2\mathrm{O}}^4·Y_{{\mathrm{O}}_2}^3\right)}^{0.25}\hfill \end{array}$$

(7)

where A is the surface area of the whole mesh surface, $A_{\mathrm{uncoated}}$ is the uncoated surface area and i refers to the technique for fabricating the coating.

3.3. The kinetic model of chromium poisoning

When chromium gaseous species are adsorbed on LSCF surface, the reactive surface decreases and consequently it is more difficult for oxygen being adsorbed and having the electrochemical reaction on LSCF surface. To model the process, a kinetic model that can predict the free adsorption area on LSCF is needed. This section will introduce the Langmuir model to describe the adsorption and desorption first. Then electrochemical and chemical reactions are considered for the adsorbate to calculate surface coverage of free adsorption sites.

Based on the Langmuir model [26,27], when gas molecules are adsorbed on the surface of LSCF, the change of the surface coverage can be described by

$$\frac{d{\mathsf{\Theta }}_i}{dt}=k_{\mathrm{a},i}{\mathsf{\Theta }}_{\mathrm{free}}-k_{\mathrm{d},i}{\mathsf{\Theta }}_i$$

(8)

where ${\mathsf{\Theta }}_i$ is the surface coverage of the adsorbate i, $k_{\mathrm{a},i}$ is the adsorption rate of i, ${\mathsf{\Theta }}_{\mathrm{free}}$ is the free surface coverage that can be adsorbed and $k_{\mathrm{d},i}$ is the desorption rate of i. $k_{\mathrm{a},i}$ can be obtained through Langmuir–Hertz equation [28] as follows

$$k_{\mathrm{a},i}=\frac{S_i{P}_i}{{\mathsf{\Gamma }}_{\mathrm{LSCF}}\sqrt{2\pi M_i\mathrm{R}T}}$$

(9)

where $S_i$ is the sticking coefficient of the gas i, $P_i$ is the partial pressure of the gas i, ${\mathsf{\Gamma }}_{\mathrm{LSCF}}$ is the density of the adsorption sites on LSCF particles, $M_i$ is the molar mass of the species i, $\mathrm{R}$ is the universal gas constant. $k_{\mathrm{d},i}$ generally can be calculated via [28]

$$k_{\mathrm{d},i}=k_{\mathrm{a},i}\exp \left(\frac{\Delta G_{\mathrm{ads},i}}{\mathrm{R}T}\right)$$

(10)

where $\Delta G_{\mathrm{ads},i}$ is the Gibbs free energy change in case species i being adsorbed on the LSCF surface.

In addition to adsorption and desorption between the gas and the surface of LSCF, there exists two reactions related to chromium poisoning. Figure 2 illustrates how chromium poisoning takes place inside the AEL made of LSCF. At the surface of the AEL, due to strontium segregation (see the left figure in Figure 2), the chemical reaction generating SrCrO₄ is dominating, which happens on the LSCF surface via [29]

$$\mathrm{C}\mathrm{r}{\mathrm{O}}_2{\left(\mathrm{O}\mathrm{H}\right)}_2\left(\mathrm{L}\mathrm{S}\mathrm{C}\mathrm{F}\right)+\mathrm{S}\mathrm{r}\mathrm{O}\left(\mathrm{L}\mathrm{S}\mathrm{C}\mathrm{F}\right)\stackrel{{\mathrm{k}}_{\mathrm{chemCr}}}{\rightleftarrows }\mathrm{S}\mathrm{r}\mathrm{C}\mathrm{r}{\mathrm{O}}_4\left(\mathrm{LSCF}\right)+{\mathrm{H}}_2\mathrm{O}\left(\mathrm{LSCF}\right)\Delta {\mathrm{G}}_{{\mathrm{SrCrO}}_4}$$

(11)

where $k_{\mathrm{chemCr}}$ and $\Delta G_{{\mathrm{SrCrO}}_4}$ are the rate of the surface coverage change due to the reaction and the Gibbs free energy change of the reaction (11), respectively. The method of determining $\Delta G_{{\mathrm{SrCrO}}_4}$ is given in Appendix B.

As CrO₂(OH)₂ diffuses deeper towards AEL/ELEC, the overpotential gradually decreases (to be more negative) which is likely to result in the electrochemical reaction being nonnegligible [3], as shown in the right figure in Figure 2. The electrochemical reaction is

$$\mathrm{C}\mathrm{r}{\mathrm{O}}_2{\left(\mathrm{O}\mathrm{H}\right)}_2\left(\mathrm{LSCF}\right)+3{\mathrm{e}}^{-}\stackrel{{\mathrm{k}}_{\mathrm{chemCr}}}{\to }0.5\mathrm{C}{\mathrm{r}}_2{O}_3\left(\mathrm{LSCF}\right)+1.5{\mathrm{O}}^{2-}\left(\mathrm{LSCF}\right)+{\mathrm{H}}_2\mathrm{O}\left(\mathrm{LSCF}\right)$$

(12)

where $k_{\mathrm{eleChemCr}}$ is the rate of the surface coverage change due to the reaction (12). Similar to the reference [7,8], a Bulter-Volmer type expression is employed to describe $k_{\mathrm{eleChemCr}}$:

$$k_{\mathrm{eleChemCr}}=k_{\mathrm{eleChemCr},0}{\left(P_{{\mathrm{CrO}}_2{\left(\mathrm{OH}\right)}_2}\right)}^{0.5}/{\left(P_{{\mathrm{H}}_2\mathrm{O}}\right)}^{0.5}2sinh\left(\frac{3\mathrm{F}}{2\mathrm{R}T}{\eta }_{\mathrm{air}}\right)$$

(13)

where ${\eta }_{\mathrm{air}}$ is the overpotential in the AEL. $k_{\mathrm{eleChemCr},0}$ is chosen as -5$\times {10}^{-8}$ ${\mathrm{s}}^{-1}$ in the model.

Considering all these processes and assuming the chemical reaction (11) obeys the thermal equilibrium, the kinetic model for chromium poisoning is summarized as follows:

$$\begin{array}{cc}\hfill \frac{d{\mathsf{\Theta }}_{{\mathrm{H}}_2\mathrm{O}}}{dt}=& k_{\mathrm{a},{\mathrm{H}}_2\mathrm{O}}{\mathsf{\Theta }}_{\mathrm{free}}-k_{\mathrm{d},{\mathrm{H}}_2\mathrm{O}}{\mathsf{\Theta }}_{{\mathrm{H}}_2\mathrm{O}}\hfill \\ \hfill & +k_{\mathrm{chemCr}}{\mathsf{\Theta }}_{\mathrm{CrOH}}{\mathsf{\Theta }}_{\mathrm{SrO}}+k_{\mathrm{eleChemCr}}{\mathsf{\Theta }}_{\mathrm{C}\mathrm{r}\mathrm{O}\mathrm{H}}\hfill \end{array}$$

(14)

$$\frac{d{\mathsf{\Theta }}_{\mathrm{S}\mathrm{r}\mathrm{O}}}{dt}=-k_{\mathrm{chemCr}}{\mathsf{\Theta }}_{\mathrm{CrOH}}{\mathsf{\Theta }}_{\mathrm{SrO}}$$

(15)

$$\frac{d{\mathsf{\Theta }}_{\mathrm{S}\mathrm{r}\mathrm{C}\mathrm{r}{\mathrm{O}}_4}}{dt}=k_{\mathrm{chemCr}}{\mathsf{\Theta }}_{\mathrm{CrOH}}{\mathsf{\Theta }}_{\mathrm{SrO}}$$

(16)

$$\begin{array}{cc}\hfill \frac{d{\mathsf{\Theta }}_{\mathrm{C}\mathrm{r}\mathrm{O}\mathrm{H}}}{dt}=& k_{\mathrm{a},\mathrm{C}\mathrm{r}\mathrm{O}\mathrm{H}}{\mathsf{\Theta }}_{\mathrm{free}}-k_{\mathrm{d},\mathrm{C}\mathrm{r}\mathrm{O}\mathrm{H}}{\mathsf{\Theta }}_{\mathrm{C}\mathrm{r}\mathrm{O}\mathrm{H}}\hfill \\ \hfill & -k_{\mathrm{chemCr}}{\mathsf{\Theta }}_{\mathrm{CrOH}}{\mathsf{\Theta }}_{\mathrm{SrO}}-k_{\mathrm{eleChemCr}}{\mathsf{\Theta }}_{\mathrm{C}\mathrm{r}\mathrm{O}\mathrm{H}}\hfill \end{array}$$

(17)

$$\frac{d{\mathsf{\Theta }}_{\mathrm{C}{\mathrm{r}}_2{\mathrm{O}}_3}}{dt}=0.5k_{\mathrm{eleChemCr}}{\mathsf{\Theta }}_{\mathrm{CrOH}}$$

(18)

$${\mathsf{\Theta }}_{{\mathrm{H}}_2\mathrm{O}}+{\mathsf{\Theta }}_{\mathrm{S}\mathrm{r}\mathrm{O}}+{\mathsf{\Theta }}_{\mathrm{S}\mathrm{r}\mathrm{C}\mathrm{r}{\mathrm{O}}_4}+{\mathsf{\Theta }}_{\mathrm{C}\mathrm{r}\mathrm{O}\mathrm{H}}+{\mathsf{\Theta }}_{\mathrm{C}{\mathrm{r}}_2{\mathrm{O}}_3}+{\mathsf{\Theta }}_{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}=1$$

(19)

$$K_{\mathrm{S}\mathrm{r}\mathrm{C}\mathrm{r}{\mathrm{O}}_4}=\frac{{\mathsf{\Theta }}_{\mathrm{S}\mathrm{r}\mathrm{C}\mathrm{r}{\mathrm{O}}_4}{\mathsf{\Theta }}_{{\mathrm{H}}_2\mathrm{O}}}{{\mathsf{\Theta }}_{\mathrm{SrO}}{\mathsf{\Theta }}_{\mathrm{CrOH}}}$$

(20)

where the subscript "H₂O", "SrO", "SrCrO₄", "CrOH" and "Cr₂O₃" denote the chemical species of H₂O, SrO, SrCrO₄, CrO₂(OH)₂ and Cr₂O₃, respectively. $K_{\begin{array}{c}\hfill \mathrm{SrCrO}4\end{array}}$ is the thermal equilibrium constant of reaction (11), defined as $K_{{\mathrm{SrCrO}}_4}=exp\left(-\frac{\Delta G_{{\mathrm{SrCrO}}_4}}{\mathrm{R}T}\right)$. $k_{\mathrm{chemCr}}$ can be obtained from the Equation (20) by getting the time derivatives at the both side of the equation as follows:

$$\begin{array}{cc}\hfill & \frac{d\left({\mathsf{\Theta }}_{{\mathrm{SrCrO}}_4}{\mathsf{\Theta }}_{{\mathrm{H}}_2\mathrm{O}}\right)}{dt}=K_{{\mathrm{SrCrO}}_4}\frac{d\left({\mathsf{\Theta }}_{\mathrm{CrOH}}{\mathsf{\Theta }}_{\mathrm{SrO}}\right)}{dt}\hfill \\ \hfill \to k_{\mathrm{chemCr}}=& \frac{K_{{\mathrm{SrCrO}}_4}\left[k_{\mathrm{a},\mathrm{CrOH}}{\mathsf{\Theta }}_{\mathrm{free}}{\mathsf{\Theta }}_{\mathrm{SrO}}-\left(k_{\mathrm{d},\mathrm{CrOH}}+k_{\mathrm{eleChemCr}}\right){\mathsf{\Theta }}_{\mathrm{SrO}}{\mathsf{\Theta }}_{\mathrm{CrOH}}\right]}{den}\hfill \\ \hfill & -\frac{{\mathsf{\Theta }}_{{\mathrm{SrCrO}}_4}\left[k_{\mathrm{eleChemCr}}{\mathsf{\Theta }}_{\mathrm{CrOH}}-k_{\mathrm{d},{\mathrm{H}}_2\mathrm{O}}{\mathsf{\Theta }}_{{\mathrm{H}}_2\mathrm{O}}+k_{\mathrm{a},{\mathrm{H}}_2\mathrm{O}}{\mathsf{\Theta }}_{\mathrm{free}}\right]}{den}\hfill \end{array}$$

(21)

where $den$ is

$$\begin{array}{cc}\hfill den=& {\mathsf{\Theta }}_{CrOH}{\mathsf{\Theta }}_{\mathrm{SrO}}{\mathsf{\Theta }}_{{\mathrm{H}}_2\mathrm{O}}+{\mathsf{\Theta }}_{{\mathrm{SrCrO}}_4}{\mathsf{\Theta }}_{CrOH}{\mathsf{\Theta }}_{\mathrm{SrO}}+K_{{\mathrm{SrCrO}}_4}{\left({\mathsf{\Theta }}_{\mathrm{SrO}}\right)}^2{\mathsf{\Theta }}_{CrOH}\hfill \\ \hfill & +K_{{\mathrm{SrCrO}}_4}{\mathsf{\Theta }}_{\mathrm{SrO}}{\left({\mathsf{\Theta }}_{\mathrm{CrOH}}\right)}^2\hfill \end{array}$$

(22)

Oxygen can only be adsorbed and have eletrochemical reactions on ${\mathsf{\Theta }}_{\mathrm{free}}$. It is supposed that the reaction current density in the AEL, $i_{\mathrm{v}},\mathrm{air}$, is proportional to ${\mathsf{\Theta }}_{\mathrm{free}}$[7,28], $i_{\mathrm{v}},\mathrm{air}$ can be rewritten as follows:

$$i_{\mathrm{v},\mathrm{air}}={\mathsf{\Theta }}_{\mathrm{free}}\frac{{\eta }_{\mathrm{air}}}{R_{\mathrm{air}}\exp \left(\frac{E_{\mathrm{a},\mathrm{air}}}{\mathrm{R}T}\right)a_{\mathrm{O}}{2}^{-0.25}}$$

(23)

Equation (23) is based on the definition of $i_{\mathrm{v}},\mathrm{air}$ in reference [18], where more details can be found.

Parameters used in the model can be found in Table C1.

4. Numerical setting

The simulation is carried out with the help of openFuelCell2 [19], which has been coupled into an open-source library, OpenFOAM. Discretization of the time derivative is achieved by the Backward Euler method. The IDA package in SUNDIALS [30,31] was used to solve Equations (14) ∼ (19).

Figure 3 illustrates the procedures of numerical calculation in this work:

• Set initial values of ${\mathsf{\Theta }}_{{\mathrm{H}}_2\mathrm{O}}$, ${\mathsf{\Theta }}_{\mathrm{SrO}}$, ${\mathsf{\Theta }}_{{\mathrm{SrCrO}}_4}$, ${\mathsf{\Theta }}_{\mathrm{CrOH}}$ and ${\mathsf{\Theta }}_{{\mathrm{Cr}}_2{\mathrm{O}}_3}$. Details can be found in Appendix D and Appendix E.
• Conduct the steady-state simulation based on the reference [18].
• Obtain the results of steady-state simulation [18] and set them as the initial values for the time-dependent chromium poisoning simulation.
• Set $t_{\max }$. $t_{\max }$ is equal to the time period of the steady operation of the SOFC stack.
• Conduct the chromium poisoning simulation. The chromium poisoning model is solved with other transient-state governing equations [18] of charge, momentum and mass transfer.
• Repeat step 5 until the time reaches $t_{\max }$.

Table 1 presents the governing equations of the chromium poisoning model and computational domains where corresponding governing equations are applied. It should be noted that governing equations of the chromium deposition are solely valid for the components made of LSCF.

All simulations are carried out on a single core of a computer (CentOS 9) with CPU i7-9700K and 32 GB RAM.

5. Results

5.1. Validation

When validating the model, it should be noticed that the model in this work considers the decreasing reactive area of LSCF owing to chromium poisoning as the only degradation mechanism. Table 2 presents three stacks used for validation, where stack designs (e.g., materials of CLs and fabrication techniques for protective coatings) and operation conditions are given. Three corresponding simulations are carried out, Ut-curves of which are shown in Figure 4.

F1002-97, with LCC12 as the CL and the WPS coating as the protective coating, was operated over ∼ 100 kh. Since the WPS technique can not generate a dense coating, LCC12 that can adsorb CrO₂(OH)₂ is used as the CL to reduce chromium poisoning in the AEL. However, the degradation rate of F1002-97 was still notably higher than the other two stacks, using the APS protective coating. As shown in Figure 4(a), the simulation points out that the degradation solely owing to chromium poisoning causes a degradation rate of ∼ 1.6 mV/kh. The increasing specific-area-resistance (ASR) after 100 kh operation can thus be calculated as ∼ 320 m$\mathsf{\Omega }$ ${\mathrm{cm}}^2$. Fang et al. [5] obtained the increasing ASR due to polarization in the AEL as ∼ 310 m$\mathsf{\Omega }$ ${\mathrm{cm}}^2$ by analyzing IV-curves and electrochemical impedance spectroscopy (EIS), which was close to the present simulation results. It is reasonable to compare the above two values. Chromium poisoning is believed to be the dominating degradation mechanism that increases the polarization in the AEL, because chromate contributed to the majority of the impurity in the AEL of F1002-97 [10].

F1004-67, with LSCF as the CL and APS coating as the protective coating, was operated for roughly 25 kh. The APS technique produces a suffiently dense coating that LCC12 as the CL is not needed. Thanks to the APS protective coating, the degradation rate decreased drastically. According to the simulation given in Figure 4(b), chromium poisoning is obviously suppressed as the simulated degradation rate is just ∼ 0.06 mV/kh. Experimental findings also support this conclusion. From a post-characterization, Menzler et al. [16] found a tiny amount of chromate in a similar SOFC stack with protective coating prepared by APS and operated for 36 kh. It was concluded that the application of the APS protective coating could almost solve the issue of chromium poisoning [15,16]. In addition, EIS was carried out on the stack F1004-67 with 20% humidified hydrogen at 700 °C. Figure 5 shows Nyquist plots and distribution of relaxation times (DRT) for F1004-67. From the Nyquist plots, it can be seen that degradation mainly arose from the increasing polarization. Besides, the DRT suggests that the polarization was dominated by the degradation in the fuel electrode, which according to previous research is characterized by the third peak (around 1 kHz) in the DRT plot [5,33,34]. In contrast, it is clearly found that the polarization in the AEL, represented by the second peak (around 100 Hz) in the DRT plot [5,33,34], contributed slightly to the total polarization as the peak barely changed. Therefore, a conclusion could be reached that chromium poisoning was inconspicuous in the stack F1004-67.

F1004-106 had the same stack design as F1004-67, but it was operated under ambient air for ∼ 2190 h. The validation here only considers 1500 h of operation with ambient air because a progressive degradation with unknown cause happened after the replacement of gas valves (in the last operation period of 690 h) [32]. The absolute humidity of the ambient air was calculated by using the database from [35] and the calculator from [36]. The simulation from Figure 4(c) shows a small change in the absolute values (∼ 0.44 mV/kh) of degradation rates before and after ambient air is used. A similar observation is found from the experimentally measured Ut-curves (the blue curve) in Figure 4(c). Moreover, Fang et al. [32] proposed the same conclusion by analyzing EIS data of F1004-106 (the F1 stack in reference [32]).

In brief, the model agrees well with experimental findings for different stack designs and absolute humidities in the inlet air.

Figure 5. EIS analysis of F1004-67. EIS was measured by a Zahner IM6 workstation. The upper figure is the Nyquist plot while the lower figure is the DRT plot produced by the Matlab toolbox DRTTOOLS [37].

Figure 5. EIS analysis of F1004-67. EIS was measured by a Zahner IM6 workstation. The upper figure is the Nyquist plot while the lower figure is the DRT plot produced by the Matlab toolbox DRTTOOLS [37].

5.2. Simulation results of F1002-97

Simulation results of F1002-97 are present here to illustrate chromium poisoning. Figure 6 gives the spatial distribution of physical variables inside the AEL. Figure 7 presents the evolution of average $k_{\mathrm{chemCr}}$, $k_{\mathrm{eleChemCr}}$ and ${\mathsf{\Theta }}_i$ inside the AEL over 100 kh.

According to Figure 6(a), ${\mathsf{\Theta }}_{{\mathrm{Cr}}_2{\mathrm{O}}_3}$ barely exists inside the AEL before 25 kh. The same result can be more clearly found in Figure 7(a), where average ${\mathsf{\Theta }}_{{\mathrm{Cr}}_2{\mathrm{O}}_3}$ is close to zero. This is because at the start of the operation, ${\eta }_{\mathrm{air}}$ is not high as there is still sufficient free surface, ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{free}\end{array}}$, for oxygen reduction reaction. As time goes by, ${\eta }_{\mathrm{air}}$ becomes more negative due to chromium deposition and hence $k_{\mathrm{eleChemCr}}$ increases exponentially as presented in Figure 7(b), which can be explained by using a Butler-Volmer type expression for $k_{\mathrm{eleChemCr}}$. At the end of the simulation, the highest ${\mathsf{\Theta }}_{{\mathrm{Cr}}_2{\mathrm{O}}_3}$ is found at the AEL/EELC where the lowest ${\eta }_{\mathrm{air}}$ locates as shown in Figure 6(d). From Figure 7(a), the average ${\mathsf{\Theta }}_{{\mathrm{Cr}}_2{\mathrm{O}}_3}$ after 100 kh in the AEL is around 0.03, which is small. However, the segregation of Cr₂O₃ on AEL/ELEC, which leads to ∼ 0.17 ${\mathsf{\Theta }}_{{\mathrm{Cr}}_2{\mathrm{O}}_3}$ locally, may still cause a nonnegligible ohmic resistance since a low electrical conducting layer is likely to form.

The distribution of ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{SrCrO}4\end{array}}$ is shown in Figure 6(b). Unlike ${\mathsf{\Theta }}_{{\mathrm{Cr}}_2{\mathrm{O}}_3}$, ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{SrCrO}4\end{array}}$ tends to segregate on the surface of the AEL, which is consistent with the experimental observation from reference [4,10,11]. Plus, in opposition to the evolution trend of ${\mathsf{\Theta }}_{{\mathrm{Cr}}_2{\mathrm{O}}_3}$, ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{SrCrO}4\end{array}}$ increases quickly within the first 25 kh while it hardly changes after 60 kh. The trend is more distinctly illustrated by the evolution of the average ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{SrCrO}4\end{array}}$ shown in Figure 7(a). This can be explained by the assumption that reaction (11) obeys the thermal equilibrium. When the operation starts, a high $k_{\mathrm{chemCr}}$ should be observed as a tiny amount of ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{SrCrO}4\end{array}}$ and ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{H}2\mathrm{O}\end{array}}$ exists on the LSCF surface. With increasing ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{SrCrO}4\end{array}}$ and ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{H}2\mathrm{O}\end{array}}$, due to $k_{\mathrm{chemCr}}$ being constrained by the thermal equilibrium, a fast attenuation of $k_{\mathrm{chemCr}}$ happens, as shown in Figure 7(b).

The spatial evolution of ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{free}\end{array}}$ is given in Figure 6(c). Similar to ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{SrCrO}4\end{array}}$, ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{free}\end{array}}$ decreases quickly at the beginning of 25 kh, probably due to the fast adsorption of CrO₂(OH)₂ at the initial stage as shown in Figure 7(a). Consequently, ${\eta }_{\mathrm{air}}$ decreases as governed by Equation (23). In addition, the simulation shows the highest ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{free}\end{array}}$ is found above the rib and close to AEL/ELEC. Advanced surface characterization is needed to further validate this distribution. At 100 kh, as a result of a rather low ${\mathsf{\Theta }}_{\begin{array}{c}\hfill \mathrm{free}\end{array}}$, ${\eta }_{\mathrm{air}}$ roughly decreases by a factor of four in comparison with the initial state, as shown in Figure 6(d). At the air inlet and air outlet, ${\eta }_{\mathrm{air}}$ is not much influenced by the chromium poisoning. The probable explanation is that electrochemical reactions mainly take place at the centre of the air electrode in the situation where the counter-flow field for the fuel and air is applied.

Figure 6. Spatial distribution of physical variables in the AEL of F1002-97 (a) ${\mathsf{\Theta }}_{{\mathrm{Cr}}_2{\mathrm{O}}_3}$, (b) ${\mathsf{\Theta }}_{{\mathrm{SrCrO}}_4}$, (c) ${\mathsf{\Theta }}_{\mathrm{free}}$ and (d) ${\eta }_{\mathrm{air}}$. Slices are in the y-z plane and uniformly distributed across the whole AEL. The first and the last slices are at the boundary faces of the geometry. The geometry is scaled according to (x,y,z) = (0.1, 1, 10).

Figure 6. Spatial distribution of physical variables in the AEL of F1002-97 (a) ${\mathsf{\Theta }}_{{\mathrm{Cr}}_2{\mathrm{O}}_3}$, (b) ${\mathsf{\Theta }}_{{\mathrm{SrCrO}}_4}$, (c) ${\mathsf{\Theta }}_{\mathrm{free}}$ and (d) ${\eta }_{\mathrm{air}}$. Slices are in the y-z plane and uniformly distributed across the whole AEL. The first and the last slices are at the boundary faces of the geometry. The geometry is scaled according to (x,y,z) = (0.1, 1, 10).

Figure 7. Evolution of average (a) ${\mathsf{\Theta }}_i$ and (b) $k_{\mathrm{chemCr}}$ and $k_{\mathrm{eleChemCr}}$ in the AEL of the stack F1002-97.

Figure 7. Evolution of average (a) ${\mathsf{\Theta }}_i$ and (b) $k_{\mathrm{chemCr}}$ and $k_{\mathrm{eleChemCr}}$ in the AEL of the stack F1002-97.

5.3. Effects of temperatures and humidities on chromium poisoning

Figure 4(b) and 4(c) demonstrate that with the help of APS protective coatings, SOFC stacks suffer little from chromium poisoning. It is thus significant to see whether the above stack designs can also efficiently prevent chromium poisoning at lower temperatures and ambient air.

Figure 8(a) studies the influence of temperature on chromium poisoning. A lower temperature leads to a smaller degradation rate owing to slower adsorption kinetics. However, it also leads to a lower voltage because conductivity decreases. To maintain the required power output, the active area needs to be larger and the number of layers in stacks should be increased. Unfortunately, both strategies lead to a challenge of the mechanical strength of the stack. Economic analysis, system-level simulation and solid mechanics modelling are needed to further demonstrate the feasibility of lowering the operation temperature.

Figure 8(b) compares Ut-curves under different absolute humidities in the air at 650 °C. The aim is to determine whether it is possible to save energy for the dehumidification of the air by maintaining an acceptable degradation rate. First of all, 2 % humidity is undoubtedly not feasible because the degradation rate due to chromium poisoning is already ∼ 0.14 %/kh. The total degradation rate in an actual SOFC stack is likely to be above the target, 0.3 %/kh [38], if other degradation processes are involved. Secondly, 1.5 % humidity is also not recommended as there is a tendency for progressive degradation after 80 kh. A humidity of around 1% could be more promising, as less energy is required for dehumidification of the air, no progressive degradation is observed and the degradation rate is low even though it is ten times higher compared to the case of 0.1 % humidity.

6. Discussion

6.1. One channel model representing the chromium poisoning in the whole SOFC stack

Due to computational cost, this work only considers the single channel of a SOFC stack to study chromium poisoning. Previously, it has been proven that this simplification is already sufficient to provide accurate predictions of IV-curves [18]. For the simulation of chromium poisoning, it is worthwhile to discuss further whether the simplification is acceptable to provide reliable results. Three variables are discussed here, temperature, pressure and the origin of CrO₂(OH)₂ since they directly influence governing equations of chromium poisoning. The SOFC stacks considered in this work were tested inside ovens. Thanks to the heat radiation, based on the temperature measured from thermocouples in the stack, it was found that the spatial variation of temperature across the whole stack during the steady operation was slight (10 ∼ 15 °C). Hence, the temperature field in the single channel is expected to be similar to that in the SOFC stack (reactive parts). Unlike temperature, the spatial distribution of pressure is difficult to measure, and usually it is obtained through numerical calculations. The simulation on the single cell stack (i.e., the geometry in reality) of F1004-67 found the range of pressure on the air side is 100 kPa ∼ 100.473 kPa [39], indicating the pressure is uniformly distributed. As for the origin of CrO₂(OH)₂, the model in this work assumes that CrO₂(OH)₂ solely comes from the internal surface of the channel via reaction (1). In principle reaction (1) could also occur on the metal frames’ surface. However, internal experiments indicate the impact of metal frames is neglectable for at least 3 kh.

In brief, although there could be some deviations, the accuracy of the single channel model on chromium poisoning is likely acceptable. Hence, its simulation results can be compared with experimental findings on the stack level.

6.2. The stack design

In this work, two stack designs are compared regarding the performance degradation due to chromium poisoning. The first employs the APS protective coating and LSCF as the CL (F1004-67). As a result of a quite dense protective coating, chromium poisoning is almost neglectable. Degradation in the stack F1004-67 arose from other processes according to post-test analysis, such as gas leakage, nickel agglomeration, manganese diffusion, increasing ohmic resistance and contact loss due to the formation of secondary phases. The second utilizes a protective coating prepared by WPS and LCC12 as the CL (F1002-97). Since the WPS protective coating has a relatively high porosity, LCC12 that can adsorb CrO₂(OH)₂ was chosen as the CL. Although the adsorption rate of CrO₂(OH)₂ by LCC12 is two orders of magnitude higher than that of LSCF (see definitions of $Sourc{e}_{\mathrm{mole},\mathrm{CrOH}}^{\mathrm{LCC}12}$ and $Sourc{e}_{\mathrm{mole},\mathrm{CrOH}}^{\mathrm{LSCF}}$ in Appendix A), this stack design is unable to prevent chromium poisoning efficiently, as shown in Figure 4(a). However, WPS technology is cheaper than APS technology. If WPS protective coatings need to be applied after considering the costs, fabrication should be optimized to lower its porosity.

6.3. Operation under higher air absolute humidity and lower temperature

As shown in Figure 8, based on the current model and only from the perspective of the degradation rate, there is a potential to lower the operation temperature to 650 °C and increase absolute air humidity to 1 %, if the APS protective coating is applied. However, it remains uncertain whether it is worthwhile to make such changes in future. The present model may underestimate the degradation due to chromium poisoning as it ignores the influence of the insulation phase SrCrO₄. Figure 9 shows the distribution of ${\mathsf{\Theta }}_{{\mathrm{SrCrO}}_4}$ inside the AEL and the CL under the above-mentioned operation conditions and stack designs. After 100 kh operation, a clear gathering of ${\mathsf{\Theta }}_{{\mathrm{SrCrO}}_4}$ on the surface of the CL can be found. ∼ 0.28 surface coverage of the insulation phase could lead to substantial ohmic resistance. A more sophisticated model is therefore needed to understand the increase of ohmic resistance.

7. Conclusion

A three-dimensional CFD model of a single channel in a Jülich F10 SOFC stack is proposed, taking into account one of the most important degradation process, chromium poisoning. The numerical simulation is capable of predicting voltage degradation during 100 kh of operating time. The model is able to predict the reduction of free adsorption sites on LSCF during chromium poisoning and the polarization degradation due to that reduction. The model is supported by comparing simulation results with experimentally measured voltage degradation and EIS data. The model can be used to study chromium poisoning under different temperatures, absolute humidity and stack designs. Based on the simulation results, three important conclusions are drawn:

• With APS protective coating, chromium poisoning in the SOFC stack is almost solved.
• Lower temperatures and less moisture in the inlet air mitigate chromium poisoning.
• It should be possible to operate a SOFC stack with 1 % humidified air at 650 °C, if APS protective coating is applied.

As a result of the segregation of ${\mathsf{\Theta }}_{{\mathrm{SrCrO}}_4}$ and ${\mathsf{\Theta }}_{{\mathrm{SrCrO}}_4}$, an increasing ohmic resistance due to chromium poisoning is expected. However, the influence of the decrease in conductivity is not considered in this model. Incorporating the decreasing conductivity into the model will produce a more reliable prediction of the performance of SOFC stacks and will be considered in future. In addition, some kinetic parameters are assumed in the present model, such as sticking coefficients and the prefactor of $k_{\mathrm{eleChemCr}}$. These parameters need to be experimentally measured or theoretically calculated for the numerical model to provide a more precise prediction of the degradation of chromium poisoning.