Microstructure Evolution in a Solid Oxide Fuel Cell Stack Quantiﬁed with Interfacial Free Energy

: Degradation of electrode microstructure is one of the key factors affecting long term performance of Solid Oxide Fuel Cell systems. Evolution of a multiphase system can be described quantitatively by the change in its interfacial energy. In this paper, we discuss free energy of a microstructure to showcase the anisotropy of its evolution during a long-term performance experiment involving an SOFC stack. Ginzburg Landau type functional is used to compute the free energy, using diffuse phase distributions based on Focused Ion Beam Scanning Electron Microscopy images of samples taken from nine different sites within the stack. It is shown that the rate of microstructure evolution differs depending on the position within the stack, similar to phase anisotropy. However, the computed spatial relation does not correlate with the observed distribution of temperature.


Introduction
Solid Oxide Fuel Cells (SOFCs) are energy conversion devices, which provide electrical power by converting the chemical energy of fuels such as hydrogen, carbon monoxide, or carbohydrates (if internal or external reforming is used). Despite their numerous advantages-high efficiency, wide variety of available fuels, and low pollution-the wider market application of SOFCs is still hindered by several issues, not the least of which the microstructure degradation affects the reliability of the ceramic porous electrodes.
A ceramic electrode is composed of several phases, and each phase has a specific transport function. An oxide phase, e.g., YSZ (Yttrium-Stabilized Zirconia) conducts ionic current, metallic phase (e.g., Nickel), and gas reagents diffuse through the open pores. The reaction of oxidation occurs at the Triple Phase Boundary (TPB). As a result, even relatively small changes to the distribution of material in the electrode cause significant changes to the electrochemical performance.
Understanding the SOFC microstructure evolution is aided by empirical studies involving long-term stack performance experiments, combined with the analysis of microstructure changes which occurred over during the experiments. Monaco et al. [1] conducted X-ray holotomography-based microstructure analyses for Solid Oxide Cells operating for periods ranging from 1000 h to 15,000 h. Brus et al. performed two independent long performance tests in which an enhancement, rather than deterioration of performance was observed [2][3][4][5].
One way to quantify the evolution of multiphase systems is by analyzing the change of its free energy. This is the basis of Phase Field Models. In PFMs, the phase distribution is diffuse, rather than discrete-they are modeled using a set of continuous functions. The interface of any two phases has a finite, non-zero thickness within which the phase volume fractions take noninteger values. Free Energy of a microstructure may be computed using a Ginzburg-Landau type free energy functional [6]. The rate of phase transition may be determined using Allen-Cahn [7] equation for nonconserved order parameters, and Cahn-Hilliard [8] equation for conserved order parameters. Depending on the scale and the physical input data, Phase Field Modeling may be considered to be either a numerical tool or a first principles' approach [9,10]. PFMs have been used in a broad range of problems, ranging from microstructure evolution and solidification [11] to multiphase flows [12].
Since porous electrode degradation due to coarsening of microstructure can be modeled as an effect of a monotonous decay of microstructural free energy, PFMs have also been applied to SOFC related problems by several research teams. Chen et al. [13] compared TPB density and tortuosity changes in an SOFC anode, as predicted by two variants of a PFM model. Cahn-Hilliard equation is generally used for the distribution of phases, and Allen-Cahn equation is employed for grain-related order parameters, as was the case in studies by Lei et al. [14] or Jiao and Shikazono [15] among others.
Jiao et al. [16] combined a PFM model with a microscale transport phenomena model to interpret data from a long-term performance study. Shimura et al. [17] performed a similar analysis for the cathode. Recently, Wang et al. [18] used distribution of phases derived from PFM of anode evolution as input into a 3D microscale model to compute the electrode polarization.
In this paper, we use Ginzburg-Landau type free energy functional to quantify the degradation observed during a previous long-term performance stack experiment. One of the purposes of our analysis is to investigate the factors affecting the degradation rates, as well as its relationship with microstructure anisotropy. While it is a common practice to solve the transient PFM employing Allen-Cahn and Cahn-Hilliard equations, we believe that the necessary conclusions can be drawn from a simpler model, since the free energy is a monotonously declining state function. While a transient model allows prognosticating the quality of postdegradation microstructures, in our case they have already been determined by experiment. Relating these results to the framework of PFM modeling will allow for a quantitative description of the observed microstructure evolution.

Experimental Analysis
The analysis presented in this paper is concerned with the results of a long-term performance experiment during which a stack of 6 cells was operated within a modular stack test bench (MSTB). Both the SOFC stack and the MSTB were produced by SOLIDPower S.p.A. of Italy. The stack consisted of six 6 cm × 8 cm anode-supported cells connected in series. The total of 10 microstructural samples were recovered from the stack, including 1 sample from a reference cell, corresponding to brand new microstructure, and 9 samples from the stack which was disassembled after a 3800 h degradation experiment. Nine samples were taken from the inlet-sides, the centers, and the outlet-sides of three cells: the top cell, the central cell, and the bottom cell. The schematic visualization of the site naming scheme is illustrated in Figure 1b. The experiment has shown that the temperature distribution within the stack, while mostly constant over the course of the study, was not uniform. The temperature readings are displayed on the experimental setup diagram in Figure 1a. The samples were imaged using FIB-SEM setup, and the resulting stacks of images were digitized into three-dimensional data sets using AVIZO (Thermo Fisher Scientific). The methodology of FIB-SEM, is presented in Figure 1c, and the sample crosssections obtained before and after the long-term performance test are presented in Figure 2. The long-term performance experiment has been previously described in another paper [2]. Three-dimensional visualizations of the samples are presented in Figure 3. Interfacial areas of phases in the samples were computed using Marching Cubes Method [19], as implemented in AVIZO. Additionally, in the paper by Brus et al. [4] we have estimated the anisotropy of the samples. The anisotropy, understood as the measure of differences in parameters along different directions, was quantified using the sum of anisotropy factors ζ i for each phase (Nickel, YSZ, Pore) such that: where τ is the tortuosity factor, defined as the ratio of the diffusion coefficient in a free space to the diffusion coefficient in a real sample in a porous electrode. τ i is the sample average tortuosity factor of phase i in the direction i (x, y, z being the standard directions), and τ ir is the averaged tortuosity factor of phase i for any direction. The tortuosity factors were estimated using Random Walk Simulation. ζ = 0 when the tortuosity factor of the sample is perfectly uniform. The anisotropy calculation methodology was previously presented in the cited open access paper by Brus et al. [4].

Mathematical Model
The Ginzburg-Landau type free energy functional for a two-component system can be described using the following equation [6].
The value ψ i , and ψ j are local phase volume fractions corresponding to a given set of two phases i, j (see Figure 4), such that ψ i + ψ j = 1. Since the phase distributions are continuous, a single point in the computational domain may have noninteger values.
is the bulk free energy of the system. The bulk free energy of an interface shared by phase i and another phase j can be computed using a double-well potential equation, which is commonly ( [6]): for a given pair of two phases i, j such that is the energy barrier between minima in the double well potential equation for a diffuse phase interface between phase i and the other phase j. When considering any two-phase interface ij (ij corresponding to Ni-YSZ, Ni-Pore, and YSZ-Pore pairs), its specific interfacial energy per unit area γ ij (J m −2 ) may be obtained by finding the interfacial phase volume fraction profile at a Double Phase Boundary , for which the interfacial energy is locally minimal. The analytical solution for the ψ(x) profile may take the following form [20]: where x (m) is the distance in the direction normal to the interface, and δ ij (m), the ij phase interface thickness equals to: The analytical solution may be used if the diffuse interface between a pair of two phases ij is one-dimensional. This assumption appears to be acceptable if a very small section of the sample is selected, and a thin interface is considered. For the ψ profile in Equation (4) γ ij (J m −2 ) may then be computed using the integral (2), yielding the following relationship [6]: The two-phase diffuse interface is visualized in Figure 4 In all of the analyzed cases, the interface gradient was assumed to be constant and perpendicular to the discrete interfaces in the digitization of the microstructure, for interfaces located far from a triple phase boundaries (TPBs). The phenomenological parameters ∆ f and 2 in Equations (2), (3), (5) and (6) can be estimated empirically. Using Transmission Electron Microscopy (TEM) Nahor et al. [21] measured Ni-YSZ interfacial energy at γ NiYSZ = 1.8-2.1 J m −2 depending on configuration. SEM measurements of SOFC atomistic structure suggest interface thickness of 0.242 nm [22]. Similarly, Chen et al. [13], use interfacial energy ratio of γ NiYSZ :γ NiPore :γ YSZPore = 2.2:1.9:1.4 and these are the values which we employ in our study. The values are used as input to Equations (5) and (6), a set of which is solved for the phenomenological parameters and ∆ f , yielding: From the relationships (7), one can obtain 2 NiYSZ = 4.86 × 10 −10 J m −1 , and ∆ f NiYSZ = 8.09 × 10 9 J m −3 . Using these values to solve Allen-Cahn or Cahn-Hilliard equations for the entire microstructure would be impractical due to high computational cost, since a sub-angstrom scale mesh would be required to accurately represent the diffuse interface. In this work, the problem is bypassed by, first-considering only two specific points in time (pre-degradation, and post-degradation cell samples), and secondidentifying several repeating spatial configurations in the computational domain, for which the local solutions are computed independently. The local area in a microstructure reconstruction discrete volume which shares an interface with a different phase is multiplied by the interfacial free energy density corresponding to the given pair of phases (γ NiYSZ = 2.2 J m −2 , γ NiPore = 1.9 J m −2 γ YSZPore = 1.4 J m −2 ). While local profiles for the interfaces located far from a triple phase junction may be viewed as essentially onedimensional for a first-principles based interface thickness, more consideration needs to be given for a near-TPB location. An analytical solution would be difficult or impossible to obtain in proximity to a junction [6]. In their "Study on three-component Cahn-Hilliard flow model", Boyer and Lapuerta [23] propose the following generalization of Equation (2) for three-phase systems: which we use in our current work. In this case, ψ i , ψ j , and ψ k such that ψ i + ψ j + ψ k = 1 are local phase volume fractions corresponding to a given set of three phases i, j, k (Ni, Pore, YSZ in the case of the discussed SOFC anode). Boyer and Lapuerta [23] propose the following f 0 (ψ i , ψ j , ψ k ) function for a three-component system.
. Then: Σ i is the spreading coefficient of phase i. In a system containing phases i, j, k, Σ i can be computed from the specific interfacial energies per unit area γ (J m −2 ) of each of the component pairs in the system: The set of (ψ i , ψ j , ψ k ) -which is used in Equation (9) to compute the specific interfacial energy for locations in the vicinity of a TPB-is obtained by minimizing F(ψ i , ψ j , ψ k ). To do so, a transient phase evolution equation F(ψ i , ψ j , ψ k , t) is solved with an arbitrarily high value of time t, at which F approaches a minimum. Since phase concentrations are treated as conserved-order parameters, Cahn-Hilliard equation is used for that purpose. The multiphase implementation of Cahn-Hilliard differential equation set, based on the work of Boyer and Lapuerta [23], can be written as follows: where µ i (J m −3 ) is the chemical potential of phase i, M 0 (m 3 s −1 ) is reference phase mobility set arbitrarily to 0.01. M 0 can be set arbitrarily, since the purpose of the equation is to locally minimize F at a given point in time, rather than to prognosticate further evolution of the microstructure. The dimensionless coefficient θ, is given by Equation (12).
where the indices i, j, k signify Ni, Pore, and YSZ phases, respectively. Since ψ Ni + ψ YSZ + ψ Pore = 1, a set of two partial differential equations is sufficient. A near-triple phase boundary area is chosen to compute a specific microstructural free energy associated with the presence of a phase junction. The Triple Phase Boundary (TPB) junction is modeled as follows: a circle is divided into three equal parts corresponding to Ni, pore and YSZ phases. Periodic boundary conditions are implemented. Then the energy minimization simulation based on the coupled set of Cahn-Hilliard equations (see Equation (11)) is run for an arbitrarily long time, until little to no change is observed. The resulting diffuse interfaces near a triple phase junction (which, in the demonstrated two-dimensional case, is a point on a plane) are presented in Figure 5. The near-TPB free energy integral F TPB (J m −3 ) is computed using Equation (9). The integration limits are restricted to ensure focus on the junction and to keep the lengths of two-phase boundaries within equal to each other. The energy in a near-TPB voxel in a real, discrete reconstruction of a microstructure is computed by adding the required interfacial energies given by Equations (8) and (11). Since three-dimensional samples are considered, it is possible to discuss the spatial distribution of the free energy within the computational domain. To do so, the local values for a given depth z of the sample (depth being the FIB-SEM milling direction) are computed by integrating the Equation (8) for a sub-volume of the computational domain corresponding to a single FIB-SEM slice: where z(m) is the depth coordinate, and ∆z(m) is equal to the voxel edge length. Figure 6 shows the sample-average specific free energy densities (interfacial energies) computed for microstructure samples obtained at different sites within the stack. As expected, it can be seen that the free energy decreased during the long-term performance experiment. The decay of free energy-that is, the difference between the interfacial free energy in the brand new cell (reference) sample, and in the post-degradation samples-can be used as an estimate of the degradation rate. The free energy decrease is largely related to the drop of the specific surface area. As seen in Figure 6, the most significant decrease was observed for the middle cell samples. On the other hand, the post-degradation values of free energy density for the top cell and the bottom cells decreased to a lesser extent. The middle cell, and the top cell appear to have little variation of F among the sampling sites. On the other hand, the downstream side of the bottom cell seems to have experienced a higher degree of decay, compared to its upstream side, and its center.

Results and Discussion
It is interesting to compare these results to the increase of anisotropy, another process linked to the microstructure evolution. The quantitative anisotropy factor was discussed in Section 2. The results of anisotropy studies, which were previously published by Brus et al. [4], showed that the middle cell experienced the greatest increase of the anisotropy factor. Similarly, the current analysis shows that F decreased the most for the middle cell. Furthermore, the downstream side of the bottom cell was shown to become more anisotropic than its center and upstream sampling sites [4]. This is analogous to the higher decrease of the free energy density in the bottom cell's downstream sample reported in the current study.
The relationship between the phase distribution of the postdegradation microstructures and their electrochemical performance was discussed in our previous paper [5] for a similar planar SOFC stack. In that case, the reduction of the cell's overpotential was observed due to the decrease of pore tortuosity-a development predicted by several PFM-based analyses reported on in the literature. However, in our study [5], it was shown that the middle cell experienced less enhancement than the extreme cells. in the samples from the brand new, reference cell. It can be seen that the distribution of local free energy within the sample (computed as the integral of the free energy Functional (9) over a single "slice", see Equation (13)) is more evenly distributed in most of the samples recovered after the long term performance experiment. The depth of the sample z(m) refers to the coordinate in the direction of FIB-SEM milling. In general, the local values in post-degradation samples do not exceed the local values in pre-degradation samples for a given z-axis position within the sample. While the sample sizes vary, the size of each representative volume used in the study is ∼1000 µm 3 which is considered a sufficient size for the purpose of microstructural analysis. The interfacial energy distribution in the middle cell samples appears to be more evenly distributed along the depth direction than in the top, bottom, and the reference cells. By contrast, the reference sample appears to have its free energy distribution skewed towards the greater depths of the sample. The preand post-degradation energy density comparison is supplemented by the decomposition of the free energy into its contributors: the interfacial energy of the Ni-YSZ interface, the interfacial energy of the Ni-Pore interface, and the interfacial energy of the YSZ-Pore interface. The results are shown in Figure 8. It can be seen that the relative contribution of the Ni-YSZ and especially the Ni-Pore interfaces decreases for the post-degradation cells. This is an expected result, since the nickel coarsening and agglomeration is the main mechanism for the microstructural free energy reduction. However, at the same time, the YSZ-Pore interfacial energy increases-as the nickel particles move and change shape, more YSZ surface is uncovered. This process is particularly apparent for the near-inlet sample of the top cell, where the increase of the YSZ-Pore interfacial energy nearly compensates for the loss of the total free energy. More decay of F is observed in the near-TPB elements, as a result of the Triple Phase boundary deterioration (previously discussed by Brus et al. [2]).The decline of F NiPore interfacial energy was most pronounced in the top cell samples, and least pronounced in the bottom cell, suggesting some relation to the temperature distribution.

Conclusions
In this paper, the microstructure data from a stack recovered after a long-term performance experiment was analyzed by comparing the microstructural free energy functional values for microstructures investigated before, and after the durability study. Interestingly, the degradation rate does not seem to be connected to the local temperature. The middle cell, rather than the top cell, experienced the most total free energy decay, despite Nickel-Pore surfacial energy component deteriorating more for the hotter cells. Additionally, location in proximity to either inlet, or outlet appear to have some influence, as the downstream location in the bottom cell experienced more decay of the free energy. It is possible that the irreversibilities related to channel concentrations of hydrogen and steam contributed to the observed tendencies. Yet, the samples from the middle cell, the cell with the most free energy decay-did not exhibit greater deterioration for the downstream location. The greatest rate of degradation was observed for the cell which was previously shown to exhibit the highest degree of anisotropy.  Acknowledgments: Matplotlib Python library was used for data visualization [24]. AVIZO (Thermo Fisher Scientific) software was used for image processing. MATLAB (Mathworks) was used for some computations.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.
Sample Availability: Samples are available on reasonable from the authors.

Abbreviations
The following abbreviations are used in this manuscript: