Numerical Analysis of Thermal Stress for a Stack of Planar Solid Oxide Fuel Cells ^†

Abstract

In this work, a 3D multi-physics coupled model was developed to analyze the temperature and thermal stress distribution in a planar solid oxide fuel cell (SOFC) stack, and then the effects of different flow channels (co-flow, counter-flow and cross-flow) and electrolyte thickness were investigated. The simulation results indicate that the generated power is higher while the thermal stress is lower in the co-flow mode than those in the cross-flow mode. In the cross-flow mode, a gas inlet and outlet arrangement is proposed to increase current density by about 10%. The generated power of the stack increases with a thin electrolyte layer, but the temperature and its gradient of the stack also increase with increase of heat generation. The thermal stress for two typical sealing materials is also studied. The predicted results can be used for design and optimization of the stack structure to achieve lower stress and longer life.

1. Introduction

As a new power generation technology, a solid oxide fuel cell (SOFC) converts chemical energy directly into electrical energy through electrochemical reactions at a high temperature [1,2,3,4]. Compared with other new energy technology, SOFCs have obvious advantages, such as high thermal efficiency with good fuel flexibility [5]. However, in the operating temperature range of 700~1000 °C, cracks and structure damage may occur in SOFC stacks due to a mismatch of thermal expansion coefficients (CTE) of the various functional components involved [6,7,8,9,10]. In the operating process, high temperature gradients and thermal stress generated by the electrochemical reactions may lead to performance degradation and even failure [11,12,13,14]. Therefore, it is of great significance to analyze the maximum thermal stress and its distribution, as well as influencing factors.

At present, experimental determination of the thermal stress is limited to non-in situ testing by scanning electron microscopy and X-ray methods after SOFC stack operating, cooling and disassembly; however, the obtained results do not accurately reflect the real evolution process during the operation [15,16,17]. The in situ study of the thermal stress is mainly based on theoretical analysis, modeling and numerical simulation. There are some previous investigations analyzing the thermal stress distributions in the SOFC stack by considering the temperature gradient, the mismatch of CTE and the sealing material behavior, etc. [18,19,20,21]. Lin et al. developed a 3D model for a 3-cell SOFC stack, taking into account for the first time the complete structure including electrodes, interconnect, Ni grid and seal, aiming to analyze effects of the supporting as well as seal material properties on the potential failure locations of the SOFC stack based on the coupled CFD/FEM method [19]. The results show that the seals may be damaged by shear stress at the room temperature, and the mismatch in CTE of the components contributes more to the thermal stress than the temperature gradient. Chang et al. obtained the temperature distribution function and thermal stress characteristics of a planar SOFC stack through the simulation by star-CD commercial software [18]. Their results show that, under the higher current and temperature gradient, PEN (Positive/Electrolyte/Negative) can increase the maximum principal stress. Selimovic et al. established the electrochemical, thermal and structural coupling models for both the steady-state and transient thermal stress of SOFC stack. The location where the maximum thermal stress occurs with reforming methane as the fuel is identified, the thermal stress values at different heating rates are also predicted [20]. Xu et al. studied the influence of gas flow channel designs and external constraints on the thermal stress; it was found that the co-flow design arrangement can reduce the thermal stress, and the additional stress generated by a fixed external force may be one of the factors causing the stack fracture [21].

The effects of the assembly force were also investigated. The findings indicate that the deformation of PEN can be reduced and good sealing performance can be maintained with a suitable compressive loading of 0.6 MPa [7]. Pekson et al. developed a full-scale model for a SOFC short stack to reveal the thermomechanical response in the furnace environment, as well as to analyze the long-term creeping behavior of the critical components under the different environments [22]. Fang et al. analyzed the influence of the stack temperature distribution on the mechanical stress and the creep strain, as well as the mechanical failure probability of the key components. The effects of the electrode materials and the oxidation process on the stack mechanical properties were studied [23]. The effect of the volume fraction of Ni employed in the anodes on the thermal stress was investigated by Song et al. by considering the physical and chemical properties of various temperature-dependent materials [24].

In the SOFC stack, the interface between the adjacent cells requires a sealant as a gasket to prevent the leakage of gas [25,26]. Glass or glass-ceramic are rigid seals that can achieve a good sealing effect during the operation without a compress load [7,10]. In the recent years, the flexible mica-based compression seals have also attracted some attention. If the seal is not glued and not compatible, individual stack components can expand and contract freely during a thermal cycling. The interface between the cell and the interconnect is easily separated by the stress and creep of the material [27,28,29], which may cause the leakage of the gas, and further the overall performance degradation [30].

A lot of research work has been presented on the thermal stress of SOFC stacks, but most of these studies focus on the analysis of the thermal stress distribution characteristics. There are few studies on the key parameters of the geometric dimensions and comparative analysis under different flow modes. In order to save calculation time and improve the convergence of the model, many previous simulations adopted the method of simplifying the heat source. In this work, a multi-physics coupling model of the entire stack structure is established to solve the flow, electrochemical reactions, heat and mass transfer processes in the stack. The purpose of this paper is to present a numerical approach for evaluating the distribution of the thermal stress inside the SOFC stack, as well as to identify the locations of the stack that is prone to rupture and failure, which play a guiding role in the design of the stack. A finite element method for a three-cell SOFC stack is applied in this paper. Based on the thermo-electrochemical processes coupled in the developed model, the distribution of the temperature and its gradient, the produced current density, as well as the gas distribution within the stack are predicted and presented for various typical operating circumstances. The impacts of the gas flow channel design and the thickness of the electrolyte layer on the overall performance and the stress distributions are assessed and analyzed. The results of this work have a better understanding of the potential rupture and failure caused by the thermal stress inside the SOFC stack, and will also be helpful to the design of the stack.

2. Model Description and Governing Equations

2.1. Model Geometry and Assumptions

In this work, three flow channel configurations were studied (co-, counter- and cross-flow), as shown in Figure 1. PEN, a sealing component, and a bipolar plate make up the repeating unit cell for the stack. The PEN active size is 91 mm x 91 mm, with a thickness of 400 μm for the anode (Ni-YSZ), 70 μm for the cathode (LSM-YSZ), and 10 μm for the electrolyte layer (YSZ). The ribs have a width of 2 mm and a thickness of 1 mm, while the width of flow channels is 5 mm. Since the average flow velocity of fuel gas inside the stack is less than 1 m/s, and the average flow velocity of air is less than 3 m/s, the calculated Reynolds numbers are 2 and 36, respectively. Therefore, the flow is set to be laminar, and both fuel and air are considered as the ideal gas.

2.2. Governing Equations

2.2.1. Electrochemical Reaction Model

The reactions that occur inside the cell can be expressed as:

$$2{\mathrm{H}}_2+2{\mathrm{O}}^{2-}\to 2{\mathrm{H}}_2\mathrm{O}+4{\mathrm{e}}^{-}$$

(1)

$${\mathrm{O}}_2+4{\mathrm{e}}^{-}\to 2{\mathrm{O}}^{2-}$$

(2)

The oxidation of hydrogen occurs in the anode activation region, and oxygen is converted into oxide ions at the cathode which further transfers to the anode side through the electrolyte layer. The redox reactions are accompanied by charge transfer reaction, in which one mole of oxygen reacts with the transfer of 4 moles of electrons.

The operating voltage of the cell is calculated by open circuit voltage (E^OCV), activation loss (${\eta }_{act}$), concentration loss (${\eta }_{conc}$) and ohmic loss (${\eta }_{ohm}$):

$$E=E^{OCV}-\left({\eta }_{act}+{\eta }_{ohm}+{\eta }_{conc}\right)$$

(3)

where the open circuit voltage can be calculated as:

$$E^{OCV}=E_0+\frac{RT}{2F}\ln \left(\frac{c_{ref}^{{\mathrm{H}}_2}{\left(c_{ref}^{{\mathrm{O}}_2}\right)}^{0.5}}{c_{ref}^{{\mathrm{H}}_2\mathrm{O}}}\right)+\frac{RT}{4F}\ln p$$

(4)

E₀ is only a function of temperature, independent of the pressure and concentration of the reactants, obtained by the following formula:

$$E_0=1.253-2.4516\times {10}^{-4}T$$

(5)

The model employs the Butler–Volmer formula to express the relationship between the activation current and the activation polarization as well as the concentration polarization [31], as follows:

$$i^a=A_v^a{i}_0^a\left[\frac{c^{{\mathrm{H}}_2}}{c_{ref}^{{\mathrm{H}}_2}}\exp \left(\frac{{\alpha }^a{n}^{a}F\left({\eta }_{act}^a+{\eta }_{conc}^a\right)}{RT}\right)-\frac{c^{{\mathrm{H}}_2\mathrm{O}}}{c_{ref}^{{\mathrm{H}}_2\mathrm{O}}}\exp \left(\frac{\left(1-{\alpha }^a\right)n^{a}F\left({\eta }_{act}^a+{\eta }_{conc}^a\right)}{RT}\right)\right]$$

(6)

$$i^c=A_v^c{i}_0^c\left[\exp \left(\frac{{\alpha }^c{n}^{c}F\left({\eta }_{act}^c+{\eta }_{conc}^c\right)}{RT}\right)-\frac{c^{{\mathrm{O}}_2}}{c_{ref}^{{\mathrm{O}}_2}}\exp \left(\frac{\left(1-{\alpha }^c\right)n^{c}F\left({\eta }_{act}^c+{\eta }_{conc}^c\right)}{RT}\right)\right]$$

(7)

where the superscript a and c refer to the anode and cathode, respectively. A_v represents the effective reaction specific surface area of the active electrode (m²/m³) with relation to the TPB (Three-Phase Boundary) length. $\alpha$ denotes the transfer coefficients, and n represents the number of electrons transferred by the electrochemical reactions. The concentration of the gas component i is represented by cⁱ. i₀ characterizes the reactivity of the electrode under the effect of the equilibrium potential, which is one of the critical factors controlling the cell overall efficiency, the expressions are as follows [32]:

$$i_0=\frac{RT}{nF}{A}_i\times \exp \left(\frac{-E_a}{RT}\right)$$

(8)

in which A_i is the pre-referential factor for the electrodes, while E_a is the reaction activation energy given in

Table 1. Parameters applied for electrochemical reactions [33].

Parameter	Cathode	Anode
$A_i$	2.35 × 1011	6.54 × 1011
$E_a\left(\mathrm{kJ}/\mathrm{mol}\right)$	137	140
$\epsilon$	0.3	0.3
$\tau$	3	3
$A_v\left(1/\mathrm{m}\right)$	1.3 × 105	1.3 × 105

.

The transfer process of the electrons and ions inside SOFC can be described by the equation of t conservation of the electric potentials:

$$i_i=-{\sigma }_i^e\nabla {\mathrm{\Phi }}_i$$

(9)

$$i_e=-{\sigma }_e^e\nabla {\mathrm{\Phi }}_e$$

(10)

where ${\mathrm{\Phi }}_e$ and ${\mathrm{\Phi }}_i$ refer to the electric potentials of electrons and ions, respectively. ${\sigma }_i^e$ and ${\sigma }_e^e$ refer to the ionic conductivity and electronic conductivity, respectively. $i$ refers to the activation current generated when electrons and ions exchange charges. Among the SOFC functional materials, the main electron-conducting materials are LSM in the porous cathode, Ni in the porous anode, as well as the collector plate. The ion-conducting materials are YSZ in the electrode and electrolyte. Their conductivity is the temperature dependent, expressed as follows:

$${\sigma }_{electron,a}=\frac{4.2\times {10}^7}{T}\exp \left(\frac{-1200}{T}\right)$$

(11)

$${\sigma }_{electron,c}=\frac{9.5\times {10}^7}{T}\exp \left(\frac{-1150}{T}\right)$$

(12)

$${\sigma }_{ion}=3.34\times {10}^4\times \exp \left(\frac{-10300}{T}\right)$$

(13)

While in the porous electrodes, the effect of porosity on the conductivity also needs to be taken into account, resulting in an effective expression for the electron/ion conductivity:

$${\sigma }_{eff}=\left(\frac{1-\epsilon }{\tau }\right)\sigma \times V_{eff}$$

(14)

where $\tau$ is the tortuosity for the electron/ion conduction (i.e., the ratio of the actual path length inside the porous electrode to the straight-line distance), $\epsilon$ is the porosity. V_eff is the volume ratio of various materials in the porous structures, i.e., the volume percentage occupied by Ni, YSZ and LSM, respectively. The parameters employed in these electrochemical reactions are listed in Table 1.

2.2.2. Gas Flow

The Navier-Stokes equation is applied for the gas flow in the channels, while the Darcy term is added to represent the momentum conservation in the porous electrodes:

$$\nabla ·\left(\rho v\right)=S_{mass}$$

(15)

$$\rho v·\nabla v=-\nabla p+\nabla ·\left[\mu \left(\nabla v+v\nabla \right)-\frac{2}{3}\mu \nabla v\right]-\frac{\epsilon \mu v}{k}$$

(16)

where v is the velocity vector and k is the permeability. For the fuel and generated H₂O at the anode and the air at the cathode side, the average density ($\rho$) and dynamic viscosity ($\mu$) of the gas mixture can be calculated from the average of the individual gas component:

$$\rho =\frac{p{\displaystyle \sum x_i{M}_i}}{RT}$$

(17)

$$\mu ={\displaystyle \sum x_i{\mu }_i}$$

(18)

where x represents the molar fraction, M represents the molar mass, and the subscript i represents the substance component.

The mass source terms for two electrodes are formulated as:

$$S_{mass}^{an}=\frac{\left(M_{{\mathrm{H}}_2\mathrm{O}}-M_{{\mathrm{H}}_2}\right)i}{2F}$$

(19)

$$S_{mass}^{ca}=\frac{M_{{\mathrm{O}}_2}i}{4F}$$

(20)

2.2.3. Gas Species Transport

The electrochemical reactions take place at the electrode TPB regions. The gas must diffuse through the porous electrodes to reach the active reaction regions, in which the Knudsen diffusion is prominent. The governing equation for the mass conservation in the diffusion of gas species is given as:

$$\nabla ·j_i+\rho \left(\epsilon v·\nabla \right){\omega }_i=S_i$$

(21)

$$S_{{\mathrm{H}}_2}=M_{{\mathrm{H}}_2}\frac{i}{2F}$$

(22)

$$S_{{\mathrm{H}}_2\mathrm{O}}=-M_{{\mathrm{H}}_2\mathrm{O}}\frac{i}{2F}$$

(23)

$$S_{{\mathrm{O}}_2}=M_{{\mathrm{O}}_2}\frac{i}{4F}$$

(24)

$$j_i=-\rho D_i^{mk}\nabla {\omega }_i$$

(25)

$$D_i^{mk}=\frac{\epsilon }{\tau }{\left(\frac{1}{D_i^m}+\frac{1}{D_i^k}\right)}^{-1}$$

(26)

where S represents the mass source term, which represents the change in the mass of the gas species due to the electrochemical reactions, $\omega$ is the mass fraction, j is the mass flow rate, and D^mk is the effective diffusion coefficient determined from Fick’s diffusion coefficient (D^m) and Knudsen’s diffusion coefficient (D^k) [6].

2.2.4. Heat Transfer

The governing equation of the heat transfer is:

$$\rho C_{p}v·\nabla T=\nabla \left({\lambda }_{eff}\nabla T\right)+Q$$

(27)

where Q is the heat source term from the electrochemical reactions, C_p is the specific heat capacity, and ${\lambda }_{eff}$ is the effective thermal conductivity determined by:

$${\lambda }_{eff}=\left(1-\epsilon \right){\lambda }_s+\epsilon {\lambda }_g$$

(28)

where the subscripts s and g represent the solid materials and gases, respectively.

The heat source term mainly consists of the electrochemical reaction generated heat and the ohmic heat generated when the current passes through the solid materials in the porous electrodes, electrolytes and connectors:

$$Q=Q_{ohm}+Q_{act}$$

(29)

$$Q_{ohm}={\sigma }_i^e{\left(\nabla {\mathrm{\Phi }}_i\right)}^2+{\sigma }_e^e{\left(\nabla {\mathrm{\Phi }}_e\right)}^2$$

(30)

$$Q_{act}={\eta }_{act}{S}_{V}i$$

(31)

where S_V is the ratio between the electrochemical reaction active area to the volume.

Table 2. Properties for different SOFC components [34].

Parameter	Porosity	Permeability (m2)	Thermal Conductivity (W·m−1·K−1)	Thermal Capacity (J·kg−1·K−1)
Anode layer	0.3	2 × 10−11	11	450
Cathode layer	0.3	2 × 10−11	6	430
Electrolyte	-	-	2.7	550
Interconnect	-	-	44.5	475
Seal	-	-	0.064	1050

lists the parameters applied:

2.2.5. Thermal Stress-Strain Relation

All the solid materials used for the SOFC stack are assumed to be linearly elastic materials, and the deformation produced by the thermal stress is minimal. The thermal stress tensor is written as:

$$\sigma =D{\epsilon }_{el}+{\sigma }_0$$

(32)

where $\sigma$ is the stress tensor, ${\sigma }_0$ is the initial stress, ${\epsilon }_{el}$ is the elastic strain and D is the elastic matrix.

The elastic matrix of the isotropic materials is determined by elastic modulus (E) and Poisson’s ratio (v):

$$D=\frac{E}{\left(1+v\right)\left(1-2v\right)}\left[\begin{array}{cccccc}1-v& v& v& 0& 0& 0\\ v& 1-v& v& 0& 0& 0\\ v& v& 1-v& 0& 0& 0\\ 0& 0& 0& \frac{1-2v}{2}& 0& 0\\ 0& 0& 0& 0& \frac{1-2v}{2}& 0\\ 0& 0& 0& 0& 0& \frac{1-2v}{2}\end{array}\right]$$

(33)

The total strain consists of the elastic strain, thermal strain (${\epsilon }_{th}$), creep strain (${\epsilon }_{cr}$) and initial strain:

$$\epsilon ={\epsilon }_{el}+{\epsilon }_{th}+{\epsilon }_{cr}+{\epsilon }_0$$

(34)

where $\epsilon$ can be expressed by displacement vector (u):

$$\epsilon =\frac{\left(\nabla u+u\nabla \right)}{2}$$

(35)

while the thermal strain can be calculated by the following formula:

$${\epsilon }_{th}=\alpha \left(T-T_{ref}\right)$$

(36)

where $\alpha$ is the thermal expansion coefficient, and T_ref is the stress-free temperature.

After heated, the stack will experience a gradual relaxation of the thermal tension during the reduction phase.

Table 3. Mechanical parameters of SOFC components [35,36].

Parameter	Young’s Modulus (GPa)	CTE (10−6·K−1)	Poisson’s Ratio
Anode	58.1	13.2	0.39
Cathode	98	11	0.3
Electrolyte	148.6	10.3	0.32
Interconnect	60	15.5	0.3
Seal	0.019	13.9	0

lists the parameters used:

Where Young’s modulus represents the ratio of the stress in the unidirectional stress state to the strain in that direction, CTE is a physical quantity that measures the degree of thermal expansion of a solid material and Poisson’s ratio is the elastic constant reflecting the transverse deformation of the material.

2.3. Boundary Condition and Validation of the Model

The inlet boundary conditions are assumed to be the same for the stack with the different flow configurations, which is also true for the initial conditions. The fuel composition is 97% hydrogen and 3% water, with an inlet rate of 1500 sccm. While the gas flow rate of air is 4500 sccm. The SOFC stack is operated at 800 °C and the inlet temperature is also 800 °C. The top and bottom surfaces of the stack are set as adiabatic boundaries, the other surfaces as the natural convection and radiation heat transfer with the environment. The convective heat transfer coefficient is 2 W∙m⁻²∙K⁻¹ [37] and the emissivity coefficient 0.3 [38], respectively.

An important setting parameter for calculating thermal stress is the free stress temperature T_f. The recent investigations have shown that the reduction of the nickel oxide to nickel under a reducing atmosphere at 800 °C releases all strains [39]. Consequently, T_f = 800 °C is employed in the current work.

The commercial software COMSOL Multiphysics 5.5 was used for the simulation in this work including geometry, mesh and solve. To assess the reliability of the model, grid-independent tests were carried out. The maximum temperature of a single-layer stack is predicted for various grid numbers. It is found that the variance range of the highest temperature is less than 1% when the number of grids surpasses 0.5 × 10⁶, as shown in Figure 2b), i.e., the simulation results are independent of the number of the grids employed.

The predicted I-V curve is compared with the relevant experimental data at the same condition [34], as shown in Figure 2. The prediction is also presented in

Table 4. Comparison between simulation and experimental data [34].

Current Density (A·cm−2)	Experimental Data (Voltage/V)	Simulation Data (Voltage/V)	Relative Error	RMSE
0.0202539	1.04981	1.05822	0.81%	0.89%
0.0852187	1.00029	1.0143	1.41%
0.146767	0.965399	0.978482	1.36%
0.211737	0.937987	0.945463	0.80%
0.340538	0.890016	0.885344	0.53%
0.404939	0.864785	0.861982	0.33%
0.46649	0.841734	0.833635	0.97%
0.53431	0.816192	0.814011	0.27%
0.598711	0.792829	0.789091	0.48%
0.663681	0.768844	0.765418	0.45%
0.728082	0.745482	0.748285	0.38%
0.792483	0.721496	0.725234	0.52%
0.858594	0.6972	0.704364	1.03%
0.921854	0.673526	0.683182	1.44%
0.990814	0.647049	0.654836	1.21%

. It is found that the associated root mean square error (RMSE) is around 1%, indicating that the simulation and experimental results are in good agreement. The following equation is used to evaluate the RMSE:

$$RMSE=\sqrt{\frac{{\left(\frac{x_1-x_{e_1}}{x_{e_1}}\right)}^2+{\left(\frac{x_2-x_{e_2}}{x_{e_2}}\right)}^2+\dots +{\left(\frac{x_n-x_{e_n}}{x_{e_n}}\right)}^2}{n}}$$

(37)

where x_n is the simulated value, while x_e is the experimental one.

3. Results and Discussion

3.1. The Effect of Flow Configurations

Figure 3(a1–a3) shows the predicted temperature distribution of the SOFC stack under three different flow modes (i.e., the co-, counter-, and cross-flow configuration). In the co-flow mode, the averaged current density is 0.758 A/cm² and the maximum temperature is 893 °C, while the current density for the counter-flow case is averaged to be 0.826 A/cm² and the maximum temperature is 918 °C. This shows that the counter-flow mode is more effective for the electrochemical reactions, but it also generates more heat. In terms of the electricity generated, the cross-flow mode is least effective, with a predicted average current density of 0.711 A/cm² and the maximum temperature of 900 °C. In all the three cases, the maximum temperature is found to be concentrated in the cell central region of the stack.

The first principal stress is the maximum normal stress at a point over a specific discretized element, which can lead to the fracture within the materials and is often employed to predict the lifetime of SOFC stacks. Figure 3(b1–b3) shows the distribution of the first principal stress for the stack with various flow arrangements under the operational circumstances. As shown in Figure 3(b1–b3), the thermal stress distribution in the counter-flow mode is similar to that in the co-flow mode. It is observed that the first principal stress in the stack is concentrated on both sides of the stack, and the stress at the rib is higher than the surrounding stress. This is mainly due to the thermal expansion of the channel ribs squeezing the PEN, resulting in the stress concentration under the channel ribs. For the co-flow, counter-flow and cross-flow conditions, the highest thermal stress in the electrolyte layer is 111 MPa, 143 MPa, and 120 MPa, respectively.

The predicted thermal stress in the cross-flow stack is highest in the air outlet flow region, which is different from the prediction in the parallel flow stacks. This is because air flowing with a higher velocity is able to take away more heat to the environment through the side surfaces. The temperature is lowest at the four corners of PEN, and the corresponding first principal stress is also lower. Under the same conditions, the average current density of the co-flow stack is 6.6% higher than that of the cross-flow stack, while its maximum principal stress is lower. Under the three flow arrangements, the temperature and the distribution of the first principal stress between the cells are not much different because the top and bottom surfaces are in adiabatic conditions. The principal stress of the cell near the gas inlet is slightly higher than that of the other cells. This is because the low-temperature fuel and air inlets are close to the first cell, resulting in a relatively large temperature gradient.

The distribution of hydrogen and oxygen species at the electrolyte-electrode interface under the three flow modes is shown in Figure 4. It is found that the concentration distribution of hydrogen at the anode is similar to each other in the co-flow and counter-flow stack, i.e., its concentration decreases gradually along the main flow direction. In the outlet flow region of the channels on both sides of the cross-flow stack, the hydrogen concentration is lower, which affects the stack overall performance as described in the previous section.

It is predicted that, in the parallel-flow stacks, the oxygen distribution in the cathode is relatively uniform, and the O₂ concentration gradually decreases along the main flow direction. In the contrast, in the cross-flow stack, uneven distribution of the gases appears in the channels, and a low O₂ concentration occurs in some regions, which greatly reduces power generation performance. It is also found that the flow rate in the channels is smaller in the central area on both sides of the cross-flow stack. Optimization of the inlet manifold structure is essential in order to achieve a uniform air distribution.

In order to make the flow rate more uniform in the channels, the numbers and positions of the inlet and outlet are varied for the stack, as shown in Figure 5. Four modes are proposed, such as: (a) one inlet and one outlet, (b) one inlet and two outlets, (c) three inlets and three outlets, and (d) three inlets and two outlets, respectively. It should be noted that the total area is the same in all the above cases. Under 0.7 V operating condition, the averaged fuel flow velocity in each channel is predicted for various design cases, as shown in Figure 5. It is revealed that the standard deviation of the averaged flow velocities is about 0.1492 for the case (a), 0.0913 for (b), 0.0742 for (c) and 0.0465 for (d), respectively. In other words, the mode (d) is the one with the best uniformity (i.e., the smallest deviation) for the flow rate distribution among the four modes.

Figure 6 shows the temperature distributions predicted for the stack with four in-and outlet design cases. The maximum temperature in all design cases is similar to each other in terms of its distribution shape. However, the predicted average current density is quite different, as shown in

Table 5. Averaged current density and first principal stress predicted for various design cases.

Mode	Case 1	Case 2	Case 3	Case 4
Current density (A/cm2)	0.71	0.76	0.72	0.78
The first principal thermal stress (MPa)	120	113	113	112

. For the design (d), the averaged current density is about 10% higher than that for the design (a), while the maximum first principal stresses are even lower. Compared with the base case in Figure 5a, the maximum thermal stress is reduced by 6.7%. It indicates that the three-inlet and two-outlet design mode in the design (d) can generate a higher current density, while its internal temperature difference is lower with a small thermal stress, caused by better heat dissipation due to better flow uniformity.

3.2. The Effect of Electrolyte Layer Thickness

The thermal stress in the co-flow stack is predicted for five electrolyte layer thicknesses, as shown in

Table 6. Electrolyte layer thickness varied in this work (μm).

Case	Case 1	Case 2 (Baseline)	Case 3	Case 4	Case 5
Anode	405	400	395	390	385
Electrolyte	5	10	15	20	25

, while the thickness of the anode is also varied to keep the total thickness unchanged (i.e., 410 μm). The case 2 is assigned as the baseline one, as shown in Table 6.

Figure 7a compares the predicted average current density and the highest temperature for various electrolyte thicknesses. As the thickness decreases, the resistance to oxide ion transfer path decreases, but the fuel utilization rate reaches 86%, which correspondingly leads to a slight increase in both the average current density of the SOFC stack and the heat generated by the electrochemical reaction. In the case 1, the current density reaches to 0.7834 A/cm², while in case 5, the current density is about 0.6882 A/cm² and the maximum stress is 101 MPa when the thickness increases to 25 μm. However, when the electrolyte is too thin, the increase in the heat production leads to an increase in the thermal stress.

Figure 7b presents the maximum first principal stress predicted for a function of the electrolyte layer thickness. As the thickness of the electrolyte layer increases from 5 to 25 μm, the current density is reduced by 12%, and the first principal stress is also reduced by 6.5%. For the longer stacks with more cells involved, the temperature and its gradient may be even higher, and the decrease in the thickness of the electrolyte layer may cause an even bigger thermal stress, which is under investigation.

Figure 7c,d also presents the temperature and current density predicted for different cases with varied electrolyte layer thicknesses in the middle cell along the main flow direction. It is clear that the gas is continuously consumed and the generated current density gradually decreases along the main flow direction. In case 5, when the electrolyte thickness is the largest, it is found that the current density increases slightly at the end of the channel. This is because the hydrogen concentration in the top and bottom cell channels is higher than that in the middle cell, which causes the hydrogen concentration in the exhaust manifold to be higher than that at the end of the middle cell, as shown in Figure 7e. This may cause a backflow to slightly increase the hydrogen concentration at the end of the channel. In other cases, this phenomenon is not obvious. When the electrolyte layer thickness is set to 5 μm, the current density in the inlet region is larger than that predicted for the other design cases with the thick electrolyte layer, however it decreases more rapidly along the main flow direction. This is due to the strong electrochemical reactions, i.e., the strong oxygen consumption in the first half of the channel, followed by a formation of a lower oxygen concentration zone in the second half of the channels. As a result, the heat production in the second half of the channel is greatly reduced, causing the temperature to increase first and then decrease along the main flow direction, as shown in Figure 7d.

3.3. The Effect of Sealing Material

The stack seal provides excellent material compliance to accommodate large misalignment strains between the components due to mismatched coefficients of thermal expansion [40,41]. However, apart from the gasket-type seals, the rigid adhesive seal design does not require a load frame. Furthermore, the rigid adhesive seal will not leak easily after a long-term usage. However, these materials are less flexible compared to compression seal designs [25,42,43]. In this work, the thermal stress distribution of the stacks with two types of seals is studied, i.e., G18 (Glass-ceramic) [44] and Flexitallic 866 (gasket) [42], respectively. The two sealing materials are applied for different sealing methods, and the Young’s modulus and Poisson’s ratio are quite different, which will greatly affect the stress. The related parameters of the two sealing materials are summarized in

Table 7. Mechanical parameters of the sealing materials [42,44].

Seals	Young’s Modulus (GPa)	Poisson’s Ratio	CTE (10−6·K−1)
G18	14.4	0.28	11.1
Flexitallic 866	0.019	0	13.9

.

Figure 8a,b show the stress distribution of G18 sealing materials during the stack operating condition. It is found that the maximum first principal stress and the third principal stress are 16.9 and 12.8 MPa, respectively. It is also clear that the maximum tensile and compressive stresses are concentrated in the air inlet region.

Figure 8c,d show the stress distribution predicted for the Flexitallic 866 sealing materials. Different from that predicted for the glass sealing materials, Flexitallic 866 has a small form modulus as a mica-based material, so the predicted stress level is much lower than that of the glass materials and the color scale range of the predicted picture is different. It is clear that the tensile and compressive stresses are concentrated at the outer region. Under these materials, the maximum first principal stress of the electrolyte layer is 116 MPa and 111 MPa, respectively.

4. Conclusions

A multi-physics coupled numerical model is developed for a planar SOFC short stack aiming to determine the distribution of the gas species, current density, temperature and thermal stress. The effects of different flow configuration and the electrolyte layer thickness are investigated for the performance, temperature and its gradient, as well as thermal stress of the stack. According to the calculation results, the thermal stress is concentrated at the edge-side ribs of the electrolyte layer, when the co-flow mode is compared with the counter-flow mode, the co-flow mode shows the current density being 8% lower but the thermal stress 22% lower. The current density in the co-flow mode is 6.6% higher than that in the cross-flow mode, while the thermal stress is 7.5% lower. For the cross-flow stack, the thermal stress is concentrated in the gas outlet regions. An uneven flow distribution is found in the different channels, and various numbers and locations are investigated for the inlet and outlet designs to improve the flow distribution. The current density increases as the electrolyte thickness decreases, resulting in an increase in temperature and its gradient, as well as the thermal stress. In general, the thermal stress of case5 in the co-flow mode is the lowest, which is about 30% lower than that in the counter-flow mode, while the current density is reduced by 17%. Two sealing materials (glass/gasket) were also investigated to analyze the effects on the thermal stress distribution in the electrolyte layer. These results show that the selection of the sealing material is important because the gasket sealing material is beneficial to alleviate thermal stress concentration. It should be noted that pure hydrogen is used as the fuel, i.e., only exothermic electrochemical reactions occur inside the stack. If syngas is used as the fuel, there will be both endothermic and exothermic reactions inside, and the temperature gradient may be higher, which is worthy of our further study.