Effects of Methane Pre-Reforming Percentage and Flow Arrangement on the Distribution of Temperature and Thermal Stress in Solid Oxide Fuel Cell

: To obtain detailed information on the temperature ﬁeld and thermal stress ﬁeld inside the solid oxide fuel cell (SOFC) fueled with partially pre-reformed methane. A three-dimensional geometric and mathematical model of the SOFC is implemented by using the ﬁnite element method in the commercial software COMSOL Multiphysics ® . The coupling characteristics were analyzed for electrode chemical reaction, multi-component mass transfer, and heat transfer process under typical operating conditions, which was further applied for predicting and analyzing the thermal stress distribution. After model validation, parametric simulations are conducted to investigate how the methane pre-reforming percentage and ﬂow arrangement affect the temperature and the thermal stress of SOFC. The simulated results show that reducing the methane pre-reforming percentages can decrease the temperature maximum and the variation range of the ﬁrst principal stress, but will increase the possibility of carbon deposition. The maximum temperature of the counter-ﬂow is about 20 K lower than that of the co-ﬂow, and the ﬁrst principal stress variation range of the counter-ﬂow is 8.6 Mpa lower than that of the co-ﬂow. The methane pre-reforming percentages have a signiﬁcant effect on the heat transfer and the thermal stress, and the counter-ﬂow can improve the temperature uniformity and reduce the thermal stress variation range.


Introduction
To promote the rationalization of energy structure and achieve the goal of "carbon peaking and carbon neutral" as soon as possible, solid oxide fuel cell (SOFC), as a new type of energy conversion device, has attracted the attention of researchers with their advantages of green environment and high energy conversion efficiency [1]. The high temperature operating conditions (873 K-1273 K) coupled with the abundance of active sites in the porous nickel-based anode allow hydrocarbon fuels such as methane and natural gas to reform inside the SOFC to produce hydrogen, improving fuel compatibility. CH 4 has become one of the main fuels of SOFC due to its abundant sources (such as biogas, natural gas, etc.), and its low price [2]. However, SOFC operates under high-temperature conditions for a long time, and factors such as uneven temperature distribution inside the cell due to strong endothermic reforming reactions, carbon deposition from hydrocarbon fuel cracking, repeated oxidation-reduction in porous nickel-based anodes, and mismatch of thermal expansion coefficients of the constituent materials further aggravate the instability of cell life [3][4][5][6][7]. In addition, there are complex multi-phase, multi-scale and multi-field coupled transfer processes in the pores and three-phase boundary sites of micro/nanoporous electrodes where electrons, ions, and gases are combined in SOFC. It is difficult to fully explore the internal mechanisms of SOFC operating conditions and the laws of uneven temperature distribution affecting thermal stress distribution by practical experimental methods [1,2]. With the development of computers, numerical simulations have gradually become a research hot spot [8,9].
The direct internal reforming SOFC can directly supply the heat released by the charge transfer and the entropy change of the electrochemical reaction to the endothermic reforming reaction, which improves the energy utilization efficiency of the system [2,4]. However, direct internal reforming of methane in SOFC also faces some technical challenges. There is a local mismatch between the endothermic methane steam reforming reaction rate and the exothermic electrochemical rate [8]. Tseronis et al. [3] conducted several parametric studies on the direct internal reforming SOFC employing numerical simulation. It is shown that the cooling effect of the methane steam reforming reaction leads to a significant local temperature drop near the anode inlet and a large temperature difference inside the SOFC. At the same time, the key components of SOFC are temperature-dependent (thermodynamic parameters such as thermal expansion coefficient will change at different temperatures), which further exacerbates the instability and complexity of the stress field with temperature changes [6,10]. Jiang et al. [11] established the three-dimensional SOFC model through the commercial software COMSOL Multiphysics ® . It is found that more than 28-35% of the maximum first principal stress is caused by the temperature gradient; furthermore, about 47-54% of maximum first principal stress is caused by different CTEs (coefficients of thermal expansion). In addition, Xu et al. [12] and Clague et al. [13] directly used the sintering temperature as the stress-free reference temperature when solving the thermal stress under operating conditions, ignoring the residual stress of the cell generated by the completion of sintering. To further optimize the reaction rate distribution inside SOFC, it is necessary to explore the reforming reaction mechanism inside SOFC. Most researchers model global homogeneous methane-steam reforming and water-gas shift reaction within porous anodes by the described global reactions and corresponding kinetic reaction methods [2][3][4]9]. Schluckner et al. [14] validated the simulation results with experimental data from an industrial-scale SOFC single cell. The experimental and numerical results for gas-phase species and the calculated electrochemical performance were almost completely consistent with the chosen global as well as elementary reaction mechanisms. Ong et al. [15] established a one-dimensional membrane electrode assembly and found that before reaching the TPBs (TPBs are thin boundaries between the electrode and electrolyte micro grains surrounded by gaseous phases that fill the pores), not all CO participates in the water-gas shift reaction to convert to H 2 , and part of the CO also participates in the electrochemical reaction. Nickel in porous anodes is a good catalyst to catalyze the breakage of C-H bonds. Elemental carbon in the form of powder and whiskers covers the surface of the anode, blocking the pores of the anode and blocking the active part of the reaction [16]. Through numerical simulation, Klein et al. [17] and Vakouftsi et al. [18] studied the carbon deposition of tubular and planar SOFC. The results show that for different temperature levels and fuel compositions, carbon may form in the anode inlet region. Compared with the direct internal reforming, the external reforming reaction can avoid the problems of carbon deposition and local low temperature, but the system efficiency is reduced due to the need for a reactor and additional heating [2,19,20]. Some degree of external pre-reforming of methane may be a promising approach. The methane is partially reformed to produce H 2 , CO, and CO 2 before the fuel is transported to the SOFC anode [2,20]. CO 2 can be used as preventing agent of carbon deposition; however, some further degradation issues can occur if the atmosphere is not sufficiently reduced [21]. Partially pre-reforming of methane increases the ratio of hydrogen to methane, narrowing the gap between the reforming reaction rate and the electrochemical reaction rate, which mitigates the cooling effect in the region near the anode inlet [22].
Based on the previous research, to explore the effects of methane pre-reforming percentages and the flow arrangements on SOFC for the objective to achieve uniform distribution of temperature and thermal stress, a detailed multi-physics three-dimensional model of SOFC was developed based on the finite element method. The model couples the SOFC electrochemical reaction, chemical reaction, gas species diffusion, electron/ion transfer process, and heat and mass transfer process. The calculation of thermal stress is based on the electrochemical model and is obtained by adding an elastic constitutive equation coupled with the temperature field. It also considers the endothermic and exothermic sources of the reaction, the changes in thermodynamic properties of materials, and manufacturing residual stress.

Model Development
The following assumptions are adopted in this model: (1) The SOFC is running in a steady state.
(2) All gas species are deemed as ideal gases, and gas flow is assumed to be laminar flow.
(3) The electrochemical reaction area is uniformly distributed. The electronic and ionic conducting phases are continuous and homogeneous in the porous electrodes. (4) Electrochemical reactions only take place on the active layer. H 2 and CO have the same electrochemical reaction mechanism. (5) Ignoring the redox reaction of Ni-NiO.

Internal Reforming Reaction Mechanism
The SOFC anode channel is fueled with methane fuel, and the methane steam reforming (MSR) and the water-gas shift (WGSR) will occur in the porous nickel-based anode, MSR is shown in Equation (1), WGSR is shown in Equation (2).
The model developed by Haberman and Young is widely used to calculate the MSR reaction rate and the WGSR reaction rate [23], and this model is also used in this study. The control equation for reaction rates is shown in Table 1. Table 1. Control equation for reaction rates [23].

Reaction Rate Expression
Methane steam reforming reaction rate is the methane steam reforming reaction rate, R s is the water-gas shift reaction rate, k rf and k sf is the forward catalyzed reaction rate constants, K pr and K ps is the equilibrium constants.
The Boudard reaction and the methane cracking reaction are considered to be the main pathways for carbon formation at high temperatures, as they show the largest changes in Gibbs energy [16]. The Boudard reaction is shown in Equation (3), and the methane cracking reaction is shown in Equation (4).
Detailed chemical models validated for carbon deposition conditions are currently lacking, the bulk of the SOFC literature has adopted macroscopic models that rely on global kinetic steps rather than molecular-scale or detailed elementary-step models. Combined with the former models, thermodynamic predictions can provide useful qualitative guidance to evaluate the risk of carbon deposition [16,17,24]. This method is also adopted in this study, the thermodynamic tendency of carbon deposition can be expressed as [17]: where Q B and Q C represent the reactional quotient of Equations (3) and (4), K B and K C represent the equilibrium constant of Equations (3) and (4), a c is the carbon activity and is taken as 1 for the calculations. A value of γ less than 1 indicates carbon deposition at the anode, and it should be noted that the calculated γ only represents the thermodynamic tendency for carbon deposition.

Electrochemical Reaction Mechanism
In the anode, CO in the fuel also participates in the electrochemical reaction. Therefore, the electrochemical reaction is shown in Equation (7).
In the cathode, the electrochemical reaction is shown in Equation (8).
Considering the various polarization losses that occur under operating conditions, the actual SOFC operating voltage V cell can be expressed as: where E is the electromotive force; η act , η conc , and η ohm are the activation overpotential, concentration overpotential, and ohmic overpotential, respectively. The concentration overpotential η conc and ohmic overpotential η ohm can be calculated by Equations (10) and (11), respectively [3].
where y bulk and y TPB are the gas species mole fraction in the channel and electrode active layer, respectively; δ is the thickness of each component; σ is the conductivity of each component. The Butler-Volmer equation is used to describe the relationship between the activation overpotential and the activation current [25], which is expressed as: where i is the current density, i 0 is the exchange current density per active surface area (A/m 2 ), a a is the electron transfer coefficient and n is the number of electrons transferred by 1 mole of fuel oxidation.
Additionally, the exchange current density for the H 2 -H 2 O oxidation process can be formulated as: where P * H 2 is about 0.7 atm; P * O 2 = 4.9 × 10 8 exp − 2×10 5

RT
; i * H 2 and i * O 2 are empirical constants. It is reported that the rate-determining step of the CO-CO 2 oxidation process is quite similar to that of H 2 -H 2 O. Additionally, the rate of CO-CO 2 oxidation is about 2-3 times lower than that of H 2 -H 2 O oxidation under the same oxidant partial pressure. If the rate of CO-CO 2 oxidation process is assumed to be three times lower than that of H 2 -H 2 O oxidation process, the modeling result is in good agreement with experimental data. The CO-CO 2 oxidation electrochemistry can be described by Equations (14) and (15) [25].
i CO 0,a = i * CO P CO /P * The total exchange current density within the anode can be formulated as: The activation overpotential can be expressed as: where Φ elec and Φ ion are the potentials of electron and ionic, respectively, Φ eq is the equilibrium potential.

Transfer Process and Governing Equations
The transfer processes that take place inside SOFCs include the transfer of gas mass, momentum, energy, and charges (ions and electrons) [1][2][3]9].
The mass conservation describes the mass transfer process of the gas species in the gas channels and porous electrodes (support layer and active layer). The mass conservation can be described by Equation (18).
∇ρ e f f U = S m (18) where ρ eff is the average density of the gas mixture, U is the velocity vector, and S m is the total mass source term. Considering the interaction between gas molecules and the collision between gas molecules and the porous wall, the extended Fick model can be used to define the effective diffusion coefficient of component i in the mixed gas. The effective diffusion coefficient D e f f i is determined by Equation (19).
where the mixture averaged diffusivity D mix,i , Knudsen diffusivity D Kn,i are shown in Equation (20).
where r pore is the pore radius. The Binary diffusivity D i,j can be calculated by Equation (21) [3]. Momentum transfer refers to the process of transferring momentum from a highvelocity fluid layer to an adjacent low-velocity fluid layer in fluid flow. The momentum conservation can be described by Equation (22).
where the µ eff is the effective dynamic viscosity of the gas mixture and S d is the source term of the momentum equation.
The energy transfer includes the enthalpy change due to the diffusion of the gas mass, the heat release from the electrochemical reaction, the heat transfer between the gas and the solid medium, etc. The energy conservation can be described by Equation (23).
where m i H i is the enthalpy change due to mass diffusion of component i, and S T is the heat source generated by internal reforming, electrochemical reactions, and ohmic effects. Charge conservation includes electron transport in electrical conductors and oxygen ion transport in electrolytes. The electron and ion conservation can be described by Equations (24) and (25), respectively.
where σ elec,eff and σ ion,eff are the electrical conductivity and ionic conductivity, respectively.
i is the activation current density. The conservation equations of the transfer process of each physical physics and the corresponding source terms are shown in Table 2. Table 2. Governing equations and corresponding source terms of SOFC model [1][2][3]9].
The porous anode : Anode/cathode active layer (Electrochemical reaction): Interconnector : S T = σ elec (∇φ elec ) 2 Electrolyte and diffusion barrier layer : The porous anode : S g,i = (mR r + nR s )M i Anode/cathode active layer : S g,i = S ge,i = A Ve R e M i Electrons conservation: −∇ σ elec,e f f ∇φ elec = S elec Anode active layer : S elec = A Ve i act,a Cathode active layer : S elec = −A Ve i act,c Ions conservation: −∇ σ ion,e f f ∇φ ion = S ion Anode active layer : S ion = −A Ve i act,a Cathode active layer : S ion = A Ve i act,c English letters: ∆H is enthalpy; i is current density; M is molecular constant; R e is electrochemical reaction rate; R ohm is ohmic resistance; S is source term. Subscripts: act is activation; a/c is anode/cathode; elec is electron; eff is effective; g is gas phase; ion is ionic.

Thermal Stress Model
The thermal stress model assumed that all components of the SOFC undergo elastic deformation when subjected to thermal loads, and mechanical theory conforms to Hooke's law, the above assumptions are based on models in the literature [10][11][12]. The total strain Crystals 2022, 12, 953 7 of 17 consisted of elastic strain, thermal strain, and initial strain, which can be calculated by Equation (26).
Thermal strain can be calculated by Equation (27).
where α is the thermal expansion coefficient of the material, T ref is the free stress temperature, which is usually determined by the preparing and processing process of the material [12]. The stress-strain relationship of elastic material under thermal loading can be calculated by Equation (28).
where σ is the stress, D is the elasticity matrix, and σ 0 is the initial stress. First, simulate the residual stress from the end of sintering to cooling to room temperature (T ref = 1400 • C), and then simulate the thermal stress under operating conditions (T ref = 20 • C). The physical parameters of the dual-phase composite are calculated based on the Chin-Lung Hsieh calculation model [26]. Parameters for the thermal stress analysis of a SOFC, as shown in Table 3.

Geometric Model
To elucidate the reaction mechanism of SOFC fueled with partially pre-reformed methane and examine the local thermal stress distribution, a planar configuration composed of anode-supported SOFC and metallic interconnects was considered. The corresponding geometric model is based on an actual production planar SOFC, as shown in  Table 4. The properties for different components of the cell are shown in Table 5. Ce0.8Gd0.2O2), a cathode active layer (CAL, La0.6Sr0.4Co0.2Fe0.8O3-δ−Ce0.8Gd0.2O2), a cathode current collecting layer (CCCL, La0.6Sr0.4Co0.2Fe0.8O3-δ), as shown in Figure 1b. The flow arrangements for the co-flow and counter-flow are shown in Figure 1c,d. The selected simulation domain is the middle flow channel unit of the cell. The geometrical parameters of the simulation domain are shown in Table 4. The properties for different components of the cell are shown in Table 5.

Boundary Conditions
Different boundary conditions are needed for the different transport phenomena modeled and these can either be a fixed value or a fixed flux of the variable at the specific boundary, as shown in Table 6. The working pressure is 1 atm; the fuel utilization , where the n f,in and n f,out are the total molar fluxes of fuel gas species at the fuel inlet and outlet side.) of cases A, B, C, and D are 41.29%, 47.52%, 54.35%, and 55.29%, respectively. The methane pre-reforming means that before the fuel is transported into SOFC, methane is partially reformed to generate H 2 , CO, and CO 2 , and the molar ratio of the generated CO and CO 2 to CH 4 is called the pre-reforming percentage. Based on 30% pre-reformed [4], the initial molar fractions of fuel with 20%, 50%, and 70% pre-reformed are supplemented. The initial molar fractions of species with different pre-reforming percentages are shown in Table 7.

Solution Methodology
The SOFC model outlined above is implemented by using the finite element method in the commercial software COMSOL Multiphysics ® . The separation method was used to solve each physical field sequentially. The PARDISO and GMRES were used for solving the fully coupled nonlinear equations.

Model Validation
The anode-supported planar cell used in this study was manufactured by Ningbo SOFC-MAN Energy Technology Co., Ltd. The tested cells had 6 cm × 6 cm (actual working area 4 cm × 4 cm, ≥0.7 W/cm 2 ) and consisted of an anode (NiO-YSZ, AAL; NiO-YSZ, ASL) of 410 µm, an electrolyte(YSZ) of 10 µm, a diffusion barrier layer (GDC) of 3 µm, a cathode (LSCF-GDC, CAL; LSCF, CCCL) of 25 µm. Humidified hydrogen and air were introduced into the anode and cathode, respectively. The flow rate of hydrogen (mole farction: 95% H 2 , 5% H 2 O) was 400 sccm, and the airflow rate was 1000 sccm. The hydrocarbon fuel (mole farction: 10%H 2 , 67.5%H 2 O, 22.5%CH 4 ) and air were introduced into the anode and cathode, respectively. The flow rate of hydrocarbon fuel was 200 sccm, and the airflow rate was 800 sccm [27]. The operating temperature was maintained at 1073 K. Figure 2 shows the overall V-I curve predicted from the model, which is compared with the experimental data under the same operating conditions. To further validate the accuracy of the model by calculating the root mean square error (RMSE) between each experimental point and the simulated point in Figure 2. The calculated root mean square error of H2−SOFC and CH4−SOFC are 5.6% and 4.5%, respectively. RMSE can be calculated by Equation (29)  To further validate the accuracy of the model by calculating the root mean square error (RMSE) between each experimental point and the simulated point in Figure 2. The calculated root mean square error of H 2 -SOFC and CH 4 -SOFC are 5.6% and 4.5%, respectively. RMSE can be calculated by Equation (29).
4. Results and Discussion 4.1. The Effects of Methane pre-reforming Percentage Figure 3 shows the distribution of mole fraction along the middle section of the anode active layer in the forward y-axis under different methane pre-reforming percentages. It is found that the distribution of the mole fractions of CH 4 and H 2 are almost the same in different cases. They all decrease gradually along the forward direction of the y-axis. Figure 4 shows the distribution of MSR reaction rate under different methane pre-reforming percentages. It is found that the maximum value of MSR reaction rates of the four cases is located in the central region of the SOFC, and the MSR reaction rate difference between the highest and the lowest values of cases A, B, C, and D are 15.5, 32.2, 56.8, and 80.3, respectively. The predictions show that increasing the pre-reforming percentage decreases the uniformity of the MSR reaction rate distribution. This is because the initial molar fraction of hydrogen in case D is much larger than that in case A. More hydrogen is involved in the electrochemical reaction, and the heat released by the charge transfer and the entropy change of the electrochemical reaction supply more heat for the endothermic MSR and increase the maximum reaction rate. While at the anode outlet, the continuous MSR consumes a large amount of CH 4 , which results in the reaction rate decrease in this region. For comparison, the reaction rate of cases A and B were significantly higher than those of cases C and D near the anode inlet.      Figure 5 shows the distribution of γ along the middle section of the anode active layer in the forward y-axis under different methane pre-reforming percentages. It is worth noting that only the γ of case A near the anode inlet (y = 0-0.0002 m) is less than 1, indicating that case A produces carbon deposition in this region. The predictions show that increasing the pre-reforming percentage can inhibit coking. It has been reported that the electrochemical reaction of SOFC helps to suppress coking [4]. Lin et al. [28] explain that the deposited carbon is re-oxidized by oxygen ions passing through the electrolyte, and another possible explanation is that the steam backflow generated in the electrochemical reaction increases the local S/C ratio [4]. As stated before, the initial concentrations of H 2 and CO 2 increase significantly from case A to case D. The enriched H 2 content increases the current density and the steam generated in the electrochemical reaction increases the local S/C ratio, which is favorable for the inhibition of coking. On the other hand, the higher H 2 content hinders the methane cracking reaction, while the CO 2 enrichment favors carbon scavenging through the reverse Boudard reaction.  Figure 6 shows the distribution of temperature under different methane Pre-Reforming percentages. It is found that the temperature in the anode inlet is the lowest, and then rises along the flow direction (y-direction). The temperature reaches its maximum value in the regions close to the outlet of the anode channel. The reason is that the convective cooling of the gas in the channel and the endothermic reaction MSR in the early stage slow down the rise of the cell temperature, and then the temperature rise accelerated the charge exchange process and drives the electrochemical reaction to release more heat. At the same time, the cooling effect of the gas at the anode outlet drops to its weakest. The combined effect of these factors leads to the highest temperature at the anode outlet (y = 0.04 m). The maximum temperature keeps increasing when increasing the Pre-Reforming percentage. The maximum temperature of case D can reach 1180 K, which is nearly 100 K higher than that of case A. The predictions show that the Pre-Reforming percentage increases and more H2 are involved in the strongly exothermic electrochemical reaction. If the studied cell is located at the top or bottom of the stack, the effect of thermal radiation should be considered, and then the actual temperature will be much lower.  Figure 6 shows the distribution of temperature under different methane pre-reforming percentages. It is found that the temperature in the anode inlet is the lowest, and then rises along the flow direction (y-direction). The temperature reaches its maximum value in the regions close to the outlet of the anode channel. The reason is that the convective cooling of the gas in the channel and the endothermic reaction MSR in the early stage slow down the rise of the cell temperature, and then the temperature rise accelerated the charge exchange process and drives the electrochemical reaction to release more heat. At the same time, the cooling effect of the gas at the anode outlet drops to its weakest. The combined effect of these factors leads to the highest temperature at the anode outlet (y = 0.04 m). The maximum temperature keeps increasing when increasing the prereforming percentage. The maximum temperature of case D can reach 1180 K, which is nearly 100 K higher than that of case A. The predictions show that the pre-reforming percentage increases and more H 2 are involved in the strongly exothermic electrochemical reaction. If the studied cell is located at the top or bottom of the stack, the effect of thermal radiation should be considered, and then the actual temperature will be much lower. Figure 7 shows the distribution of the temperature gradient under different methane pre-reforming percentages. It is found that the maximum value of the temperature gradient of the four cases appears close to the anode inlet, and the maximum temperature gradient of case D can reach 59,883 K/m, which is approximately 6.7 times larger than that of case A. The predictions show that the temperature gradient keeps increasing when increasing the pre-reforming percentage. As stated before, the MSR reaction rate of case D is minimal near the anode inlet, resulting in an extreme mismatch between the endothermic reforming reaction rate and the exothermic electrochemical reaction rate in this region. It is also found that the temperature gradient of the electrode corresponding to the channel is significantly larger than that of the electrodes corresponding to the ribs on both sides. This is because the gas can diffuse directly from the channel to the corresponding porous electrode above the channel, while the gas diffusion to the electrodes corresponding to the ribs on both sides requires a longer diffusion path. As a result, the reaction at the electrode corresponding to the channel is more intense and the temperature gradient is larger. different methane Pre-Reforming percentages. Figure 6 shows the distribution of temperature under different methane Pre-Reforming percentages. It is found that the temperature in the anode inlet is the lowest, and then rises along the flow direction (y-direction). The temperature reaches its maximum value in the regions close to the outlet of the anode channel. The reason is that the convective cooling of the gas in the channel and the endothermic reaction MSR in the early stage slow down the rise of the cell temperature, and then the temperature rise accelerated the charge exchange process and drives the electrochemical reaction to release more heat. At the same time, the cooling effect of the gas at the anode outlet drops to its weakest. The combined effect of these factors leads to the highest temperature at the anode outlet (y = 0.04 m). The maximum temperature keeps increasing when increasing the Pre-Reforming percentage. The maximum temperature of case D can reach 1180 K, which is nearly 100 K higher than that of case A. The predictions show that the Pre-Reforming percentage increases and more H2 are involved in the strongly exothermic electrochemical reaction. If the studied cell is located at the top or bottom of the stack, the effect of thermal radiation should be considered, and then the actual temperature will be much lower.  Figure 7 shows the distribution of the temperature gradient under different methane Pre-Reforming percentages. It is found that the maximum value of the temperature gradient of the four cases appears close to the anode inlet, and the maximum temperature gradient of case D can reach 59,883 K/m, which is approximately 6.7 times larger than that of case A. The predictions show that the temperature gradient keeps increasing when increasing the Pre-Reforming percentage. As stated before, the MSR reaction rate of case D is minimal near the anode inlet, resulting in an extreme mismatch between the endothermic reforming reaction rate and the exothermic electrochemical reaction rate in this region. It is also found that the temperature gradient of the electrode corresponding to the channel is significantly larger than that of the electrodes corresponding to the ribs on both sides. This is because the gas can diffuse directly from the channel to the corresponding porous electrode above the channel, while the gas diffusion to the electrodes corresponding to the ribs on both sides requires a longer diffusion path. As a result, the reaction at the electrode corresponding to the channel is more intense and the temperature gradient is larger. The first principal stress is applied to the analysis of both ceramic and metallic systems and can distinguish between maximum tensile and compressive stresses (neglecting the shear stresses perpendicular to the plane direction). Therefore, this study focuses on comparing the first principal stresses, with positive and negative values representing tensile and compressive stresses, respectively [12]. Figure 8 shows the distribution of the first principal stress under different methane Pre-Reforming percentages. It is found that the distribution of the first principal stress is almost the same in different cases. It is also found that the electrode directly above the channel and the electrode near the wall are mainly subjected to tensile stress, while the electrode near the junction of the channel and the rib is mainly subjected to compressive stress. This is mainly because the upper and lower interconnects are constrained by fixed stresses and the channels cannot provide support, resulting in stress gradients in the x−direction of the electrode. The simultaneous presence of tensile and compressive stresses can easily induce the curling of material, which can lead to fracture. Figure 9 shows the variation range of the first principal stress under different methane Pre-Reforming percentages. It is found that the variation range of the first principal stress in case D (383.3 Mpa) is about 9.6% higher than in case A (349.7 Mpa). The increased range of stress fluctuations accelerates the possibility of cell crack generation and crack expansion. The first principal stress is applied to the analysis of both ceramic and metallic systems and can distinguish between maximum tensile and compressive stresses (neglecting the shear stresses perpendicular to the plane direction). Therefore, this study focuses on comparing the first principal stresses, with positive and negative values representing tensile and compressive stresses, respectively [12]. Figure 8 shows the distribution of the first principal stress under different methane pre-reforming percentages. It is found that the distribution of the first principal stress is almost the same in different cases. It is also found that the electrode directly above the channel and the electrode near the wall are mainly subjected to tensile stress, while the electrode near the junction of the channel and the rib is mainly subjected to compressive stress. This is mainly because the upper and lower interconnects are constrained by fixed stresses and the channels cannot provide support, resulting in stress gradients in the x-direction of the electrode. The simultaneous presence of tensile and compressive stresses can easily induce the curling of material, which can lead to fracture. Figure 9 shows the variation range of the first principal stress under different methane pre-reforming percentages. It is found that the variation range of the first principal stress in case D (383.3 Mpa) is about 9.6% higher than in case A (349.7 Mpa). The increased range of stress fluctuations accelerates the possibility of cell crack generation and crack expansion. resulting in stress gradients in the x−direction of the electrode. The simultaneous presence of tensile and compressive stresses can easily induce the curling of material, which can lead to fracture. Figure 9 shows the variation range of the first principal stress under different methane Pre-Reforming percentages. It is found that the variation range of the first principal stress in case D (383.3 Mpa) is about 9.6% higher than in case A (349.7 Mpa). The increased range of stress fluctuations accelerates the possibility of cell crack generation and crack expansion.

The Effects of Flow Arrangements
To further study the operation mechanism of SOFC fueled with partially pre-reformed methane, a model with different flow arrangements (co−and counter−flow) under the same operating conditions was simulated for case C (methane Pre-Reforming percentage 50%). The operating conditions are shown in Table 6. Figure 10 shows the distribution of MSR reaction rate for case C with different flow arrangements. It is found that near the anode outlet, the reaction rate of MSR in counterflow is significantly lower than that in co-flow, but the uniformity of reaction rate is better than that in co-flow. The main factors affecting the reaction rate of MSR are the amount of reaction gas and the heat required for the reaction. Due to the opposite direction of air and fuel flow in the counter−flow, the electrochemical reaction near the anode outlet is weak and does not provide sufficient heat for the endothermic MSR reaction. The methane conversion rate of the co-flow (52.9%) is higher than that of the counter-flow (41.7%), and the hydrogen conversion rate (59.9%) is lower than that of the counter-flow (69.9%). This is because, at the anode outlet, the high MSR reaction rate in the co-flow aggravates the consumption of methane and supplements the amount of hydrogen.

The Effects of Flow Arrangements
To further study the operation mechanism of SOFC fueled with partially pre-reformed methane, a model with different flow arrangements (co-and counter-flow) under the same operating conditions was simulated for case C (methane pre-reforming percentage 50%). The operating conditions are shown in Table 6. Figure 10 shows the distribution of MSR reaction rate for case C with different flow arrangements. It is found that near the anode outlet, the reaction rate of MSR in counterflow is significantly lower than that in co-flow, but the uniformity of reaction rate is better than that in co-flow. The main factors affecting the reaction rate of MSR are the amount of reaction gas and the heat required for the reaction. Due to the opposite direction of air and fuel flow in the counter-flow, the electrochemical reaction near the anode outlet is weak and does not provide sufficient heat for the endothermic MSR reaction. The methane conversion rate of the co-flow (52.9%) is higher than that of the counter-flow (41.7%), and the hydrogen conversion rate (59.9%) is lower than that of the counter-flow (69.9%). This is because, at the anode outlet, the high MSR reaction rate in the co-flow aggravates the consumption of methane and supplements the amount of hydrogen. Figure 11 shows the distribution of temperature of each component for case C with different flow arrangements. It is found that the maximum temperature of the co-flow is located in the anode outlet of the cell, and the maximum temperature is 1124 K. While for the counter-flow, the maximum temperature moves to the central region of the cell, and the maximum temperature is 1103 K. The predicted results reveal that the maximum temperature obtained in the counter-flow is about 20 K lower than that of co-flow, and the uniformity of temperature distribution in the counter-flow is better than that of co-flow. This is because the air and fuel flow in opposite directions in the counter-flow, resulting in a weaker electrochemical reaction between the fuel inlet and the air inlet. The electrochemical reaction is enhanced in the center of the cell, more heat is released, and the temperature reaches its maximum. In addition, it is also found that the temperature distribution trends of each component are almost the same.
arrangements. It is found that near the anode outlet, the reaction rate of MSR in counterflow is significantly lower than that in co-flow, but the uniformity of reaction rate is better than that in co-flow. The main factors affecting the reaction rate of MSR are the amount of reaction gas and the heat required for the reaction. Due to the opposite direction of air and fuel flow in the counter−flow, the electrochemical reaction near the anode outlet is weak and does not provide sufficient heat for the endothermic MSR reaction. The methane conversion rate of the co-flow (52.9%) is higher than that of the counter-flow (41.7%), and the hydrogen conversion rate (59.9%) is lower than that of the counter-flow (69.9%). This is because, at the anode outlet, the high MSR reaction rate in the co-flow aggravates the consumption of methane and supplements the amount of hydrogen. Figure 10. Distribution of MSR reaction rate for case C. Figure 11 shows the distribution of temperature of each component for case C with different flow arrangements. It is found that the maximum temperature of the co-flow is located in the anode outlet of the cell, and the maximum temperature is 1124 K. While for the counter−flow, the maximum temperature moves to the central region of the cell, and the maximum temperature is 1103 K. The predicted results reveal that the maximum temperature obtained in the counter-flow is about 20 K lower than that of co-flow, and the uniformity of temperature distribution in the counter-flow is better than that of co-flow.  This is because the air and fuel flow in opposite directions in the counter−flow, resulting in a weaker electrochemical reaction between the fuel inlet and the air inlet. The electrochemical reaction is enhanced in the center of the cell, more heat is released, and the temperature reaches its maximum. In addition, it is also found that the temperature distribution trends of each component are almost the same.  Figure 12 shows the distribution of the temperature gradient of each component for case C with different flow arrangements. It is found that the maximum temperature gradient of co-flow locates in the anode inlet of the cell, and the maximum temperature gradient is 35,100 K/m. While for the counter−flow, the maximum temperature moves to the anode outlet and the maximum temperature gradient is 40,100 K/m. The temperature gradient in the SOFC stack is mainly due to the local mismatch between the heat absorption of the endothermic reforming reaction and the heat release of the exothermic electrochemical reaction. As described in the reaction rate analysis above, the MSR reaction rate in the anode outlet is relatively small. In addition, it is also found that the temperature gradient distribution trends of each component are almost the same.  Figure 12 shows the distribution of the temperature gradient of each component for case C with different flow arrangements. It is found that the maximum temperature gradient of co-flow locates in the anode inlet of the cell, and the maximum temperature gradient is 35,100 K/m. While for the counter-flow, the maximum temperature moves to the anode outlet and the maximum temperature gradient is 40,100 K/m. The temperature gradient in the SOFC stack is mainly due to the local mismatch between the heat absorption of the endothermic reforming reaction and the heat release of the exothermic electrochemical reaction. As described in the reaction rate analysis above, the MSR reaction rate in the anode outlet is relatively small. In addition, it is also found that the temperature gradient distribution trends of each component are almost the same. Figure 13 shows the distribution of the first principal stress for case C with different flow arrangements, and Figure 14 shows the distribution of the maximum value of the first principal stress for each component. It is found that the distribution trend of the first principal stress is almost the same. The variation range of the first principal stress of the co-flow is 367 Mpa, and the variation range of the first principal stress of the counter-flow is 358.4 Mpa. The predicted results reveal that the variation range of the first principal stress of the counter-flow is about 8.6 Mpa lower than that of the co-flow. It also found that the tensile stress of the diffusion barrier layer is the largest and the tensile stress of the anode active layer is the smallest with two flow arrangements. From the diffusion barrier layer to the electrolyte layer, the stress value drops sharply from 350 Mpa to 33 Mpa. The first reason is that the thermal expansion coefficient of the electrolyte is the lowest among all components and Young's modulus is the largest, and it is more sensitive to external effects. The second reason is that the properties of the electrolyte and the cathode material are discontinuous. The third reason is that the diffusion barrier layer is thinner, the thermal expansion coefficient is higher, and the tensile stress on the cathode side is the largest. These factors lead to a sudden stress change in the electrolyte to the diffusion barrier layer. This area is the most vulnerable part of the PEN due to the irregular distribution of stresses.  Figure 12 shows the distribution of the temperature gradient of each component for case C with different flow arrangements. It is found that the maximum temperature gradient of co-flow locates in the anode inlet of the cell, and the maximum temperature gradient is 35,100 K/m. While for the counter−flow, the maximum temperature moves to the anode outlet and the maximum temperature gradient is 40,100 K/m. The temperature gradient in the SOFC stack is mainly due to the local mismatch between the heat absorption of the endothermic reforming reaction and the heat release of the exothermic electrochemical reaction. As described in the reaction rate analysis above, the MSR reaction rate in the anode outlet is relatively small. In addition, it is also found that the temperature gradient distribution trends of each component are almost the same.  principal stress is almost the same. The variation range of the first principal stress of the co-flow is 367 Mpa, and the variation range of the first principal stress of the counter-flow is 358.4 Mpa. The predicted results reveal that the variation range of the first principal stress of the counter-flow is about 8.6 Mpa lower than that of the co-flow. It also found that the tensile stress of the diffusion barrier layer is the largest and the tensile stress of the anode active layer is the smallest with two flow arrangements. From the diffusion barrier layer to the electrolyte layer, the stress value drops sharply from 350 Mpa to 33 Mpa. The first reason is that the thermal expansion coefficient of the electrolyte is the lowest among all components and Young's modulus is the largest, and it is more sensitive to external effects. The second reason is that the properties of the electrolyte and the cathode material are discontinuous. The third reason is that the diffusion barrier layer is thinner, the thermal expansion coefficient is higher, and the tensile stress on the cathode side is the largest. These factors lead to a sudden stress change in the electrolyte to the diffusion barrier layer. This area is the most vulnerable part of the PEN due to the irregular distribution of stresses.

Conclusions
A detailed multi−physics three-dimensional model is developed in this study. The effects of the methane Pre-Reforming percentages and the flow arrangements on the temperature and the thermal stress were investigated. The conclusion is as follows: The methane Pre-Reforming percentages have a significant effect on the heat transfer and the thermal stress. Reducing the methane Pre-Reforming percentages can improve 348  principal stress is almost the same. The variation range of the first principal stress of the co-flow is 367 Mpa, and the variation range of the first principal stress of the counter-flow is 358.4 Mpa. The predicted results reveal that the variation range of the first principal stress of the counter-flow is about 8.6 Mpa lower than that of the co-flow. It also found that the tensile stress of the diffusion barrier layer is the largest and the tensile stress of the anode active layer is the smallest with two flow arrangements. From the diffusion barrier layer to the electrolyte layer, the stress value drops sharply from 350 Mpa to 33 Mpa. The first reason is that the thermal expansion coefficient of the electrolyte is the lowest among all components and Young's modulus is the largest, and it is more sensitive to external effects. The second reason is that the properties of the electrolyte and the cathode material are discontinuous. The third reason is that the diffusion barrier layer is thinner, the thermal expansion coefficient is higher, and the tensile stress on the cathode side is the largest. These factors lead to a sudden stress change in the electrolyte to the diffusion barrier layer. This area is the most vulnerable part of the PEN due to the irregular distribution of stresses.

Conclusions
A detailed multi−physics three-dimensional model is developed in this study. The effects of the methane Pre-Reforming percentages and the flow arrangements on the temperature and the thermal stress were investigated. The conclusion is as follows: The methane Pre-Reforming percentages have a significant effect on the heat transfer and the thermal stress. Reducing the methane Pre-Reforming percentages can improve 348

Conclusions
A detailed multi-physics three-dimensional model is developed in this study. The effects of the methane pre-reforming percentages and the flow arrangements on the temperature and the thermal stress were investigated. The conclusion is as follows: The methane pre-reforming percentages have a significant effect on the heat transfer and the thermal stress. Reducing the methane pre-reforming percentages can improve the MSR reaction rate uniformity and lower the temperature maximum and the variation range of the first principal stress, but it will increase the possibility of carbon deposition.
The flow arrangement can adjust the temperature field distribution trend and the thermal stress variation range. The maximum temperature of the counter-flow is closer to the center of the cell, and the maximum value is about 20 K lower than that of the co-flow, but the maximum temperature gradient is higher than that of the co-flow. The first principal stress variation range of the counter-flow is 8.6 Mpa lower than that of the co-flow.