Solid oxide fuel cell interconnect design optimization considering the thermal stresses

The mechanical failure of solid oxide fuel cell (SOFC) components may cause cracks with consequences such as gas leakage, structure instability and reduction of cell lifetime. A comprehensive 3D model of the thermal stresses of an anode-supported planar SOFC is presented in this work. The main objective of this paper is to get an interconnect optimized design by evaluating the thermal stresses of an anode-supported SOFC for different designs, which would be a new criterion for interconnect design. The model incorporates the momentum, mass, heat, ion and electron transport, as well as steady-state mechanics. Heat from methane steam reforming and water–gas shift reaction were considered in our model. The results examine the relationship between the interconnect structures and thermal stresses in SOFC at certain mechanical properties. A wider interconnect of the anode side lowers the stress obviously. The simulation results also indicate that thermal stress of coflow design is smaller than that of counterflow, corresponding to the temperature distribution. This study shows that it is possible to design interconnects for an optimum thermal stress performance of the cell.


Introduction
The solid oxide fuel cell (SOFC) is considered to be one of the most promising new energy technologies for high energy efficiency, environmental friendly and fuel diversity [1,2]. The performance of an SOFC is strongly related to properties of the electrode materials, seals or interconnect materials and impurities in the fuel [3][4][5][6][7]. During the past decades, there have been significant research aiming to increase the performance of SOFC by exploring new materials or improving synthesis techniques and optimizing the cell structure [8][9][10][11]. However, the SOFC technology still faces many challenges, before it could be commercialized, since cost reduction and an extended lifetime are required.
Numerical simulation as an economic approach to design materials, cell or stack structure was used to develop the SOFC technology further. The effect of microstructure on effective ionic and electrical electrode conductivities was investigated by numerical optimization, which determines the effective conductivities and degradation of materials [12]. The mass and heat transport, and electrochemical reaction were considered for SOFC optimization by developing numerical models [13][14][15]. The mixing ratio of fuel, gas flow rate and reforming reaction, utilizing various fuels were studied to archive an optimized SOFC system. Beside the material structure, heat and mass transfer processes and electrochemical problems, the thermo-mechanical phenomena attracted significant attentions in recent years [16][17][18]. As the most important part of thermo-mechanical phenomena, thermal stresses occurring to various components in SOFC strongly impact the cell lifetime by inducing cracks, contact loss and structural instability [19].
Thermal stress analysis of a single cell can be used to evaluate the possibility of cracks or flaws that was not able to be detected earlier and easier in experimental works [20]. For a planar anode-supported SOFC, thermal stress and thermal fluid behavior can be analyzed using a 3D integrated numerical model [21]. Maximum principal stress within the cell increases with a high current and temperature gradient, which was in good agreement with the experimental work. Since the flow distribution also had influence on the temperature distribution, thermal stresses of planar SOFC operating at coflow and counterflow were compared [22]. The results indicates that the temperature gradient near the fuel inlet for counter-flow pattern is much larger compared to that of co-flow pattern. Thermal stress in an single cell is different compared to the stress distribution in the stack, i.e., future studies are needed. Jiang et al. [23] built a thermo-electrochemical-structure model combining both the finite-volume and the finite-element approaches to investigate the thermal behavior and the thermal-stress of a SOFC with the bonded compliant seal design. It was found that, the thermal stresses depended on the location and the cell voltage. The second largest stress region within the cell was near the inlet. Except the effects of materials properties, the operating conditions such as the load of cell in a stack should be investigated for practical modeling. Effects of the applied assembly load on the thermal stress distribution in the integrated planar SOFC stack with a compressive sealing design were characterized by a 3D multi-cell model [24]. The results showed that the expansion mismatch rather than the applied compressive load dominated the thermal stress distribution of the cell components. There was also evidence that ceramic components suffered significant stress when subjected to an idealized operating duty cycle [25]. Stresses were generated due to differential thermal expansion of the layers indicating that there is a high probability of failure during these phases of the duty cycle, i.e. when cell is cooled from sintering to room temperature or is heated from room temperature to operating temperature. However, the just mentioned studies mainly focused on the thermal stresses distribution, considering the multiphysics process or operating conditions but neglecting the interconnect structural effects. As a component for connecting the electrodes and loads, the interconnect is not only essential to provide paths for electron transport in single cell (or to the neighboring cell in stack) and to protect the electrode material from damages in an ambient environment, but can also maintain the thermal mechanical properties matching to the adjacent electro-active components. Therefore, the selection criteria of interconnect materials and optimization of the structure is more important and stringent than other components of the cell. The ceramic interconnect possesses a good stability and retains a fine compatibility with other components. However, the conductivity of ceramics is not appreciable below 600°C and the poor sinterability of ceramic interconnects is also a challenge for implementation [26]. While, the metallic alloys are low cost, ease of fabricating and high mechanical strength, made it more attractive than ceramic oxides as interconnects in intermediate temperature SOFC stacks [27].
Normally, the interconnect optimization focused on the electrical performance, degradation processes or temperature distribution, since the structure and the shape of interconnects were related to concentration polarization and electric resistance [28][29][30]. However, the thermal stress depends on the temperature distribution and the load conditions of a complete single cell would be affected by the interconnect structure. With contact area increasing, power density and temperature gradient would be reduced when a decreasing size of collecting pins tended to gain a better temperature homogeneity and power density [31]. Beside the contact area and the size of collecting pins, addition of interconnect ribs is another option for modifying interconnect structure. While increasing the number of interconnect ribs and reducing the gap between them at the cathode side, the lateral conduction distance can be decreased, and an enhancement of the cell performance of more than 30 % can be achieved [32].
Besides, altering the width of ribs is easier both in numerical and experimental design. 3D numerical model was built to analyze the effect of rib width on cell performance [33]. The relationship, between the contact resistance and the optimal rib widths was given as a guidance for an engineering optimization. When considering the size of ribs, the increased width gives a better conduction of the electrical current and reduced ohmic losses, while narrow ribs are needed to facilitate a more uniform distribution of reacting gases and thus promoting the electrochemical performance. This implicates that a tradeoff of rib size to the cell performance is very significant.
In this work, in order to optimize the interconnect structure, based on calculations of the thermal stress field, a 3D comprehensive bi-layer model was developed. Elaborating the reliable strategies for the optimization coupling electrochemical and mechanical properties is of crucial importance for new SOFC designs. Different cell configurations with a geometry design based on the finite element method were applied in this model. The influences from the shape and size of interconnect as well as the components for SOFC involved the thermal stress distribution were investigated in terms of ion/electronic, momentum, mass and heat transport.

Mathematical model
Conservation and constitution laws were applied to each domain to obtain a comprehensive model containing all the above mentioned phenomena. Stress analysis of an SOFC was performed using the commercial software COMSOL Multiphysics (version 5.0). A half-cell model with bipolar channels operating with humidified hydrogen, carbon monoxide and methane as fuel. The material structure parameters and geometry were selected based on a cell developed and tested at Ningbo Institute of Material Technology & Engineering (NIMTE) [34]. The geometry was built based on the interfacial zone or layer thickness in those experiments [35,36] and the interfacial layers in the model were assumed to be 15 and 20 lm thickness for active anode and cathode layer and 0.1 m length for the cell, as shown in Fig. 1. Note that, the coflow case was used as a demonstration for our model geometry. The parameters of the geometry are presented in Table 1. Besides, two different interconnects designs are also shown in Fig. 1. The design 1 shows symmetrical interconnects with small contact area with the electrodes, while shape 2 has large contact area with electrodes.
4 Ion and electron transport For the ion and electron transport processes, the critical voltage distribution can be described as the potential difference between the anode and cathode current collectors [37]. The operating cell voltage (E) is reduced due to the internal resistance and polarizations, as shown in Eq. (6) Here, E OCV is the open-circuit voltage, g act , g ohm , g conc the activation, ohmic and concentration overpotentials respectively.
Here, / is the potential and E eq the equilibrium voltage. The index ''a'' and ''c'' stand for the anode and cathode, respectively. The electrode materials are noted as ''s'' and the electrolyte material as ''l''. When a hydrogen-steam mixture is used as fuel, the reversible open-circuit voltage, E OCV in Eq. (6) corresponds to the Gibb's free enthalpy. DG f;x;T is the Gibb's free enthalpy of the oxidation reaction of hydrogen [38]: The impact from carbon monoxide and hydrocarbons should be taken into account for a complete model. However, Eq. (10) is sufficient for the study performed in this work with focus on thermal stress.
The activation polarization/current density relationship is described by the Bulter-Volmer equation: Here, p is the partial pressure at the triple phase boundary (TPB) or bulk fluid within the gas channels (b) and n e is the number of electrons transferred per reaction. In the high current region with sufficiently large g a where exp an e;a Fg a RT Equation (11) can be rewritten as The activation polarization is calculated by using the exchange current density. For the electron and ion transfer problem in the anode and cathode, the charge conservation equation and the governing equations for the ion and electron transport by Ohm's law are implemented as following: where r is the ion/electron conductivity and / is the potential, i s is the total volumetric current density. In the anode, i s is the volumetric local current densities of electrochemical reactions of fuel. The electronic conductivity in the electrodes (r s , within the anode and the cathode) and ionic conductivity in the electrolyte (r l ) are calculated as described in Ref. [39]. The actual length that ions and electrons are transported in the electrodes is increased when the effects of the microscopic structures in the material are considered. This is accounted for by using the effective conductivities of each material which were determined by the structure-dependent tortuosities and the volume fractions/porosities [40].
As a medium for ion and transfer from the cathode to anode, the charge conservation in electrolyte can be written as where Q l is the source term, in the electrolyte which is equal to zero because there is no production or consumption of ions in electrolyte. When the heat source Q h for the ion/charge transport is taken into account, the energy conservation equation from the Joule effect can be expressed as 5 Momentum, mass and heat transport The temperature distributions depends on the flow of gases in the channels and porosity of the electrodes that are related to the momentum, mass and heat transfer processes. The momentum transport in the porous materials as well as in the fuel and air channel are solved simultaneously [41], as Here, F is the volume force vector, p the pressure, k the permeability of the porous material, ũ the velocity vector, e the porosity and w the viscous stress tensor. The viscosity (l) and density (q) for the gas mixtures are dependent on local temperature and the mole fractions and are calculated as described in Ref. [42]: where the b k is the species dependent parameter and ''k'' stands for the number of species dependent parameters in the viscosity equation. x i mole fraction of species, M i the molecular weight of species i. The mass balance for a gas species, i, can be expressed as in which J i is the mass diffusion flux and x i is the mass fraction, R i is the extra source term for production and consumption of species. Note that, R i ¼ 0 in the gas channels.
v is the stoichiometric coefficient for species, i is the volumetric current density for the electrochemical reaction, and n i is the number of participating electrons in the electrochemical reaction.
In this work, the Maxwell-Stefan model is used because of its simplicity and good accuracy: D T i is the thermal diffusion coefficient which is neglected due to its very small effect here.
A concentration gradient in the porous electrodes is considered when the effects of mass transfer under normal operating conditions for the redox reaction are calculated. In the porous medium, the molecular diffusion and Knudsen diffusion are commonly used to describe the mass diffusion mechanisms. The molecular diffusion is used to describe the situation when the pore size is significantly larger than the mean free path of the gas molecule diffusion [43]. When both Knudsen and molecular diffusion are considered, the effective diffusion coefficients can be calculated as Here, e is the porosity and s is the effective tortuosity factor. The effective diffusivity coefficients are calculated assuming electrodes with corrected factor e s ¼ 0:16. D k;ij is the Knudsen diffusion coefficient of the component i with the component j in a gas mixture, which is defined as [44] where r e is the effective pore radius.
where M i ; M j are the molecular weight of species i, j, respectively.
Heat generation within the cell from the electrochemical reactions (electron and ion transport as well as electrode reactions) and the internal steam reforming reaction within the anode are included as source term in the model [45]: The temperature is assumed to be locally equal for the solid material and the gases, based on the minor temperature difference between the gas-and the solidphase. The energy conservation equation can be written as Here, Q h is the heat generation or consumption, k eff is the effective thermal conductivity and c p,g is the gas-phase specific heat. The heat source term Q h , due to electrochemical reactions, ion and electron transport and losses through the activation can be defined as where DS r is the reaction entropy change, r ref is the reforming reaction rate and DH ref is the enthalpy change of the reforming reactions.

Thermal stress couplings
The materials used in this model are the ones most frequently used for SOFC. Table 2 shows material mechanical properties, which are obtained from literatures [39,43,46]. It is noted that the mechanical properties are determined by the content of materials as well as the synthesis process, which is the reason for the difference in properties listed in Table 2.
The material properties of the interfacial layers are assumed to be between those of the electrode and the electrolyte respectively. All of the materials are assumed to behave as linear elastic and isotropic. It should be noted that the boundaries of this single cell model are assumed to be free, i.e., load effects are not considered here. SOFC components expand with temperature, causing thermal strains to develop in the material when the deformation is constrained. The overall strain results from the summation of elastic and thermal stresses, i.e., the initial contributions are neglected: The form of strain as follow e el ¼ e xx ; e yy ; e zz ; c yz ; c xz ; c xy À Á ; here e xx , e yy , e zz , c yz , c xz , c xy are the longitudinal and shear components for strain, respectively. The thermal strain depends on the temperature, T, the stress-free reference temperature, T ref , and the coefficient of thermal expansion (CTE), a: Determining the stress free temperature is critical as it directly affects the magnitude of the thermal stress induced in the material. For SOFCs it is widely accepted that T ref is the sintering temperature, at which different layers are joined.
The stress-strain relationship for the linear material was calculated as where r 0 is the initial stress, which is treated as the residual stress in model. The elasticity matrix (D) for isotopic material is defined as, .0436) were used as fuel in this model. Note that the ion and electron transport was considered as a testing part in our model, i.e., the H 2 mixture and not the 30 % pre-reformed natural mixture was considered in this section, which is sufficient for the study on thermal stress. The air inlet is defined as air, including oxygen and nitrogen. The boundary conditions for the outlets are defined as convective fluxes. The inlet gas temperature is defined by the operating temperature (1,000 K) and the outlet is defined as a convective flux. The boundaries at the top and the bottom walls of the cell are defined as the symmetry conditions, because it is assumed that the cell is surrounded by other ones with identical temperature distribution. Also the boundaries at the side walls are defined as symmetry conditions, because the channel section is assumed to be surrounded by identical channel sections. The potential at the anode current collector is set to zero and the one at the cathode current collector as the cell operating voltage (0.7 V). All other boundaries and interfaces are electrically insulated. The governing equations are segregated in five different groups: 1. Velocity field, pressure distribution, and pressure corrections. 2. Temperature distribution. 3. Ion and electron distribution. 4. Mass fraction distribution on the air side and fuel side. 5. Thermal stress distribution.
The segregated solver is applied for 7,136,654 degrees of freedom and the solution tolerance is defined to 0.001 for each segregated group. The calculation time is around 80 h on a single computer with 32 GB RAM and a CPU with 4 GHz. Note that it is hard to give an exact value for the calculation time since the model is built in several steps, where each step starts its calculation from the previous one.

Results and discussion
Note that the color legend and the operating parameters are kept the same for all cases. The bi-layer SOFC model geometry for a half-cell (symmetry exists in the single cell) includes interconnects with the other components of a SOFC, as illustrated in Fig. 1a. The gas flow direction of coflow case is pointed out with an arrow, while the fuel flow is in the opposite direction for the counterflow case. The experimental data from Kyushu University are used for some adjustment of the simulation conditions as well as validation for the model results, as shown in our previous work [45]. For optimization of the shape, two types were designed to investigate the support and contact effects of interconnects, as shown in Fig. 1b. The rib in shape 1 has an inverted trapezoidal, with the same width of ribs contacting to the PEN structure (support positive layer-active positive layer-electrolyte-support negative layer-active negative layer) as the standard case, but width of ribs linked to the top interconnects increased by 0.1 lm. With the shape 1, the trapezoidal ribs that contact to the PEN structure is wider than the standard case (0.1 lm), as shown in shape 2. Those two types were designed, considering the possibility that the contact area and the mechanical strength of the structure would affect the thermal stress distribution.
The temperature distribution for the counterflow case of standard shape is presented in Fig. 2. The temperature increase along the main flow direction is caused by the strongly exothermic electrochemical reactions as well as from the different polarizations. Note that the temperature decrease since the endothermic reforming reactions are considered. The maximum temperature (1,124 K) is located in the inside of the cell near to the position the fuel inlet (0.04 m to the fuel inlet).
Coupling the mass and heat transport, reforming reaction, solid mechanical and the first principal stress of SOFC with different designs are shown in Fig. 3. The positive values tell that tensile stress and negative value for compressive stress. The active electrodes and electrolyte component involved in a tensile stress while the support electrode for a compressive stress. This stress difference results from the properties mismatch and it may lead to curling of the components. It is noted that maximum tensile and compressive stresses for coflow are smaller than that for counterflow, corresponding to the temperature distribution. The greater temperature gradient in counterflow also results in an obviously tensile stress for interconnect at fuel inlet, as seen in Fig. 3a. Similar stress distribution can be seen in shape 1 and 2, shown in Fig. 3c and d. The stress decreases with the new interconnect shape design, and the shape 2 is slightly smaller than shape 1, which may be due to the subtle difference in temperature distribution of heat conduction for the two interconnect structures.
The effects of mismatch mechanical properties can be analyzed by comparing the stress distribution at each components' interface. Since the mismatch between the interconnect and the cathode is larger than for the anode (Table 2), the thermal stress at cathode side is larger than at anode side, as shown in Fig. 4. It should be noted that, the horizontal direction represents the z axis in the model (direction vertical to gas flow) and the vertical direction represents the x axis (direction parallel to gas flow). The palliation effects of thermal stress of shape 1 and 2 design was obvious. The interconnect shape design is more effective to reduce the tensile stress in the PEN structure. However, the thermal stress distribution would be extended in the shape 2 design case, because the interconnect has a higher heat conduction compared to the gas, leading to a higher temperature on the surfaces of support layers.
The thermal stress on the interface close to the inlet and outlet for the different designs is shown in Fig. 5. Since the temperature gradient at air inlet is smaller, a tiny change of stress can be seen in Fig. 5a, c, e. The slight increase of stress may be due to the slightly increased temperature. However, the tensile stress of the interconnect at the anode side was obviously reduced, and this decrease can also be seen at the active layers and in the electrolyte, compared to the fuel inlet in Fig. 5b, d, f. An increased compressive stress can also be observed within the support layers. Figure 6 shows that the stress has a relationship with the location, i.e. for different components. For the interconnect at the cathode side, the stress at the air inlet is bigger than that at the fuel inlet, and the minimum stress locates at the middle position of the cell. However, the maximum stress belongs to the fuel inlet for the interconnect at the anode side, and stress at air inlet is close but slightly higher than that at the middle position of cell. When considering the thermal stress distribution for the PEN structure, the middle position dominates at the electrolyte while the fuel inlet and air inlet are comparatively large at the cathode and anode support layers, as depicted in the embedded figure.
It is observed in Fig. 7 that the thermal stress is almost the same when the height of rib at anode side is increased, but causes a rise at the interconnect close to the cathode side, which is illustrated in the figure embed in Fig. 7. This may be explained by the symmetry stress effect of structure. Figure 8 shows that the strong influence of the rib width on the thermal stress. Wider ribs and ribs covering bigger fraction of the cell may reduce the interface resistance to current flow by increasing the electrode-interconnect contact area and reducing the current path through the possibly high resistance electrode material. However, the gas channel would be narrowed, as a consequence, the electrochemistry reaction is affected. To illustrate, considering the three designs in Fig. 8, wider interconnect at the anode side would lower the stress. Conversely, this may be useless for reducing the stress at the cathode side.
The thermal stresses of the coflow case for an electrolyte with a wider rib shown in Fig. 9. One can see that the stress obviously rise with wider interconnect ribs. The stress changes increases along the flow because of the rise of temperature. This difference may occur due to the superposition of interconnect and electrodes side (ribs) changes, which affect gas concentration at the three-phase boundary (TPB) region.

Conclusions
Interconnects are critical to minimize the overall stack resistance and weight as well as to enhance the stability when thermal stresses are considered. A comprehensive model to investigate the relationship between the cell structure, temperature distribution and thermal stresses in SOFCs at certain mechanical properties was built using the finite element method. Proper design of interconnect in conjunction with single-cells were implemented based on the thermal stress optimization performed in this work.
The palliation of thermal stress of different designs was obvious. The interconnect shape design was the most efficient to reduce the tensile stress in PEN structure. Besides, thicker interconnect ribs causes a rise at the  interconnect close to cathode side while this effect for thermal stress at the anode side is smaller. It is also observed that, wider interconnect of anode side lowers the stress. We show, by numerical calculations, that the thermal stress can be reduced through structural correlations. This optimization can be used as a new strategy for interconnect optimization in SOFCs which is beyond single criterion for electrochemical performance.  Mole fraction of species j (mol/mol) Greek symbols a Coefficient of thermal expansion (K -1 ) b Transfer coefficient e Porosity e el Elastic stresses contribution e th Thermal strain contribution g Over potential (or polarization) (V) U Potential (V) k Permeability (m 2 ) l Dynamic viscosity (Pa s) q Density (kg/m 3 ) r l,s Ionic/electronic conductivity (X -1 m -1 ) r Stress v Poisson's ratio s Tortuosity D Numerical difference