Thin Solid Film Electrolyte and Its Impact on Electrode Polarization in Solid Oxide Fuel Cells Studied by Three-Dimensional Microstructure-Scale Numerical Simulation

: In this work, a three-dimensional microstructure-scale model of a Solid Oxide Fuel Cell’s Positive-Electrolyte-Negative assembly is applied for the purpose of investigating the impact of decreasing the electrolyte thickness on the magnitude, and the composition of electrochemical losses generated within the cell. Focused-Ion-Beam Scanning Electron Microscopy reconstructions are used to construct a computational domain, in which charge transport equations are solved. Butler–Volmer model is used to compute local reaction rates, and empirical relationships are used to obtain local conductivities. The results point towards three-dimensional nature of transport phenomena in thin electrolytes, and electrode-electrolyte interfaces.


Introduction
As highly efficient energy conversion devices, capable of converting the chemical energy of fuels, directly into electrical energy, Solid Oxide Fuel Cells (SOFCs) have received much interest from the contemporary scientific community [1][2][3]. One major source of thermodynamic losses generated in an SOFC device is the resistance to the conduction of oxide ions in the ceramic phase of the cell. Thus, it is usually beneficial from a thermodynamic standpoint to decrease the thickness of the electrolyte. Fabrication methods, such as thermal spray, chemical vapor deposition, physical !vapor deposition or pulsed laser deposition allow introducing electrolytic films thinner than 10 µm (as thin as 1 µm) [4,5]. Another interesting technique, known as Plasma-Enhanced Atomic Layer Deposition has allowed fabrication of electrolytes thinner than 0.1 µm. A number of research teams have managed to manufacture working anode-supported SOFC cells with electrolytes as thin as 10 nm, close to the theoretical limit related to dielectric breakdown strength [6]. Baek et al. [6] achieved such low thickness by applying freestanding electrolyte membranes, manufactured using silicon micromachining techniques. Nonetheless, cells with such a thin electrolyte face a number of problems, including durability issues, such as cracks and pinholes, as well as local irregularities stemming from the grain shape [6].
The understanding of the impact of local microstructure inhomogeneities may be improved through the application of a three-dimensional, pore-scale model of charge transport phenomena.
In a non-continuous model of a porous electrode, each node within the computational domain is assigned a specific transport function equivalent to the transport function of the corresponding phase within the porous ceramic-metal electrode. The required digital distribution of phases within the microstructure may be obtained using statistical modeling, or three-dimensional reconstructions of a real electrode, obtained using nanotomography. The latter approach has been made possible by the application of Focused Ion Beam Scanning Electron Microscopy (FIB-SEM) to analyze SOFC electrodes, starting with the ground-breaking works by Wilson et al. [7]. In FIB-SEM nanotomography, layers of sample are removed using a Focused Ion Beam. Subsequently, the uncovered cross-sections are imaged using a scanning electron microscope. A non-continuous electrode model was applied in works by Suzue et al. [8] (for a stochastic reconstruction), Shikazono et al. [9], Kanno et al. [10], Kishimoto et al. [11] Carraro et al. [12] and others.
So far, modeling efforts regarding the impact of electrolyte thickness on electrochemical performance have been relatively sparse, as most works focus on the thermo-mechanical properties of the device [13,14]. A study by Park et al. [15] involved a continuous-electrode, cell-scale model of an SOFC unit cell, including cases with electrolyte thickness ranging from 80 µm to 100 µm. Despite the sensitivity to local microstructure irregularities, as well as the potential to introduce hierarchical microstructures, and mesoscale electrode-electrolyte interface structures [16], the efforts to model electrochemical transport phenomena at microstructure scale have been rare. Iwai et al. [17] used a pore-scale model with non-continuous phase distribution to compare loss generation for several electrolyte thicknesses ranging from 1 µm to 10 µm, demonstrating increase in standard deviation of current densities within electrolytes thinner than 5 µm. Additionally, some influence of cathodic grain size was predicted by the model. Kishimoto et al. [18] furthered this analysis to discuss the impact of fabricating a grooved electrolyte.
The literature survey unravels the knowledge gap and the necessity to investigate the impact of thin electrolyte on microstructure-scale potential distributions in an SOFC. The purpose of the research described in this work is to apply a microstructure-scale, non-continuous model to analyze the impact of cross-electrolyte effects on the performance of an SOFC positive-electrolyte-negative (PEN) assembly more extensively. The cross-electrode interactions are present if there exists a difference between the performance predicted by a model which assumes uniform current flow through the electrolyte, and the model in which no such assumption is made. Previously we have constructed and validated a model of this type, using it for the purpose of analyzing heterogeneity [19,20], and aging phenomena [21]. However, in these studies, the electrodes were investigated individually. In our current paper, we present a model, where the computational domain includes both the anode, and the cathode. The decomposition of thermodynamic losses is presented alongside with distribution of charge transfer and electric potential within a cell. A number of Positive-Electrolyte-Negative assembly simulations with different electrolyte thicknesses are studied, and effects of cross-electrode phenomena are assessed in particular.

Mathematical and Numerical Model
The modeled chemical reaction includes the half reactions of the combustion of hydrogen: The charge transfer model is based on the Poisson conservation with the source term being the local, volumetric charge transfer rate: where φ el (V) and φ ion (V) are the electron-and the ion-conducting phase potentials. σ el and σ ion (Ω −1 m −1 ) are the conductivities in electron-conducting and oxide phases, while i A m −3 is the volumetric charge transfer rate, computed using the Butler-Volmer model: i 0,tpb (A m −1 ) is the specific Triple Phase Boundary (TPB) exchange current density, i 0,dpb (A m −2 ) is the specific Double Phase Boundary (DPB) density, l tpb (m m −3 ) is the local Triple Phase Boundary (TPB) length density, and A dpb (m 2 m −3 ) is the local Double Phase Boundary (DPB) length density. Superscripts 'ano' and 'cath' indicate values specific to the anodic and the cathodic reaction sites respectively. α fwd and α bcw are, respectively, forward and backwards charge transfer coefficient. The specific exchange current densities and charge transfer coefficients for the cathode are taken from studies by Matsuzaki et al. [22], Miyoshi et al. [23], and Kim et al. [24]. For the anode, the data by de Boer [9,25] and Holtappels et al. [26] is used. The model assumes non-continuous distribution of phases within the computational domain, which is based directly on the digital reconstructions from FIB-SEM nanotomography. The TPB and the DPB densities are non-zero only in proximity to the reaction sites. For the sake of simplicity, gas composition gradient in the active layer is neglected. The cell analyzed for the purpose of this study is composed of a Ni-YSZ anode, a YSZ electrolyte, and an LSCF-GDC cathode. The microstructure reconstructions discussed in this paper are taken from our previous study [21,27]. The parameters of each electrode are discussed in Table 1. In our previous research, the transport phenomena on each of the electrodes were solved individually, while the potential drop on the electrolyte was computed separately. Such an approach was sufficient when thicker electrolytes were analyzed. To perform simulations involving thin electrolytes, the model was expanded and tailored to analyze cathodic, and anodic transport phenomena within a single three-dimensional computational domain. Fourier boundary condition is implemented. For each data point, a potential difference is designated arbitrarily. Open Circuit Voltage (OCV) is set as the reference potential value at the electron-conducting boundary of the anode. The scheme of the boundary conditions is visualized in Table 2. The set of equation is discretized using the Finite Volume Method (FVM), with each finite volume in the computational domain corresponding to a voxel within the source FIB-SEM data. The solution is achieved using Successive Over-Relaxation (SOR) method with local linearization of the source term. The methodology is visualized in Figure 1. The resulting model equations are solved using in-house numerical code written in C++.
zb,an z2 zb,cat More details regarding the mathematical, and the numerical model have been provided in our previous papers [19,21].

Results
Validation of the model was carried out by comparison of simulated cell potential to experimental data. The experimental data-microstructure reconstructions and current-voltage relationships-were obtained from a SolidPower S.p.A stack, consisting of 9 anode-supported cells connected in series [27,28]. Each cell consisted of a 240 µm Ni-YSZ anode, 10 µm Yttria-stabilized Zirconia (YSZ) electrolyte and a 50 µm Gadolinium Doped Ceria (GDC) Lanthanum Strontium Cobalt Ferrite (LSCF) cathode. The microstructure data is presented in Table 1. The results of the comparison are presented in Figure 2. For each data point in the chart, the horizontal coordinate corresponds to the voltage measured experimentally, while the vertical coordinate refers to the voltage predicted by the model. A point located on the line represents perfect fit. While some discrepancy between the experiment and the simulation is observed, the agreement is satisfactory, given that reaction and conduction parameters are taken from open literature, rather than being fitted to the experimental data set. The simulation was performed for virtual electrodes with different electrode thicknesses, ranging from 0.10 µm to 10 µm. Since a thin portion of the anode is considered, concentration losses are neglected. The simulated polarization curves are presented in Figure 3. Each line in Figure 3 represents combined thermodynamic losses on positive and negative electrodes for the given current. The distribution of curves demonstrates that as the electrolyte becomes thinner, the losses decrease, although not at linear rate. This indicates the need to explore this phenomenon to a greater extent. Three-dimensional current and potential distributions are presented in Figure 4. It can be seen that most of the voltage loss due to reaction and conduction irreversibility occurs on the side of the cathode, which is relatively thin, considering the thickness of the active layer (indicated by potential gradient). The irregular distribution of ion-conducting phases has a greater impact on the cathode-side, due to high tortuosity of the GDC phase, and lower conductivity of the LSCF phase for this particular electrode. On the other hand, the YSZ phase has high volume ratio and relatively low tortuosity, which-combined with high TPB density-results in lower electrochemical losses.   The generated overpotential can be divided into the components related losses occurring due to the resistance to ionic current conduction (ohmic overpotential), and the reaction irreversibilities (activation overpotential) for each electrode. The results of overpotential decomposition, together with the active layer potential and current distributions are presented in Figure 5. The subfigures presented on the left-hand side show the charge transfer rate between ion and electron-conducting phases (blue line), together with the electric potential of the ion-conducting phase (red dash-dot line), and the electron-conducting phase (red dotted line) for the total overpotential η = 0.125 V. It can be seen that the charge transfer rate distributions are jagged and uneven-this is related to the non-continuous distribution of reaction sites within the computational domain, which is based on the FIB-SEM sample. In the electrodes, the potentials distributions are similar to hyperbolic tangent functions, while in the electrolyte, they appear to be linear. This holds true for all the discussed distributions. As the electrolyte thickness is decreased, the current generated on the electrodes increases. Subfigures presented on the right-hand side show decomposition of total overpotential for different current densities. The maximum current density marked on the charts refers to the highest current density obtained within the presented simulations. For the analyzed current range, the losses appear to be divided more or less evenly between the reaction irreversibilities (activation overpotential) and conduction irreversibilities (ohmic overpotential). Cathodic losses are roughly twice as high as the losses on the anode. Most of the change in polarization appears to be the result of changing electrolyte overpotential. While, as expected, the share of losses generated on the electrolyte shrinks as its thickness decreases, the composition of losses generated on the electrodes does not change significantly. It is clear that as the electrolyte becomes thicker, less current is generated, and more of the overpotential is devoted to overcoming the resistivity of the electrolyte. This is the reason why a comparative analysis of the problem requires normalization, for example by means of isolating and subtracting the losses generated on the electrolyte. Figure 6 displays the potential losses generated on the electrolyte for cells of different thicknesses. As the cell becomes thinner, the losses do not decrease linearly. For every case, halving the electrolyte thickness results in a smaller decrease of overpotential, than what would be expected from a linear relationship. The three-dimensional model is able to predict this behavior since it captures the horizontal component of the current flow through the electrolyte, occurring due to the local irregularities in the microstructure at the electrode-electrolyte interface. Thus, the pathways for the ionic current are longer than what would be expected from a one-dimensional model assuming that the ionic current in the electrolyte flows in a straight line. As suggested before, the overpotential-current relationships derived from the three-dimensional model were corrected to remove the influence of electrolyte by subtracting the averaged potential drop on the electrolyte itself. The corrected polarization curves are presented in Figure 7. Interestingly, the following tendency emerges: the thinner the electrolyte, the smaller the losses on the electrodes. The differences among the polarization curves excluding the losses on the electrolyte suggest the presence of cross-electrode phenomena, as the electrode-only losses decrease to a greater extent than what would be expected from just the removal of the current-resisting material. Most of the difference is present on the cathode side, while the anodic losses remain mostly unaffected by the decrease of electrolyte thickness, increasing slightly as the distance to the cathode shrinks.  Figure 8 depicts the share of overpotential components for different electrolyte thicknesses. Intuitively, the total share of losses generated on the electrodes is smaller for thicker electrolytes. However, the model predicts that the share of losses on the cathode, in relation to electrode-only losses, increases. Thus, it appears that the majority of electrode performance improvement in the presence of a thin electrolyte occurs on the cell's cathode. These tendencies are observed for both the ohmic overpotential, and the activation overpotential. The results from Figure 7, and Figure 8 indicate that the electrolyte thickness indirectly impacts the overpotential by value higher than its ohmic resistance.

Conclusions
In this paper, a microstructure-scale transport model was used to compute potential fields and current distributions in a positive-electrolyte-negative assembly of a Solid Oxide Fuel Cell for several cases, in which different electrolytes of different thicknesses were considered. The model has shown that reducing electrolyte thickness below 10 µm yields diminishing returns, as its relationship to overpotential is non-linear. This is likely an effect of a horizontal component to current pathways, resulting from microstructure irregularities. Although, the results show little effect of electrolyte thickness on the composition of electrode overpotential, some differences were predicted nonetheless.
The current-overpotential relationships indicated that the electrode-only losses would decrease as the electrolyte became thinner. As the current density increases, so do the differences, although they do not seem to exceed 5%. This behavior could be explained by the local phase distribution affecting the pathways of ions penetrating from the cathode to the anode. Additionally, it was predicted that the bulk of the increase happens on the cathode, as the anodic losses did decrease by a small margin. These tendencies were observed for both the activation overpotential, and the ionic conduction overpotential. The relatively low magnitude of the electrode-only losses may be related to the high degree of homogeneity in the microstructures used for the purpose of this study. If the electrodes were more anisotropic, for example featuring a periodic inhomogeneity, the cross-electrode effects would likely be more prominent. The results suggest that transport on a thin electrolyte is not one-dimensional, pointing towards the importance of using a three-dimensional microstructure scale model to tackle this research problem. When electrodes are analyzed separately, which is common in the literature, the obtained results may differ compared to the simulation conducted on the entire PEN structure, particularly when the electrolyte is thinner than 10 µm. Acknowledgments: This research was supported in part by PL-Grid Infrastructure. Matplotlib Python library was used for data visualization [29].

Conflicts of Interest:
The authors declare no conflict of interest.