Microstructure Evolution of Solid Oxide Fuel Cell Anodes Characterized by Persistent Homology

A Solid Oxide Fuel Cell consists of solid porous electrodes separated by ceramic electrolyte. The electrodes are composed of a ceramic metal composite, in which every material has a specific mass or charge transport function, such as electron conduction, ion conduction, gas diffusion. The pathways for mass and charge flow are not straightforward, as they are dependent on the microscopic arrangement of conducting materials and open pores. Furthermore, the electrode-specific half reaction may only occur if reaction sites, such as triple phase boundaries or double phase boundaries are available. The overpotential generated in the microstructure accounts for a significant part of the voltage losses generated in an SOFC device. Currently microstructure-scale modeling relies on direct association of the electrochemical performance to physical attributes of a microstructure, such as phase volume fraction and phase tortuosity, grain size, availability of reaction sites. Microstructure-oriented analysis is generally based on a material reconstruction - a digital or analytical model of a micro-scale representative volume within the studied material. The digital material reconstruction may be based on either experimental measurements obtained using microscopic techniques, or an artificially generated model i.e. a collection of packed primitives, such as spheres. In the case of SOFCs the latter include stochastic microstructure reconstruction methods such as sphere-packing simulation Schneider et al. (2006); Nishida and Itoh (2011), in which the electrode manufacturing process (powder mixing and sintering) is approximated. When a reconstruction is available, the properties of a microstructure may be approximated using experimental measurements complimented by analytical methods. stereological methods which extrapolate the three-dimensional distribution of phases in the microstructure from a limited number of cross-sections obtained by means of scanning-electron microscopy (SEM) Yeong and Torquato (1998). Recently, advancements in tomographic techniques, such as X-Ray tomography, Focused Ion Beam Scanning Electron Microscopy, and their derivatives (i.e Plasma Focused Ion Beam SEM, Broad Ion Beam SEM) have allowed for preparation of three-dimensional digital reconstructions of SOFC microstructure Wilson et al. (2006) corresponding directly to the distribution of material in a sample. The microstructure reconstruction may then be used to prepare a computational domain for solving mass and charge transport equations Suzue et al. (2008); Kanno et al. (2011). Alternately, the principles of the continuous electrode theory Newman and Tiedemann (1975) may be employed, and the electrode may be modeled as continuous medium, in which all the transport processes occur in parallel. The microstructure properties are then accounted for using macroscopic parameters, such as phase volume fraction, phase tortuosity, reaction site density Iwai et al. (2010). The microstructure reconstruction may be used to estimate the aforementioned parameters: the volume expansion method Iwai et al. (2010), the midpoint method may be used to estimate triple phase boundary density Kanno et al. (2011), the marching cubes method may be used to estimate interphase surface area Lorensen and Cline (1987), diffusion simulation methods, such as the random walk method may be used to estimate sample tortuosity Iwai et al. (2010). Alternately, empirical correlations may be implemented to estimate the transport parameters using grain size and volume fraction. However, the assumptions at the foundations of these methods may render them inaccurate in many applications Tjaden et al. (2018).

The disadvantages of the quantitative approach related to the complicated process of obtaining data, together with the significantly limited possibilities of observing direct changes in the structure of the topological microstructure, present the necessity to propose an alternative method of describing the electrode material. Thanks to the use of tomography techniques to represent the SOFC anode material, it is possible to apply topological data analysis (TDA) methods, in particular persistent homology (PH). This solution allows for the observation of the characteristic properties of a given scale of the analyzed object and presenting them in the form of a persistence diagram (PD) or after appropriate data processing as a persistence image (PI). Persistent homology is a relatively young and still under development topological tool Edelsbrunner et al. (2000), but it has already found application in many fields of science and engineering. TDA methods are particularly popular in medical research. The SARS-Cov-2 virus analysis by Chung and Ombao (2021) used persistent homology to represent the topological properties of the closed and open form of the coronavirus and to compare the differences between human and feline virus. The results of the research presented in the form of a persistence diagram allowed for the identification of characteristic properties for individual variants of SARS-Cov-2, which, according to the authors, may allow for the development of a faster and more accurate automated method for recognizing a given type of virus in the future. No less important analyzes based on TDA methods are also performed in other branches of medicine, i.e. neurobiology. In the publication of Bendich et al. (2016) researchers used persistent homology to describe the aging changes of brain tree arteries. On the basis of the obtained results, it was observed that even a significant filtration of the PH data allows for a satisfactory description of the aging factor impact. Research using persistent homology has found application in hepatocellular ballooning in liver biopsies analyzes (Teramoto et al. (2020)), where persistence images were used, and in oncology (Bukkuri et al. (2021)), allowing better observation of the cancer architecture. Another field of science that uses the development of TDA methods is geology. In the work of Robins et al. (2016) various two-phases systems were studied using persistent homology. These analyzes showed the correlation of the percolation threshold not only for 0th homology (connected components) and an improvement in determining the percolation radii was found if the data was previously filtered. The team of Suzuki et al. (2020) confronted the results obtained with topological data analysis with conventional methods, where in some cases the use of persistent homology provided additional information that was impossible to observe with the classical approach. TDA methods have also found application in flow analysis. In the research of the Suzuki et al. (2021) authors used persistent homology, in particular 1st order homology (loops or open holes) to model physical phenomena occurring in porous material. Flow studies based on topological methods showed a high compliance with classical simulations and, importantly, there are no contraindications to apply this approach in the case of mass and energy transport processes. The possibilities of topological data analysis in the field of fluid mechanics were presented by Kramár et al. (2016) using persistence diagrams in both the vorticity during Kolmogorov flow and the Rayleigh-Bernard Convection temperature field studies. These analyzes confirmed the stability of persistent homology in the case of limiting the simulated values, which was observed by preserving the characteristic properties for the Rayleigh-Bernard Convection. TDA methods have found application in many other fields of science (Fedotov et al. (2022); Hirata et al. (2020); Ichinomiya et al. (2020)), offering the possibility of insight into the structure of various systems. Based on the available scientific sources, it is clear that persistent homology seems to be the right tool for the analysis of the porous material of the solid oxide fuel cell anode.

This research aims to use computational topology in the new field, which is solid oxide fuel cell microstructure analysis. Apart from the new application of computational topology, we also use a unique set of data that includes a three-dimensional digital representation of anode microstructure before and after long-term operation. The research points to a new direction in fuel cell analysis and lays down the background to improve the utilization of deep learning for fuel cell study.

Phase volume fraction, phase tortuosity factor, and reaction site density (either the triple phase boundary, or the double phase boundary) are often a necessary and sufficient input for a mathematical model of voltage losses in a small SOFC device. This has been demonstrated in one of our studies, where a model using microstructural data was validated against the results of a 300 W SOFC stack experiment Brus et al. (2017). Multiscale model, such as the model developed and validated by Mozdzierz et al. Mozdzierz et al. (2019) can further account for phenomena occurring apart from the microstructure, such as the irreversibilities of flow in an electrode channel.

Previously, we have performed a number of solid oxide fuel cell long term performance studies, in which the emphasis was put on the microstructure evolution Brus et al. (2015a, b). A 100 W 6-cell SOFC stack (an in-series or parallel connection of multiple cells), manufactured by SolidPower SPa of Italy, was operated over a period of 3700 h. The voltage generated by the system was being measured continuously, revealing not only a lack of deterioration of stack electrochemical parameters, but also a slight voltage increase of about 2.58%. After the long-term operation experiment was completed, the stack was disassembled, and the cells were investigated using Focused Ion Beam Scanning Electron Microscopy (FIB-SEM). Similar measurements were performed on a reference cell (reference cell being a brand-new cell of the same composition, but one which did not participate in the long-term performance experiment). FIB-SEM study allowed us to obtain three-dimensional digital reconstructions the microstructures within the cells. The reconstructions were then analyzed using tomography data processing software (Avizo by Thermofischer Scientific/FEI) as well as in-house code to quantify the microstructural parameters such as Triple Phase Boundary (TPB) density and tortuosity of phases. The TPB density was determined to be a relevant parameter for the analysis, since it is the main reaction site within an SOFC anode, and its deterioration has often been perceived as one of the main causes of the SOFC degradation. On the other hand, tortuosity and volume fraction increase the resistivity to current conduction and gas diffusion. The studies have revealed Brus et al. (2015a, b) that for all of the analyzed samples the TPB density decreased in relation to the values measured in the reference cell. But so did the tortuosity of the ion-conducting phase and the open pores. The degree to which the cells’ microstructural parameters have changed were found to vary depending on the position within the stack. Additionally, the parameters of cells investigated after the long term performance experiment were found to exhibit higher anisotropy - the tortuosity measured in different directions varied more significantly for samples recovered after the stack experiment Brus et al. (2020). Another long term stack experiment involving the SolidPower stack of the same type has shown a similar tendency. In another study, the microstructure reconstruction was used to construct a three-dimensional computational domain for simulating microscale mass and charge transport in the electrodes Prokop et al. (2019). Empirical relationships for material and substance parameters were used, along the common Bulter-Volmer model in order to solve the Poisson differential equation set for conservation of mass and charge. A numerical solution was obtained for microstructure samples taken from the reference cell, and the cells investigated after the degradation experiment. The numerical model predicted a decrease in voltage loss for the post-degradation samples based on their microstructures, thereby supporting the hypothesis that the increase of terminal voltage was the result of microstructure evolution. The digital microstructure reconstructions were later used in a pilot study to determine the main cause of the microstructure evolution. By quantifying the interfacial free energy (a thermodynamic state function, which is expected to monotonously decrease for spontaneous processes) it was shown that the performance enhancement was likely the result of an expected nickel coarsening or a similar internal processes. Prokop et al. (2021).

In order to present the process of the creation of the persistence diagram, it is necessary to introduce a number of concepts related to the TDA methods. Homology is the topological invariant (a parameter that is immutable during the transformation of a given topological space) that represents a set of $n$-dimensional holes (Obayashi et al. (2018)). It characterizes the topological properties of a given scale of the analyzed object, allowing for comparison of two different topological spaces (Pereira and de Mello (2015)). The application of persistent homology requires the discretization of continuous space, which is performed by using the so-called complexes. There are many different types of methods used to represent a topological space. For the purpose of the research contained in this paper, the simplicial complexes are considered. This method consist in building cells based on simplices representing the simplest geometric objects with flat faces. A simplical complex $S$ should be understood as a set of simplices consisting of: points (0 - simplex), segments (1 - simplex), triangles (2 - simplex), tetrahedrons (3 - simplex) and other multidimensional simplices presented in Fig.3 (Zhu (2013)) therefore, it can be concluded that for a given topological space $\mathbb{M}$, homology assigns a set of complexes $H_{k}(\mathbb{M})$ describing the topological properties of $\mathbb{M}$, where $k=0,1,\dots$ is the order of the homology. The dimension $H_{k}(\mathbb{M})$ is known as $k^{\mathrm{th}}$ Betti number and is commonly denoted as $\beta_{k}$ (Edelsbrunner (2012)). The purpose of this parameter, which is a topological invariant, is to count the $k$-dimensional features of the studied space. Equivalent topological spaces have the same value of the Betti numbers. In this paper, the three-dimensional space will be examined, hence the $k$ = 0, 1 and 2 are considered, denoting connected components, loops (cycles) and voids. The main goal of persistent homology is therefore to calculate homology on many scales while observing topological properties of the system that persist across those scale. The term ’persistent’ is associated with the appearance and vanishing of the observed properties of a given order of homology as a result of carrying out subsequent stages of topological space filtration. When considering a certain set of topological data (see Fig.4 (a)) in order to extract the characteristic properties of a system, complexes are constructed. For the purpose of filtering presented data set, the $r$-ball model was applied such that:

where $T_{r}$ is the set of all points distant from any element of the $P$ by $r$, $B_{r}(x_{i})=\{y\in\mathbb{R}^{N}:||y-x_{i}||\leq r\}$ is a ball centered at $x_{i}$ with radius $r$ and $N$ denotes the dimensionality of space Obayashi et al. (2018). As the radius increases segments (1 - simplex) are created which delimit a certain two-dimensional space (Fig.4 (c)). Thanks to this process $\mathrm{PH}_{1}$ homology (cycle) is distinguished. Further radius enlargement results in a creation of a face (2 - simplex), which causes the loop to disappear completely (Fig.4 (d)). In order to determine the persistence of the topological properties, the whole process of space filtration is monitored where the lifespan of individual homology is described by a segment from birth to death (understood as appearance and disappearance of the feature).

Constructing a persistence diagram from a barcode is relatively easy - instead of drawing a segment of the length from appearance of a topological property to its disappearance, a point with coordinates of birth and death is plotted on the graph. Therefore, persistence diagram $\mathrm{PD}_{k}$ should be understood as representation of the set of all $k^{\mathrm{th}}$ homology elements. Usually in the PD diagram x-axis corresponds to the birth of a given topological property, while the y-axis represents its death. A characteristic element of PD is the diagonal $f(x)=x$ dividing the plot area into to regions, where, due to the nature of data representation, the part located under the line is always empty (it is impossible for a topological property to die before its appearance). Points lying on the diagonal and in its direct proximity are categorized as noise (due to their negligibly short life), while those positioned at a considerable distance from function $f(x)=x$ represents proper topological properties.

The representation of topological data in form of PD, although relatively convenient, does not always work well in further computer analyzes (Adams et al. (2017)). Therefore, there was a need to simplify the presentation of data used in subsequent stages of research based on methods , i.e. machine learning or decision tree algorithms. Such a solution appears in the so-called persistence image (PI), which is the conversion of the input PD to a finite-dimensional vector. This method not only retains the connection to the original persistence diagram and the topological properties it represents, but is also resistant to noise, efficient to compute, and has the ability to adjust the representation for a given system. The first step in converting the topological data representation of the system into PI is to perform linear transformation $N:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ of the original PD, such that: $N(x,y)=(x,y-x)$ Adams et al. (2017). Then so-called persistence surface (PS) is created by applying a weighted sum of functions defined with a distribution $g_{u}$ centered at each point in the PD. In the research presented in this paper normalized symmetric Gaussian was used. For $g_{u}:\mathbb{R}^{2}\rightarrow\mathbb{R}$ with mean $m=(m_{x},m_{y})$ and variance $\sigma^{2}$, the $g_{u}$ is defined as Adams et al. (2017):

and persistence surface $S_{N}$ as:

The introduction of the weight function $w$ is due to the necessity to ensure the stability of the persistence surface. The last transformation leading to PI is done by integrating $S_{N}$ over sub-areas $p_{i}={(x,y):p_{j}\cap p_{k}=\emptyset,\bigcup_{k}p_{k}=\mathbb{R}^{2}}$ of a discretized domain, which can be represented as Adams et al. (2017):

The representation of topological data in form of persistence image is influenced by a number of parameters selected depending on the system under consideration. The important issue is the selection of the appropriate grid (often recognized in the literature as resolution) for the persistence image. Research shows (Adams et al. (2017)) that the data presented by means of PI relatively well preserve the representation of characteristic topological properties even for significantly reduced resolution. Another crucial aspect worth consideration is defining the appropriate weight function for a given problem. It was documented that the points with the longest persistence are not always the most essential in the analysis of a studied system (see Bendich et al. (2016)). This parameter should be selected based on the type of problem under consideration. Applying a different distribution is also possible and determining the correct value of the variance should be examined. Transformation from persistent diagram into persistent image are presented in Fig. 5. As persistnce landscape and image require additional parameters, that have to be selected for given application, in this work only persistence diagrams will be presented and compared.

When applying the TDA methods to represent microstructural changes in anode material, a valid question arises about the accuracy of the technique and its sensitivity to changes during the aging process. While geological and, above all, material science studies using persistence diagrams are particularly valuable from the perspective of the SOFC fuel electrode, they do not fully answer the possible problems that should be considered when analyzing this material. The possible lack of projection of the loss of continuity of a given phase or the appearance of a new tunnel in the case of SOFC cells may be associated with a misrepresentation, resulting in erroneous predictions as to the operating parameters of the entire device. For these reasons, it was decided to conduct synthetic microstructure analyzes consisting of two phases in form of spheres and a third phase filling remaining space, as shown in Fig 7. Three series of tests were performed, analyzing the scale and significance of changes in the realization of the same random fields for system with the same parameters (volume fraction, particle size, sample size), with increasing the diameter of the spheres in subsequent samples and with the varying particle sizes inside the sample drawn from normal distribution. Selected $\mathrm{PD}_{0}$ and $\mathrm{PI}_{0}$ for the test series are presented in Fig. 8. Research shows that the representation of synthetic microstructures is relatively stable to the influence of various factors. The differences for the first series are, as expected, minimal - there are minor disproportions between the diagrams (especially $\mathrm{PD}_{0}$), but the character typical for the order of the PD is retained in each case. It should be emphasized that the discrepancies are clearly understandable - they stem from the differences in the distribution of the spheres within the sample. The results obtained for second series, although more diversified, allow to distinguish even in a visual approach the direction of changes taking place in the microstructure while maintaining the general characteristics of the representation. It is evident that with the increasing radius of the spheres more persistent features appear. The greatest contrast between the representation is observable for the series that introduced varying particle sizes within the sample. However, even in this case, it is possible to distinguish the direction of increasing "dispersion" of birth-death points in relation to the reference case. The data representation in the form of PI also shows similarities and characteristic trends of topological changes even for third series. Nonetheless, it should be emphasized that selected linear weight function marginalizes distant single points of high persistence in favor of a large number of features in close proximity of diagonal, which results in less visible differences between images. The conducted tests for synthetic microstructures should not be extrapolated without reflection for real SOFC anode materials, but nevertheless they demonstrate the PD representation abilities and indicate possible directions of analyzes of the obtained data.

To test differences between new and after long-term operation, we have chosen cell no. 5. The selection of the cell is dictated by the possibility of comparing observations based on TDA with data obtained from the quantitative approach methods presented in our previous works Brus et al. (2015a, b). $0^{\mathrm{th}}$ order diagrams for Nickel phase are juxtaposed in Fig. 11. In this case, a greater number of features being above 0 in the $Y$ axis are observed in comparison to the representation of the reference sample. When analyzing the microstructure in relation to the sampling site, there is a noticeable trend in the decreasing number of birth-death points with positive vertical coordinate along the fuel flow direction, with an improvement in persistence of the features for diagrams of each sample.

$\mathrm{PD}_{0}$ diagrams for YSZ are presented in Fig. 12. The representation of the connected components of YSZ phase of the cell subjected to the degradation process is similar to the corresponding diagram of the reference sample - relatively small differences between individual PDs are observed, but they do not have the specific tendencies.

The most significant differences in $\mathrm{PD}_{0}$ of the pore phase, which can be seen on Fig. 13, are visible similarly to nickel in the section of the diagram located above zero on the $Y$ axis, with the downstream microstructure resembling reference sample the most. There was no characteristic trend of changes in features in the diagrams for the pore phase depending on the place of sampling.

It should be emphasized that all $\mathrm{PD}_{0}$ representations retained the general character of the points distribution typical for the respected phase. Comparing the above observations to the previous analyzes based on the quantitative methods, the tendencies of topological changes for the nickel phase was consistent, while the study of the continuity parameter for samples taken upstream and center sites did not show such significant differences as those presented on the diagram. The share of the YSZ percolated cluster determined on the basis of the quantitative methods was 100% for both reference sample and all microstructures of the cell at the top of the stack. These data show partial consistence with the observations of the topological properties representations - all samples were similar, while some birth-death points appearance negate the complete phase continuity for each of the microstructures. This difference is caused by the assumption of the continuity of the phase located at the sample boundary. The lowest consistency is found in the data for pores, which in the case of quantitative approach methods indicate an improvement in phase continuity, although indeed the downstream microstructure is the most similar to reference sample. It should be emphasized, however, that the visual assessment of the data contained in $\mathrm{PD}_{0}$ of the pore phase was inaccurate due to the inconsistency of changes between the individual diagrams.

PD diagrams of $1^{\mathrm{st}}$ order of nickel phase are presented in Fig. 14. For the diagrams corresponding to the nickel phase, some changes in the trend of the distribution of the birth-death points depending on the place of material collection were observed, while the nature and direction of these changes is difficult to determine.

In the case of YSZ, as seen in Fig. 15, the representations of topological properties are extremely close to each other, and the scale of changes taking place as a result of the degradation process is unnoticeable.

Persistent diagrams of 1st order for pore phase are presented in Fig. 16. The greatest changes in $\mathrm{PD}_{1}$ are observed for the pore phase, where in relation to the reference sample there is a greater number of tunnels appearing in further stages of filtration and this tendency correlates with the direction of fuel flow.

$\mathrm{PD}_{2}$ of the YSZ phase is characterized by a much greater number of the presented features than for nickel, as shown in Fig. 18. Comparing the diagrams of the microstructure subjected to the aging process with the material from the brand new cell, a general increase in the number of birth-death points and appearance of a relatively large group of individual properties with very high persistence, which were born in the final stages of filtration, are observed.

For the $\mathrm{PD}_{2}$ representation of the pore phase, shown in Fig. 19, a certain trend in the distribution of topological properties can be noticed depending on the place of sampling. Starting from the reference sample and then carrying out the analysis of microstructures in the direction of the the fuel flow, the "flattening" along the diagonal of the main grouping of birth-death points is observed, with a simultaneous increase in single points of high persistence for each subsequent sample.

The article used computational topology to analyze microstructural data from solid oxide fuel cells before and after long-term operation. It is the first report that uses persistent homology for microstructure analysis of solid oxide fuel cell electrodes. Summarizing the presented observations, it can be concluded that even a strongly imperfect visual analysis of topological data allows for the initial selection of information about microstructure and distinguishing trends of topological changes taking place in the material. In the case of $\mathrm{PD}_{1}$ and $\mathrm{PD}_{2}$ analyzes, there is no benchmark, which for $\mathrm{PD}_{0}$ were the data obtained using the quantitative methods. However, it should be emphasized that ultimately the topological data analysis methods should cooperate with deep-learning algorithms, which would allow the identification of characteristic changes impossible to observe by the human eye.

Research project supported by the program "Excellence Initiative – Research University" for the AGH University of Science and Technology.