Meso-Scale Phase-Field Modeling of Microstructural Evolution in Solid Oxide Fuel Cells

Microstructural evolution occurs in solid oxide fuel cells (SOFCs) during operation, which cause severe electrochemical performance degradation as well as structural failure such as crack formation. Nickel (Ni) particle coarsening is believed to cause microstructural evolution in the Ni-based anode of SOFCs. Furthermore, the accumulation of expanded pores during the microstructural evolution is responsible for crack formation. Based on the diffuse-interface theory and the phase-ﬁeld method, integrating both Cahn-Hilliard equation and Ginzburg-Landau equation, a meso-scale phase-ﬁeld model was established to investigate microstructural evolution and crack formation in the Ni / Yttria Stabilized Zirconia (YSZ) anode of SOFCs. Then, the model was applied to explore the coupling effects of microstructural evolution on SOFCs performance degradation. Simulation results show that particle size, porosity, and particle size ratio of Ni-YSZ are most inﬂuential parameters in microstructural evolution. It was foundthat the triple phase boundaries (TPBs) area was decreased by 24% and power density was decreased by 11.03% after 1000 hours of operation. In addition, the power density was decreased by 27%. This work is expected to provide us a comprehensive understanding of SOFCs’ microstructural evolution as well as a tool for SOFC’s performance degradation analysis. of Ni particle coarsening in electrode microstructure. Randomly gen- erated microstructure was considered and the time step size of 1 hour was used for the simulation. Developed phase-ﬁeld model was vali- dated against the experimental work by Tanasini et al. 24 the modeling work can be used to predict the morphology of three phases associ- ated with the electrode. The volume fraction of the phases is used to ﬁnd out the TPB area reduction. Statistical material properties due to microstructural changes obtained from the phase-ﬁeld simulation results are used as input parameter to the previously developed cell level model to quantify the performance degradation. The simulation result shows that the average particle size is in- creased during temporal microstructural evolution. The Ni particle size is increased by 26%, which leads to the TPB area reduction. Moreover, TPB phase contiguity is also affected signiﬁcantly due to 24% reduction if Ni-YSZ particle size. Temporal evolution of the microstructure leads to a decrease in the TPB area by 24%, which has a profound effect on the microstructural performance level. The power output is found to be decreased by 11.03% after 1000 hours of operation. Additionally, crack formation due to the accumulation of pores during phase diffusion is also demonstrated. This work is expected to provide us a comprehensive understanding of SOFCs’ microstructural evolution as well as a tool for SOFC’s performance degradation analysis.

Solid oxide fuel cells (SOFCs) could be one possible solution to the depleting energy resources of the modern world due to its fuel flexibility, low carbon emissions, and high efficiency. With up to 60% energy efficiency and a life expectancy of 40,000 hours, the SOFC has emerged as an ideal candidate to meet the energy challenges of the modern world. [1][2][3] However, the commercialization of SOFC is hindered due to its high operating temperature, manufacturing cost, and structural instability. [4][5][6] Structural evolution of SOFC often leads to the structural failure and shortening of SOFC lifetime. 7,8 Moreover, structural evolution results in compromising SOFCs' electrochemical performance, mostly in the electrodes.
A SOFC consists of both electrodes (i.e., anode and cathode) and a ceramic electrolyte. The anode is consisted of ion conducting phase and electron conducting phase to facilitate the fuel oxidation reaction. 9 The cathode reduces the oxygen into oxygen ion, which is then diffused to anode through the ceramic electrolyte, such as yttria stabilized zirconia (YSZ). 10 Electrochemical reaction of SOFC mainly occurs in the triple phase boundary (TPB) area, 11 where the pore phase, the ion conducting phase, and electron conducting phase join to convert chemical energy of fuel gases to electrical energy. 1,8,12,13 TPB area is the key microstructural property controlling the electrochemical activities of SOFC. 14 It is desired to have a stable microstructure with optimized TPB area. 15,16 However, recent studies reveal that TPB area is often diminished due to microstructure evolution. 10,[17][18][19][20] The microstructure evolution is controlled by the Ni particle coarsening, which in turn leads to the performance degradation and crack formation. Particle coarsening process was first reported by Ostwald et al. 21 In detail, Ostwald ripening refers to the late stage of phase separation. During the phase separation, either droplets or particles evolve to minimize the surface free energy via minimizing their surface-to-volume ratio, which occurs by different mechanisms such as condensationevaporation mechanism. 22 However, Gibbs-Thompson effect is often used to explain the particle coarsening phenomena 23 in the solid. Various experimental investigations explored the effects of Ni particle coarsening on the performance degradation of SOFC. For example, Tanasini et al. 24 revealed the performance degradation is directly related to Ni particle coarsening in the anode electrode. Experimental study by Simwonis et al. 25 showed that electrical conductivity of anode decreases due to Ni coarsening. Nelson et al. 26 reported TPB contiguity is affected by Ni particle coarsening. Faes et al. 27 investigated the interplay between particle size and TPB by measuring Ni-YSZ particle size and the TPB area using microscopy techniques. Experimental investigation by Iwata et al. 28 showed that the pore radius decreases due to Ni particle coarsening. Recent work by Cronin et al. 29 demonstrated that both pore percolation and TPB area are affected by Ni particle coarsening process.
Those aforementioned experimental investigations provide insight into the SOFC degradation process caused by Ni particle coarsening, however, the numerical model is required for theoretical explanation of experimentally observed phenomena as well as evaluating SOFCs' performance in the long term. In this modeling work, the coarsening of Ni particle is considered as a capillarity driven phenomena. 30 According to the Gibbs-Thomson effect, materials with high curvature regions have higher chemical potential than that of the lower curvature regions, this difference in chemical potential controls the transportation of materials from higher curvature regions to lower curvature regions this process is known as minimization of total free energy. In SOFCs, the modeling of the microstructural evolution is dependent on local interface mean curvature of Ni particle, which is driven by the process of minimization of total free energy of the system. 25 The total free energy function is governed by the dominant transport mechanism (i.e., surface-and volume-diffusion) of Ni particle. 31 Additionally, high operating temperature of SOFC leads to the enhanced mass transport of Ni particle. 32 During its operation, the surface diffusivity (or mobility) of Ni particle is much higher than that of the YSZ particle, which essentially makes YSZ stationary or non-evolving compared to the Ni particle, [33][34][35] serving as the supporting structure in which Ni coarsens and evolves via diffusion along Ni-pore interfaces. Moreover, the evolution of the microstructure leads to the pore size growth and the accumulation of expanded pores during the surfaceand volume-diffusion, which causes crack formation. 36 However, it is quite challenging capturing microstructural evolution due to explicit tracking of the interface moving boundaries.
To overcome the aforementioned issue, researchers have been applying the phase-field method for illustrating the variety of interfacial phenomena among different immiscible interfaces, microstructural evolution, and crack formation. 30,37,38 The phase-field method, which allows explicit tracking of the phase interfaces, has emerged as a powerful computational method for modeling transport phenomena and complex microstructural evolution process. 9 In addition, the modeling approach allows visualizing the microstructural evolution in two-and three dimensional. 9 The phase-field method separates an immiscible binary mixture into a domain pure in each component. It applies different orientation fields to defining the crystallographic orientations of Ni-and YSZ phases. The interfacial dynamics of the anode is defined by a variable called Order Parameter (OP) (i.e., phase-field variable). For instance, an OP C = 0 to 1 can be used to indicate the domains of different constituents' particles, i.e., C = 1 indicates one particle phase while C = 0 denotes the other particle phase in an immiscible binary mixture. The interfacial region is tracked by 0 < C < 1. The phase-field modeling approach has been successfully applied to explore the Ni particle coarsening associated with SOFCs evolution. For example, modeling work by Li et al. 9 revealed that TPB area is significantly reduced due to the microstructural evolution of SOFC. Chen et al. 30 quantitatively characterized the microstructural evolution by examining the TPB density, interfacial area, and tortuosity versus time using phase-field method. Abdeljawad et al. 19 established an integrated phase-field model coupling the analysis of the microstructural coarsening process with SOFC electrochemical performance evaluation of SOFC. Many researchers, including but not limited to Aranson et al., 39 Karma et al., 40 Mieche et al., 41 and Marconi et al. 42 have also explored crack formation using the phase-field model. Despite previous efforts, further investigation is yet to be done on long-term degradation of SOFC caused by microstructural evolution, which has been identified as one of the major sources of SOFCs structural failure and performance degradation. 19,30 Experimental investigation of SOFC electrodes degradation needs more further theoretical analysis. 24 In particular, a more comprehensive understanding is needed to investigate the interplay between the dynamics of microstructural evolution at the meso-scale and electrochemical performance at the macro-scale.
In this work, a meso-scale computational framework was established to investigate the quantitative effect of Ni particle coarsening on the performance of SOFC. An integrated phase-field model was developed to couple microstructure evolution and crack formation in the anode. The developed meso-scale model captures the dynamics of Ni particle coarsening and crack formation in Ni-YSZ anodes. Moreover, it tracks the evolution of microstructural features that influence the electrochemical performance. The time dependent microstructural properties obtained from the phase-field model was applied as effective properties to our previously developed SOFC electrochemical model 1 to investigate the electrochemical performance of SOFC.
A graphical representation of performance degradation induced by SOFC microstructural evolution can be referred to Figure 1.

Model Development
In this work, SOFC anode consists of three constituents, Ni particle, YSZ, and pore. The mass of each phase inside anode is assumed to be unchanged. The modeling approach is based on the minimization of the total free energy of the anode that is represented by F tot . YSZ is assumed to be stationary due to slow diffusion compared to Ni that evolves via diffusion along Ni-pore interface. Three OPs are utilized to represent Ni phase (ζ), YSZ phase (ξ), and crack (ϑ), respectively. Based on the phase-field method, an OP changes continuously from 0 to 1 based on diffuse-interface theory. The range over which the OP changes from 0 to 1 is defined as interfacial region or interfacial thickness. For example, ϑ = 1, ϑ = 0, and 0 < ϑ < 1 represent the intact, fully broken, and the transitional interfacial thickness of the material. Minimization of total free energy leads to the interfacial diffusion. Interfacial diffusion can be formulated by using the Cahn Hilliard equation as follows, where, M is the mobility function that depends on the OPs of the system, ∇ 2 is the Laplace operator, and μ is the chemical potential of the system. The chemical potential is driven by total free energy that is based on dominant transport mechanism. The three phases of Ni-YSZ anode is considered immiscible, therefore, the effect of elastic energy is assumed to be negligible. F tot represent the total free energy of anode and can be expressed as follows, where the first term on the right represents gradient energy, ϕ donates interfacial thickness, f (ϑ) is the fracture OP of anode, and f (σ) is the elastic strain energy. f(ζ ) is a generic double wall function that can be represented as follows, Two boundary conditions (i.e., contact angle at the TPB and no flux at the Ni-YSZ interface) are implemented with the Cahn-Hilliard equation. No flux BC (i.e., ∂μ × ∂ζ = 0) can be solved using the smoothed boundary method (SBM). 43,44 The BC of contact angle is implemented as elaborated in the literatures. 45,46 The Cahn-Hilliard equation for anode within SBM framework can be written as follows, where M is the mobility function of two OPs (i.e., ζ and ξ), κ = M(ζ, ξ)∇μ.
The chemical potential of anode is denoted by μ and can be represented by the following equation, 36 where ϕ represents interfacial thickness, g(ϑ) = ϑ 3 (4 − 3ϑ), is the coupling function between the fracture OP and elastic stress field. Karma et al. 40 showed the choice of g(ϑ) does not affect the transport mechanism as long as g(ϑ) = 0. ϑ = 0 represents the material is fully fractured. E is the Young's modulus, υ is the Poisson's ratio, and is the volume fraction of Ni phase (i.e., ratio of Ni volume to total anode volume).
Multiplying with the OP ξ, Eq. 4 can be rewritten as follows, Applying the no flux boundary condition, i.e., ∂μ × ∂ζ = 0, into Eq. 6, The unit normal vector of YSZ particle can be formulated as follows, n = ∇ξ |∇ξ| [8] Immiscible properties of Ni-YSZ controls the crystallographic orientation of the structure and is used to formulate the TPB contact angle θ, where the negative sign is due to the outward normal of the Ni particles.
The energy equilibrium at the TPB corresponds to an extremum of the total free energy, i.e., ∂ F tot = 0. Hence Eq. 2 can be written as follows, Eq. 9 can be formulated as follows, The microstructural evolution governing equation can be defined as follows, In Eq. 12, M is the mobility function and can be formulated as follows, [13] where, g(ξ) is a function to control the Ni particle mobility near the TPB area, can be represented by a polynomial function. 30 The pores in anode are treated as pre-existing crack or initial cracks. Following the previous work, 36 the governing equation for crack formation using fracture OP can be described as follows, [15] where, ω is the relaxation constant 47 of fracture order parameter, α is gradient energy co-efficient, H is the energy barrier, 48 σ st is the elastic energy density, and σ is the threshold value strain. 36 The microstructural properties are expressed in non-dimensional units in order to allow the investigation with a wide range of particle size and particle size ratio. Suitable parameters of the mobility function are necessary for the effective calculation of the Ni particle coarsening process. Due to diffuse interface nature of the microstructure, it is necessary to get an asymptotic analysis of phase-field method for the normalization of time. Following the references, 49,50 an asymptotic analysis has been developed to measure the characteristic time (time normalization), where, L represents characteristic length scale of sample, k B is the Boltzmann constant, T is the temperature, N v is the atomic number density, δ s is surface diffusin depth, γ s is surface energy, and D s is surface diffusivity.

Numerical Setup and Solver Flow
The simulations were carried out by solving governing Eqs. 12-15 in FEniCS, an open source finite element software. Numerical calculations were on the order of 884736 elements and required about 16 GB of RAM and 140 hrs of time when solved in parallel over seven cores on an Ubuntu Linux workstation with 24 3.4 GHz Intel Xeon processors. The numerical tolerance for convergence and tolerance for the Newton solver were both set to be 1 × 10 −6 . The Newton iterative solver was utilized due to nonlinear nature of the governing equations. An open-source multiplatform tool ParaView was used for data analysis and visualization. The numerical solver flow can be referred to Figure 2.

Model Validation
The work focuses on the modeling of microstructural evolution in SOFC anode. Randomly disordered Ni-YSZ anode was adopted for the simulations due to the difficulty of obtaining initial crystallographic orientation. Ni particle size of 0.55 μm and Ni-YSZ particle size ratio of 0.9 were used as initial values. Time step size for simulation represented one-hour operation in real time. The mobility for Cahn-Hilliard equation was taken as one. Simulation results show 26% increment in the Ni particle size after 1000 time steps, which agrees with the experimental result reported by Tanasini et al. 24 Additionally, the particle size ratio of Ni-YSZ was decreased by 16%, while Tanasini et al. 24 reported 14.6% decrement based on their experimental study. Figure 3 shows the comparison between the experimental and modeling work for both particle size and particle size ratio. Comparison of TPB evolution shows a close match with the literature from Li et al. 9 As for the volume ratios of three different phases, we used 30% Ni phase, 30% YSZ phase, and 40% porosity to validate the simulation result obtained from the modeling work. Figure 4 shows close match between the two simulation results for TPB fraction (i.e., TPB fraction is defined as the ratio between TPB volume and anode volume). The TPB area is calculated using the ParaView. The key model parameters and values are listed in Table I.

Results and Discussion
In this work, the developed meso-scale phase-field model was applied to investigate the anode microstructural evolution and crack  formation. The SOFC anode microstructure evolution is driven by the Ni particle coarsening process, which leads to TPB area reduction, long term performance degradation, and crack formation of SOFC. The developed model was also used to obtain the time dependent microstructural properties during microstructure evolution, which is later used to quantify the performance degradation of SOFC.
The temporal microstructural evolution in two-dimensional is shown in Figure 5. A randomly disordered microstructure is generated in a normalized domain, where blue color represents the Ni particle, green color represents the YSZ, and red color is the pore,  to explore the microstructural evolution process. During SOFC operation, the anode evolves over time due to coarsening of Ni particle. Due to the nature of diffusive interface, YSZ-pore phase interface and Ni-pore interface can be tracked with obtaining their microstructural parameters (e.g., particle size ratio, TPB fraction) by solving governing Eqs. 12-15. This particle coarsening process is driven by the minimization of total free energy, which causes average particle size to increase and the percolation connectivity to decrease. Two different stages of microstructures at T = 40 (i.e., 40 hrs operation) and T = 1000 (i.e., 1000 hrs operation) are shown in the figure.
The temporal three dimensional structural simulations are completed for an anode microstructure of 25 μm × 25 μm × 25 μm with a mesh grid of 96 × 96 × 96. The Ni-YSZ microstructure was generated following a previous study, where three phase volume fraction was prescribed as 27% Ni phase, 30% YSZ phase, and 43% pore phase. 30 Figure 6 shows temporal three dimensional microstructural evolution for Ni, YSZ, and pore after 1000 hours of operation. Each of the constituents' phases is shown in Figure 6. The simulation results show the internal structure of the connectivity of each phase, which can be utilized to quantify the growth of the average characteristic length and discriminate the percolated cluster from isolated clusters in the evolved anode.
Ni particle coarsening decreases TPB density and increases the resistance of anode. Figure 7 visualizes the reduction of the TPB area during the microstructural evolution process. The result shows, rapid reduction of TPB area at the initial stage of the coarsening process, which agrees with the experimental observation where the main performance change occurs at the initial stages. Due to high interfacial energy at the initial stages, the TPB area decreases rapidly, but it slows down over the time period. The TPB area reduction is consistent with the temporal phase coarsening, which indicates the performance degradation is directly related to the particle agglomeration and coarsening of Ni particle. The overall TPB area reduction is found to be approximately 24%, which is close to the experimental work of Cronin et al. 29 The evolution of pore phase has been illustrated in Figure 8. Numerous pores exist in anode to facilitate gas transportation. Preexisting pores in Ni-YSZ grow due to interfacial diffusion during anode microstructure evolution which results in the increase of pore radius. By examining the evolution of pore phase, it can be concluded that the pores are growing in size and evolving into a more round shape, thus minimizing total free energy. The coupling effect of diffuse interface and stress evolution plays a vital role in the formation of cracks. Radius of the pore increases during the minimization of interfacial area. The interfacial diffusion during the evolution process leads to the accumulation of these expanded pores which causes crack formation. In the meanwhile, as Ni particle diffuses during microstructural evolution, the tensile stress field in the inner region evolves, which causes stress concentration at the crack tip. The stress concentration eventually leads to the formation of crack and subsequently structural failure. The simulation results are in good agreement with annealing experiments on Ni-YSZ anodes with Ni and YSZ particles of similar sizes, where a 25% decrease in total Ni-pore interfacial area was observed. 24 The temporal pore phase evolution is presented in Figure 9, which helps to visualize the crack formation during microstructural evolution in three dimension. Only the pore phase after the evolution is shown for an anode microstructure of 25 μm × 25 μm × 25 μm with a mesh grid of 96 × 96 × 96 to demonstrate crack formation. In the modeling work, the initial crack appears after 780 dimensionless time steps. Then, the crack formation develops quickly, as shown in Figure 9, at two different dimensionless time steps, T = 40 and T = 1000.
The microstructural features and effective material properties extracted from meso-scale phase-field simulations can be utilized to investigate the influence of Ni coarsening on the electrochemical performance degradation of SOFCs. For the electrochemical analysis in this section, we assigned an absolute value to the size of non-evolving YSZ particle as the characteristic length scale. In order to make the results as general as possible, all dimensions were scaled with respect to the characteristic length and no absolute physical dimensions of  the simulation box had been assigned. Peak power density plots are shown to quantify the performance degradation of SOFC due to microstructural evolution. Figure 10 compares the peak power density plot at the initial stage and after the degradation for an anode. The result shows, the peak density is 9550 W m −2 at the initial stage of the operation. However, the final peak power density decreases by 4.3% and reaches to 9150 W m −2 after the initial 500 hours of operation. After 1000 hours of operation the peak power density is further reduced to 8680 W m −2 . Overall, 1000 hours of operation results in 9.11% decreased power output. However, the performance degradation appears to slow down during 500-1000 simulation time steps. The performance degradation is caused by the combined effects of continuous loss of TPB area, decrease in electronic conductivity, and the crack formation. Figure 11 shows the performance degradation of anode supported SOFC due to crack formation. The formation of crack causes the interruption of current flow, which in turn increases the polarization resistance of the anode microstructure. The increased polarization resistance due to crack formation causes a sudden fall in the power density curve. 54,55 The increase of polarization resistance is obtained from the experimental study. 55 The peak power density at the initial stage of anode supported SOFC is 11345.34 W m −2 . But the formation of cracks causes increased resistance and results in a 27% reduction in the peak power density and the obtained peak power density is 8210 W m −2 .

Conclusions
A meso-scale phase-field model is developed to investigate performance degradation and crack formation during anode microstructural evolution. A phase-field model was developed to illustrate the effect Figure 11. Power density curve before and after crack formation for anode.
) unless CC License in place (see abstract). ecsdl.org/site/terms_use address. Redistribution subject to ECS terms of use (see 207.241.231.81 Downloaded on 2018-07-20 to IP F625 of Ni particle coarsening in electrode microstructure. Randomly generated microstructure was considered and the time step size of 1 hour was used for the simulation. Developed phase-field model was validated against the experimental work by Tanasini et al. 24 the modeling work can be used to predict the morphology of three phases associated with the electrode. The volume fraction of the phases is used to find out the TPB area reduction. Statistical material properties due to microstructural changes obtained from the phase-field simulation results are used as input parameter to the previously developed cell level model to quantify the performance degradation.
The simulation result shows that the average particle size is increased during temporal microstructural evolution. The Ni particle size is increased by 26%, which leads to the TPB area reduction. Moreover, TPB phase contiguity is also affected significantly due to 24% reduction if Ni-YSZ particle size. Temporal evolution of the microstructure leads to a decrease in the TPB area by 24%, which has a profound effect on the microstructural performance level. The power output is found to be decreased by 11.03% after 1000 hours of operation. Additionally, crack formation due to the accumulation of pores during phase diffusion is also demonstrated. This work is expected to provide us a comprehensive understanding of SOFCs' microstructural evolution as well as a tool for SOFC's performance degradation analysis.