Solid Oxide Fuel Cell Performance Analysis through Local Modelling

Solid Oxide Fuel Cells (SOFC) are an emerging technology among different fuel cell types since they are successfully used in stationary cogeneration units to produce heat and electricity. Different scale applications are proposed as alternative energy sources for residential usage and industrial power plants, reducing the greenhouse gas emissions which characterize fossil-fuel-based processes. Their spread is favoured by the development of proper simulation tools that allow system design optimization and control in real-time operations. For this purpose, model building and validation, through comparison with experimental observations, are fundamental steps to guarantee the simulation validity. A single-anode-supported planar SOFC with two possible cathodic current collector designs is tested in common operating conditions, evaluating the performance through EIS analysis and characteristic curves. These provide a preliminary validation for the proposed 2D steady state simulation code. This model, implemented in Fortran, makes it possible to forecast the main SOFC local properties on both the anodic and cathodic sides. The key point of the code is the electrochemical kinetics, based on a semi-empirical approach where requested parameters, derived from fitting of experimental results, are introduced in physically based equations. In this way, the influence of specific cell design on system performance is evaluated.


Introduction
Solid Oxide Fuel Cells (SOFCs) are promising technologies for more eco-friendly energy production with minimum polluting emissions. They are based on the direct conversion from reactant chemical energy to electricity through an electrochemical process. High operating temperatures guarantee optimized conditions in terms of cost and obtained performance, favouring the system kinetics. Differently from low-temperature fuel cells [1], ceramic oxides substitute expensive noble metals as electrocatalysts. Moreover, they are more resistant to pollutants, so that, for example, a syngas, containing not only H2 but also CO and CO2, can be used without catalyst poisoning [2]. In particular, the electrode materials have to include several properties: good catalytic activity (for H2 oxidation and O2 reduction at the anodic and cathodic sides, respectively); high electronic and ionic conductivity making it possible to maximise the electrode volume use; and mechanical and chemical compatibility with other cell components [3]. This results in selective, resistant and long-lasting catalysts, which provide minimum overpotential with a high number of Triple Phase Boundary (TPB) active sites, where redox reactions occur due to the co-presence of gas, electrons and ions [4]. Meanwhile, the electrolyte should be a thick material with high ionic conductivity and negligible internal electron transport. In view of these features, the common industrial solution is a layer of doping fluorite-structure ceramics as electrolyte, between a porous Ni-based cermet as the anode and LaSr perovskite as the cathode [5]. Ni particles are deposited on a solid oxide bed, guaranteeing good activity for hydrogen dissociation and wide active area. Moreover, the metal-ceramic combination favours a considerable electrical and ionic conductivity [6]. Meanwhile, for the cathode, perovskites are well-known catalysts for the oxygen reduction mechanism, adding the advantages of a Mixed Ionic Electronic Conductor (MIEC) [7]. High working temperature not only makes it possible to optimise the materials used, but also to obtain a more efficient system for energy production. Indeed, SOFC is a cogeneration unit, which can provide both electricity and heat from exhaust gases [8][9][10]. Nevertheless, still this technology is not fully ready to be competitive on the market. An incentive to SOFC spread at industrial level is the development of specific models, which simulate both electrochemical kinetics and overall cell performance. Their application makes it possible to monitor the system, to predict and identify the abnormal plant situations. Indeed, knowing process inputs, real-time quality information is provided. Simulation is also used for sensor validation: namely, fault detection through the continuous comparison between measured and modelling values, and for interferential measures, when direct laboratory tests are not possible or too expensive. In the literature, different models are proposed, more or less detailed in view of their specific applications. In a 0D approach, the cell is reduced to a point, so it allows a quick performance evaluation with a minor computational effort. It is used for the complete power plant simulation, where SOFCs are paired with other unit operations [11]. A 1D model studies system property changes along the main flow direction [12], whereas 2D or 3D tools consider gradients on cell planes and in all spatial directions respectively [13][14][15]. They provide local behaviour information, which could be difficult to measure through experimental tests. It is possible to achieve punctual control of cell operation, which allows the detection of hotspots, channelling effects or irregular use of the active surface. Still, in this case, the simulation resolution is very complex and time consuming. The model reliability is strictly connected to proposed electrochemical kinetics, which depends on specific SOFC properties, such as electrocatalytic materials, geometrical design, electrode micro-structure, used interconnects, and so on [16][17][18]. Careful model building and tuning are fundamental to guarantee the empirical agreement. For this purpose, the simulation should be based on both theoretical remarks and experimental comparison. Since the crucial point is the identification of different overpotentials, specific electrochemical tests are performed to detect the main system dependences. A first evaluation of cell behaviour is provided by characteristic curves, which show the current density influence on obtained voltage. Increasing the electric load, there is a gradual tension drop due to different polarization losses, which depend on intrinsic material features and system feed conditions, such as temperature and gas composition [19]. For each specific I-V pair, the ohmic, activation and concentration overpotential terms are estimated through Electrochemical Impedance Spectroscopy analysis (EIS). The EIS measurements are fitted through Equivalent Electric Circuit models (EEC), where every element is representative of a physicochemical process occurring in the SOFC. They give quantitative information about mass transport, reaction rates and the resistances of each cell component [20]. However, due to the distribution and coupling of electrochemical phenomena, EEC models may be complex and strongly influenced by the composition, morphology and combination of cell materials. Therefore, no unique modelling is proposed in literature, but rather several solutions depending on the cell type and testing conditions [21][22][23].
In this context, the study proposes a 2D model of a planar H2-feed SOFC that focuses on evaluation of the main operating parameters on both the anodic and cathodic sides. Thus, this simulation gives local information on cell behaviour and can suggest optimized working conditions depending on the provided inputs. Maps, which show the distribution of gas compositions, temperature, current density and polarization losses, are detected on cell planes. The electrochemical kinetics is based on a semi-empirical approach. Starting from a theoretical formulation, requested parameters derive from EIS analysis and I-V curves in order to consider the specific characteristics of the studied system. The electrochemical process is schematized through equivalent electric circuit to underline different polarization contributions. In view of this, an experimental campaign is performed on a single anode-supported SOFC. Two solutions for the cathodic current collectors are tested, focusing on different occurring overpotentials and how to correct the inhomogeneous current distribution. In particular, the influence of contacting layer between Pt-mesh and electrode is evaluated. The balance of plant and SOFC design are minimized, aiming at the study of performance due to only electrochemical processes avoiding the influence of other cell components.

The SIMFC Code and its Kinetic Core
Previously validated for Molten Carbon Fuel Cells (MCFCs) [24], the SIMulation Fuel Cell (SIMFC) is in-house developed Fortran code, which allows a quick forecast of planar cell behaviour. It consists of a 2D local model, based on material, energy, momentum and charge balance resolution in steady state condition, assuming to have H2-rich anodic feed. The basic formula are reported in [25]. The differential equations are solved introducing the finite difference approach to simplify problem in a system of ordinary equations. Results are obtained assuming to reach the same potential in every cell point. The performances are evaluated on different cell planes, which are divided in an optimised number of sub-elements, where balances are applied. Since each variable is calculated at specific local conditions, a detailed description of transport and reaction mechanisms occurring inside the cell is obtained. Evaluation of I-V curves and maps on both anodic and cathodic sides of FC chemical-physical properties are the main results.
The core of the code is the electrochemistry kinetics, which has the major influence on system performance. The work focuses on the simulation of an anionic conductive-electrolyte SOFC, where fed O2 is reduced at cathode (Equation 1), producing O 2-ion flow to anode, where H2 oxidation occurs (Equation 2). The global reaction is shown by Equation (3).
The cell performance is evaluated considering the Open Circuit Voltage (OCV), the maximum theoretical obtainable value, and subtracting the different polarization losses (Equation 4).
The first term derives from the Nernst equilibrium law, and it is correct considering possible occurring leakage phenomena (Equation 5) [26]: where the reversible voltage E 0 is function of temperature according to Equation (6) [27]: Different polarization contributions are evaluated by a semi-empirical approach: the specific cell parameters are obtained through experimental data fitting.
The ohmic overpotential, due to current transfer resistances, depends on cell materials and temperature, assuming a thermally activated mechanism for the charge migration (Equation 7) [25].
The activation polarization is related to the electrochemical process; indeed, a surplus of energy is required to overcome the energy barrier before the reaction begins. This is derived from the Butler-Volmer formulation [27], neglecting the concentration gradients between bulk and active sites since the activation loss is relevant at low current load above all. Thus, the hyperbolic sine equation is introduced for both electrodes, guaranteeing a more precise formulation (Equation 8) [28].
The exchange current density, J0, the theoretical value at equilibrium conditions, depends on both temperature and gas composition. Assuming work at atmospheric pressure, Equations (9) and (10) are used [29].
The concentration overpotential considers the resistance due to gas transport: when high currents are applied, the diffusion of reactants and products from and to electrochemical reaction sites is too slow to maintain the initial bulk composition at active sites and so a drastic drop of cell performance occurs. If air is fed as oxidant, the cathodic concentration loss can be neglected [30]. The anodic term is derived from the Butler-Volmer equation, considering that the exchange current density diverges in this case. 1D material balance is developed along the electrode thickness to detect the gas profile (Equation 11; for the complete explication, refer to [31]).
The diffusion coefficients are the combination of molecular contributions, consisting of interactions among different particles, and the Knudsen term due to molecule-wall collisions. The first is evaluated through the Chapman-Enskog approach for a binary mixture (Equation 12), corrected in order to consider the presence of more compounds (Equation 13) [32].
Meanwhile, the Knudsen term is equal to (Equation 14) [33]: Considering two series mechanisms, the effective diffusion coefficient is calculated with Equation (15).

The Reference Equivalent Circuit
Experimental data are fitted through the Equivalent Electric Circuit (EEC) reported in Figure 1. Such a configuration was used in a previous study to model the same types of anode-supported cells [34]. An inductive element (L) is used to fit the respective high-frequency contribution, commonly due to connections between the sample and the instrument. The second circuit element (R-s) represents the series resistance, depending on the electrolyte and current collector/electrode contacting losses. Two R//CPE meshes (R1-anode//CPE1-anode and R2-anode//CPE2-anode, where CPE stands for Constant Phase Element) are added to fit the high-frequency contributions of the cell polarization. Such elements describe gas diffusion, charge transfer and ionic transport processes occurring at the anode functional layer. The medium-frequency range is fitted through a Gerischer element, which the resistive cathodic contribution is calculated from (R-cathode). This takes into account the surface exchange kinetics of oxygen and its ionic diffusivity within the bulk of the electrode. The low frequency part of the cell impedance is described using a bounded Warburg element, modelling gas diffusion in the anode substrate (R3-anode).

EIS Analysis and Kinetics Parameter Identification
EIS analysis is an accurate technique, which allows a detailed knowledge of SOFC performance; indeed, the major loss sources are distinguished through the deconvolution method reported in Section 4.3. For every working condition, the Nyquist plot (Zre vs. Zim, where Zre and Zim are the real and imaginary parts of the impedance, respectively) can be evaluated. As Figure 2 shows, it reports the impedance recorded at OCV at different temperatures and gas compositions. All the spectra are characterized by an inductive contribution, visible in the negative part of each plot, in the highest frequency range (i.e., above 10 kHz). For lower frequencies, spectra evolve in the positive part of the Nyquist plots, indicating the presence of resistive and capacitive effects. From this portion, it is possible to identify the R-s contribution of the EEC as the high-frequency intercept of the real part axis. Furthermore, the low-frequency intercept is the total cell polarization (R-tot); the sum of ohmic, activation and diffusion terms, as described in Section 2.2. Thus, according to the literature [30,35], open circuit impedance analysis allows the determination of model parameters. Changing the working conditions, their influence on electrochemical kinetics is evaluated. At first, EIS data are studied in the case of a Pt-mesh+LSCF cell, where the contacting paste guarantees a homogeneous current distribution on the cell cathodic plane, thus reducing the resistive losses attributable to R-s at the tested temperatures: 1073 K and 1023 K (Figure 2a,b respectively). In particular, impedance measured at 1023 K exhibits higher values of both R-s and R-tot due to thermal dependency of the material conductivity and catalytic activity, as considered in the proposed model. The EIS data for the Pt-mesh+LSCF cell are analysed to detect the values P1-P2 of ohmic resistance and the anodic reaction activation energy Eatt,an. Increasing the amount of oxygen at the cathode side from 10% to 100% results in a shrinkage of the spectra due to the reduction of the mid-and lowfrequency contributions. The same effect is underlined changing the H2 percentage from 47% to 97%. The effects of gas composition on the impedance are investigated in detail for the cell Pt-mesh+LSCF, as displayed in Figure 2c,d for the cathode and anode sides, respectively. In both cases, the increase in reactant concentration (i.e., oxygen at the cathode, hydrogen at the anode) favours the electrochemical kinetics, reducing the cell polarization without affecting the R-s contribution. These measurements permit the estimation of anodic and cathodic kinetic orders A and C, whereas H2O dependence is not directly evaluated, since tests are performed without modifying this variable. The pre-exponential coefficients of activation polarization, γan and γcat, are derived from characteristic curves as well as the leakage loss. Other required parameters are retrieved from the literature [25]. Then, these values are adapted for a Pt-mesh cell, considering that the cathodic cell plane is not actually isopotential [34]. As Figure 2a,b show, the internal resistance increases in this second case, since the absence of cathodic contacting paste disadvantages electrical conduction. The cathodic process is also found to be influenced by an irregular current distribution on the cell plane. Thus, some different kinetic parameters are detected in the Pt-mesh cell case. Table 1 reports the proposed electrochemical kinetic values for both tested cells. These parameters are in accordance with previous literature studies, where the activation energy ranges between 100-120 kJ/mol for anodic and 110-160 kJ/mol for cathodic processes, respectively [36]. The O2 kinetic order is commonly considered equal to 0.25 since the cathodic rate-determining step is identified as the oxygen ion formation and incorporation into the electrolyte [37]. In contrast, there are no well defined values for H2 and H2O dependences. Parameter A varies from 0.1 to 2 and B from −0.5 to 1 [29]. The pre-exponential coefficients are specific to the electrocatalytic materials used (i.e., Ni/8YSZ-cermet as anode and LSCF as cathode), yet the detected values are comparable with the range in the literature (γan ≅ 5.5 × 10 8 -5.5 × 10 10 A/m 2 , γcat ≅ 7 × 10 8 -7 × 10 10 A/m 2 ) [30].
The requested parameters for concentration overpotential are presented in previous works [31,33], while the micro-structural cell features are derived from post-experimental characterization of different layers (refer to Section 4.4).

SOFC Local Performance Evaluation
In the first step, characteristic curves are analysed to evaluate the system's behaviour in the tested conditions (refer to Section 4). As expected, the use of only Pt-mesh as a current collector causes worse performance at the same operating conditions (Figure 3). A proper tuning of model parameters, which considers this effect, allows the forecast of both configurations. The results show a good match between simulated and experimental values of I-V curves; the error is always lower than 2%. This confirms the validity of the introduced assumptions. When the overall behaviour is simulated, the cell can be considered as an isopotential plane for both anodic and cathodic sides. Indeed, the worst performance detected in the second case study can be modelled only by adjusting kinetic parameters, without a major modelling inaccuracy. Thus, at this simulation level, a more detailed approach would only involve greater computational effort. However, the characteristic curves do not permit local control of the system. Thus, the simulation provides further information compared to experimental results through the maps of the main chemical-physical cell properties. Considering tested co-flow configurations in common operating conditions, the distributions of gas composition and current density are obtained ( Figure  4). As expected, the electrochemical process develops along the flow direction: the highest conversion is obtained in the first part of the cell, where the current density is major. Therefore, the possible fuel starvation in the last cell sections can be easily evaluated through these maps, which makes it easier to choose proper anodic and cathodic flow rates to avoid a probable cause of degradation [38]. Another issue is the development of hotspots, which can cause cell damage due to the exothermic electrochemical process. In view of this, the proposed code also makes it possible to evaluate the temperature distribution on both sides ( Figure 5). In this case, the gradient is around ten degrees, given the low cell dimension, yet it is no more negligible at the industrial scale, where the possibility to forecast temperature maps can prevent the system failure. It is noteworthy that the difference in temperature gradient changes with the cathodic current collector. The higher obtained limit in the Pt-mesh+LSCF cell is 1033 K, whereas, when only Pt-mesh is present, a 1037 K value is reached, above all due to ohmic resistance increase ( Figure 6). In the first case, it is about 0.13 Ωcm 2 , and it rises to 0.38 Ωcm 2 in the second cell type.

Anode-Supported Cells and Current Collectors
Commercially available Anode-Supported Cells (ASCs) were purchased from CeramTec AG (Marktredwitz, Germany). Cells consisted of a 5 × 5 cm 2 Ni/8YSZ-cermet anode (where 8YSZ: 8 mol % Y2O3 doped ZrO2), a 5 × 5 cm 2 8YSZ electrolyte with a GDC (Ce0.8Gd0.2O2-δ) interlayer at the cathode side followed by a LSCF (La0.6Sr0.4Co0.2Fe0.8O3-δ) 4 × 4 cm 2 cathode (Figure 7). Metallic meshes (mesh width of 0.12 mm, 3600 meshes per cm 2 , Heraeus Deutschland GmbH & Co. KG) were used as current collectors between cell electrodes and the electronic characterization device. A platinum mesh was used at the cathode side, with a nickel one at the anode to avoid alloying with the nickel contained in the electrode. The current collectors were spot-welded, with the platinum wires used as potential probes and current carriers in a 4-wire circuit configuration.

LSCF Contacting Paste
Cells were tested with and without LSCF contacting paste at the cathode side. The use of a contacting paste aims at improving cell performance by increasing the contact area between two connected elements [39]. It was prepared by mixing a commercially available LSCF powder (FCM, USA) with a solvent blend (94 wt.% terpineol, 6 wt.% polyethylene glycol) to form an ink (60 wt.% powder, 40 wt.% solvent). This mixture was homogenized in an Exakt 80E three-roll mill (Exakt Advanced Technologies GmbH, Norderstedt, Germany) and brushed over the cathode to obtain a homogeneous layer. The platinum mesh was positioned immediately after the contacting paste application in order to be embedded, extending the contact area.

Cell Testing
Cells were tested using the setup configuration described in [40]. Each cell was heated up to 1173 K with a ramp speed of 5 K/min in flowing dry air at the cathode side (0.5 slpm/min) and dry nitrogen at the anode (0.5 slpm/min). At this temperature, the anode gas flow composition was stepwise turned to 1 slpm/min of 3% wet hydrogen to perform the electrode reduction, while the air flow at the cathode was increased to 1 slpm/min. After one hour in such conditions, the furnace temperature was changed to 1073 K and 1023 K to evaluate the cell performance by I-V curves and EIS analysis. Moreover, some specific EIS measurements were set in OCV conditions, to underline the influence of anodic and cathodic compositions on electrochemical kinetics. H2 molar composition was set at 97%, 67% and 47%, by diluting with N2 and fixing 3% of moisture, whereas O2 was set at 10%, 21%, 50% and 100%.
The electrochemical characterization was carried out using a Solartron 1287 potentiostat/galvanostat coupled with a Solartron 1260 frequency response analyzer (Solartron Analytical, Farnborough, UK). I-V curves were measured by applying an increasing current load by steps of 0.125 mA/cm 2 up to cell voltages not below 0.6 V. EIS spectra were recorded applying a 10 mV AC signal in the frequency range 0.1 Hz-100 kHz and deconvoluted through a Complex Non-Linear Regression Least-Squares (CNRLS) method [41][42][43][44][45]. Fitting was carried out using the software ZView® (Scribner Associates Inc., Southern Pines, USA), accepting an error below 1% and a relative uncertainty for each element lower than 30%.

Post-Experimental Characterization
Cells were studied post-experimentally, preparing their cross-sections. This procedure consisted in mounting in epoxy resin, cutting and polishing with diamond suspensions up to 250 nm grain size. Sections were observed and analyzed using a Scanning Electron Microscope (SEM) Zeiss Evo 40 (Carl Zeiss AG, Oberkochen, Germany), equipped with an EDX detector PentaFET (Oxford Instruments, High Wycombe, UK) . Micro-structural parameters (i.e., the porosity and thickness of each layer) were evaluated through image analysis (Fiji Imagej 1.47 h) on SEM micrographs. Figure 7 shows the cross section of a cell used in this study. The cathode is characterized by inhomogeneous porosity and irregular top surface, while the Ni/8YSZ-cermet anode exhibits a porosity gradient decreasing towards the electrolyte. Table 2 reports the detected parameters valid for both configurations.

Conclusions
SOFC behaviour simulation is a useful tool for both the proper system design and the control during real-time operation. Data at the local scale are essential to detect possible abnormal working conditions and so to forecast the beginning of degradation mechanisms. In this context, a 2D model for planar SOFC simulation is developed to evaluate the main chemical-physical properties, such as gas composition, temperature and current density. A single anode-supported SOFC is tested in common operating conditions, evaluating the performance through EIS analysis and I-V curves. As the developed equivalent electric circuit shows, temperature, H2 and O2 inlet compositions have the main influences on polarization losses. Two different configurations are considering, adding Pt-mesh with and without contacting paste as cathodic current collector. In the first case, a higher electric conductivity is obtained, and the cell has a better performance. In contrast, the use of only Pt-mesh provokes a drastic voltage drop, due to an irregular current distribution. Different ohmic and cathodic polarization contributions are introduced in the SOFC simulation to consider this aspect. The good agreement between experimental and simulated data validates proposed modelling assumptions. Thus, isopotential planes can be supposed without a relevant error if the global cell performance has to be evaluated. Meanwhile, the inhomogeneous current load can be no more negligible in a microscale modelling approach, where local behaviour at the electrode-current collector interface is studied.