Thermal Stress in Full-Size Solid Oxide Fuel Cell Stacks by Multi-Physics Modeling

: Mechanical failures in the operating stacks of solid oxide fuel cells (SOFCs) are frequently related to thermal stresses generated by a temperature gradient and its variation. In this study, a computational fluid dynamics (CFD) model is developed and further applied in full-size SOFC stacks, which are fully coupled and implemented for analysis of heat flow electrochemical phenomena, aiming to predict thermal stress distribution. The primary object of the present investigation is to explore features and characteristics of the thermal stress influenced by electrochemical reactions and various transport processes within the stacks. It is revealed that the volume ratio of the higher thermal stress region differs nearly 30% for different stack flow configurations; the highest probability of potential failure appears in the cell cathodes; the more cells applied in the stack, the greater the difference in the predicted temperature/thermal stress between the cells; the counter-flow stack performs the best in terms of output power, but the predicted thermal stress is also higher; the cross-flow stack exhibits the lowest thermal stress and a lower output power; and although the temperature and thermal stress distributions are similar, the differences between the unit cells are bigger in the longer stacks than those predicted for shorter stacks. The findings from this study may provide a useful guide for assessing the thermal behavior and impact on SOFC performance.


Introduction
As a high-temperature (600-800 • C) power generation device, a solid oxide fuel cell (SOFC) can directly convert chemical energy stored in hydrocarbon fuel gases into electrical energy with high conversion efficiency and the minimal emission of pollutants [1].However, its short working life has limited its commercial applications.According to reports, coking, poisoning, particle coarsening/agglomeration, and stress generated in the fuel electrodes are the primary problems leading to SOFC performance degradation.The stresses involved in SOFCs consist of the following: the residual stress/strain caused by the manufacturing process (including redox strain), thermal stress/strain, high-temperature creep strain during long-term operation, and chemical expansion coupled with electrochemicalmechanical effects [2].The thermal stresses due to high temperatures and mismatching between the coefficients of thermal expansion (CTEs) of the electrode materials during the stable operation of SOFC stacks are significant contributors to the cracks in the electrodes [3].These cracks can result in component breakdown and thus affect the durability of the stack.Therefore, a reasonable thermal stress distribution is essential to enhance stability and prolong the operating life of the stacks [4].Currently, most numerical studies of the thermal stresses adopt the method of one-way coupling, wherein the temperature field obtained from multiple physical fields is used to simulate the thermal stress field.
Given that current experimental studies take a longer time and have a complex equipment arrangement which is often difficult to implement, numerical simulation methods have become important for analyzing the thermal stresses in SOFC stacks [5].
A few of the published work have been conducted for various aspects of SOFCs in general, from micro-to macroscopic scales using various numerical analysis methods, as outlined in [6].Some of these studies focus on SOFC stack operating processes and overall performance including the thermal stress distribution but mainly take the loopedmaterial components or even assume the positive/electrolyte/negative (PEN) as a grouped cell component [7].Others have looked at single channels, ignoring the effects of constraints on other SOFC components (such as seals, metal frames, etc.) [8].It is also a fact that the thermal stress distribution in SOFC stacks depends not only on the operational conditions of the cells, encompassing the gas properties, but also on the overall stack designs/configurations [9].It is therefore essential to consider all details of the cell components included in the stacks (particularly industrial-sized ones) in a systematic manner (the so-called full-size stack model and simulation method).
Lin et al. [7] first considered the complete stack model by coupling computational fluid dynamics (CFD) or the finite element method (FEM).They investigated the effects of stack support conditions, the viscosity of glass ceramic sealants, the temperature gradients, and the mismatched coefficients of thermal expansion (CTEs) between the components on seal performance.Nakajo et al. [10] employed a similar approach to evaluate the stress and contact pressure distribution of various components in the single-cell stack model.Lin [11] developed a three-dimensional model of the SOFC stack using a commercial finite element analysis (FEA) tool to study the effect of external loads on the thermal stress distribution of planar designed SOFC stacks with compression seal designs.Zheng et al. [12] developed an SOFC stack model using a multi-physics coupling method to evaluate the effects of different flow patterns, manifold configurations, electrolyte layer thicknesses, and sealant materials on the thermal stress distribution of planar SOFC stacks.Xu et al. [13] established a two-dimensional mathematical model to study the thermal response of tubular methanolfueled SOFCs.The results indicate that the harmful temperature gradients caused by high current density cannot be resolved by excessive air supply.The key to effectively control the current density is to optimize the operating potential.
In simulating thermal stresses in large-scale SOFC stacks, the simplification of multiphysics fields or geometric structures is mainly employed.As for simplifying multi-physics fields, Peksen et al. [14] studied a 36-cell stack to determine the distribution of thermal stress, deformation degree, and transient stress response without considering the heat generated by the actual electrochemical reactions and electron transfer process within the stack components.Yuan et al. [15] developed a Fortran language code to study the impact of various air intake configurations on the performance as well as the thermal stress of a 20-cell stack in cross-flow mode, and it was revealed that changing the air flow direction significantly reduced the thermal strain while leaving the power output of the SOFC stacks unaffected.However, the coupling effect of the electrochemical reactions and transport process with the temperature field was simplified, which does not consider the structural details of the stack manifolds and seals, etc. Wang et al. [16] studied the influence of different fuel cycling modes on the maximum temperature of the SOFC system, and the application of heat exchangers can effectively improve the system's electrical efficiency.
In terms of simplifying geometric structures, the homogenization method has achieved certain research outcomes.For different research objectives, the PEN components, channels, and ribs, or even the entire unit cell, can be assumed to be a single entity.Navasa et al. [17] established a 3D SOFC stack model consisting of 100 unit cells based on the homogenization method and the volume averaging technique, significantly reducing computational time while ensuring computational accuracy.Meanwhile, a detailed distribution of internal thermal stresses can be retrieved by sub-models.Using the same approach, Rizvandi [18] investigated the effects of different operating parameters and structures on the thermal stress distribution in the stack.Miao et al. [19] conducted the first study of the localized Energies 2024, 17, 2025 3 of 25 fracture phenomena at the full-scale stack level by employing fracture mechanics to examine the localized mechanical failures.This model can be utilized for optimizing stack design and operating conditions, aiming for reducing the risk of localized mechanical ruptures and enhancing stack reliability.
In summary, the analysis of the thermal stresses generated during the operation of SOFCs has been mainly conducted through numerical simulation.However, there are several key issues that remain unresolved: (1) Only gas flow within the stack is considered, aiming to optimize the stack structure and to improve gas distribution uniformity.On the other hand, only the average temperature field is considered in some cases, aiming to explore the distribution of temperature or stress distribution within the stack, without accounting for the electrochemical reactions.(2) Only the electrochemical reactions are considered but simplified by reducing the three-dimensional electrochemical reactions occurring within the electrodes into twodimensional ones.
(3) The full-size of the stack and the coupling of the multi-physics processes are considered, yet the stack size is relatively small and cannot fully reflect the characteristics appearing inside the industrial-sized stacks.(4) Large-sized stacks and multi-physics coupling are all considered, while the real geometric details may be omitted.
These simplifications significantly reduce the computational time and improve the convergence performance of the developed models, but they may introduce some uncertainties in the simulation results and even incorrect conclusions.In reality, the structures of large-sized SOFC stacks are complex, and the distribution of the physical parameters is very uneven.Therefore, it is necessary to establish an SOFC model that can reflect the actual geometric structure and details of the stack by coupling with multi-physics transport processes, aiming to identify the main contributions of the operating/configuration conditions, as well as their effects on the thermal stress distribution.
Based on the discussion above, it is essential to develop the suitable models and simulation methods for evaluating the internal thermal stress and its distribution and identifying stack regions and locations with too big thermal stresses, as well as for investigating the thermal stress difference between the cell layers.For this investigation, a mathematical framework is formulated and applied for capturing the thermal stresses in a full-size SOFC stack assembly consisting of 10 cells including all components and structural details.
The developed model uses the finite element method (FEM) to address all coupled governing equations for mass, momentum, charge, and heat transport coupled with the electrochemical reactions that occur in each single cells.The computational domain includes all components such as the PEN components, the inlet/outlet manifolds, the seals, the frames, and the current collectors, etc.To verify the precision of the model, the predicted results of the model will be compared with the experimental data.
Creative modeling techniques are employed to implement large-scale fully coupled multi-physics modeling without simplifying geometric details and physical processes, aiming to assess the impact on the thermal stress.The predicted results are the identified regions of the highest thermal stress generated during the steady-state operation of hydrogen-fueled SOFC stacks.Furthermore, the thermal stress differences between the unit cells in the stack are analyzed and compared with the predicted results for a shorter stack to understand the influence of stack size.Predicted failure probabilities under the aforementioned operating conditions are provided, aiming to offer new insights on identifying the risk levels of potential component failures and improving the design of SOFC stacks.

Modeling Methods and Validation
This section details the creation of a linked multi-physics process model operating under steady-state conditions.The model comprises the geometric model, governing equa-tions, and boundary conditions.Additionally, the subsequent assumptions are employed in this study: (1) The gases in the model flow laminarly (based on the estimated Reynolds number) are ideal and incompressible.(2) The thermal radiation inside the SOFC is neglected.
(3) Homogeneity is exhibited by the porous electrode.(4) The constant temperature of 800 • C is maintained by the external ambient.
(5) The anode, electrolyte layer, and cathode materials are isotropic and linear elastic materials.(6) The thermophysical characteristics of all materials composing the SOFC stack remain constant, irrespective of the local temperature.(7) The anode, cathode, and electrolyte layer interfaces, along with the connecting body and sealing material, form a continuous structure that permits collective deformation without fracturing.(8) The external loads imposed during the preparation and assembly of the stack and the residual stresses in the material are negligible.

Stack Parameters and Geometric Model
The investigated structure of the SOFC stack is shown in Figure 1.The air flows into/out from the stack through the inlet/outlet manifolds, while the fuel gas is isolated from the air by the seals.Each unit cell includes a metal bipolar plate, anode, electrolyte layer, cathode, and seal.

Electrochemical Model
The controlling equation of the ion/electron transport process in the porous electrode is as follows: where i, σ, and ϕ represent the resulted current density, electron/ionic conductivity, and potential.The subscripts i and e represent the ions and electrons, while ▽ represents the gradient operators.The stack's operating voltage is given by From a vertical perspective, the first layer of the cells at the bottom is designated as the base unit cell, with the number of the unit cells gradually increased upward, and the final cell is designated as the top one.For the co-/counter-flow configuration of the stack, based on the symmetry characteristics of the stack physical field as well as its structure, only half of the stack geometry is included with the middle plane designated as the symmetry; while for the cross-flow arrangement of the stack, the complete stack geometry is included.The cell dimensions are 91 × 91 mm 2 , and the number of flow channels in the stack is 13, with a width ratio of 5:2 between the flow channels and the ribs [20], as shown in Table 1.The controlling equation of the ion/electron transport process in the porous electrode is as follows: ) where i, σ, and ϕ represent the resulted current density, electron/ionic conductivity, and potential.The subscripts i and e represent the ions and electrons, while ▽ represents the gradient operators.The stack's operating voltage is given by where E is the operating voltage, and E OCV is the open-circuit voltage.Inside the cell, η act , η ohm , and η conc denote the activation polarization, ohmic polarization, and concentration polarization.
The open-circuit voltage is typically determined using the standard electrode potential, gas constant R, temperature T, pressure p, the Faraday constant, and the concentration of substance i [21]: where E 0 is solely dependent on the temperature and remains unaffected by changes in pressure and gas concentration.This value can be computed using the subsequent equation: In this paper, the electrode reaction kinetics is elucidated through the Butler-Volmer equation.This equation delineates the connection between activation polarization and local current density, as represented by the expressions for the anode and cathode provided below [22]: where α a and α c denote the electron transfer coefficients at the anode and the cathode, the parameter n represents the quantity of electrons transferred by 1 mole of the gas reactants engaged in the electrochemical reactions, A v represents the effective specific surface area active within the porous electrodes, the superscripts a and c represent the anode and cathode, and the concentration of the gas component i is denoted by c i .i 0 is the exchange current density, as expressed by the following equation [23]: The pre-exponential factor for the electrodes, A i , and the activation energy of the reaction, E a , are represented by the following: Based on Ohm's Law, the ohmic polarization of a SOFC stack can be described as where h is the thickness, e is the electrolyte layer, and int is the internal connector.
The SOFC anode and cathode activation polarization can be found by the subsequent equation: where x is the molar fraction.Temperature-dependent electrode conductivity is expressed as follows: The electronic conductivity of the interconnect material is 8.7 × 10 5 S•m −1 .
Adjusting the conductivity of electrons/ions is crucial given the influence of the porous microstructure of the electrodes.The specific correction formula used is as follows: where θ represents the porosity, τ denotes the tortuosity, and V eff signifies the percentage of each material per unit volume.Table 2 lists the parameters involved in the concerned electrochemical reactions [24].Based on the assumption of laminar gas flow, the continuity equation in the gas flow channels is as stated below [25]: ∇•(ρv) = S mass (17) where ρ denotes the density of the mixed gases, and v denotes the velocity vector.The momentum conservation in the porous electrodes is depicted by adjusting the typically employed Navier-Stokes equation with the inclusion of a Darcy coefficient that considers the porosity of the electrodes: Amongst these, k signifies the permeability, whereas µ signifies the dynamic viscosity of the mixed gases.
The density and dynamic viscosity of the mixed gases can be acquired by the following: where M i denotes the molar mass of the substance.S mass stands for the mass source term composed of the consumption and production of the gas components during the electrochemical reactions.The mass source terms inside the porous electrode of the anode and cathode are formulated as follows:

Gas Species Transport Equation
The electrochemical reactions within the SOFC stack take place at the triple-phase boundary (TPB) at the interface between the porous electrode and the electrolyte layer.The conservation equation for every gas species can be formulated as follows [26]: where j i is the mass flow rate, computed by The effective diffusion coefficient of the gas species i is represented by D mk i [26].

Heat Transfer Equation
The universal energy conservation equation is implemented as follows [12]: which describes the balance of energy in a system, accounting for various forms of energy and their interconversion within the system.C p is the specific heat capacity of the mixed gases, λ eff denotes the effective thermal conductivity, and Q denotes the heat source term inside the SOFC stack.What is determined is the thermal conductivity effectiveness of the porous electrodes [12]: where λ s and λ g denote the thermal conductivity of the solid and gas.
The coefficients employed in the governing equations discussed above are detailed in Table 3 [27].Anticipated to be minimal is the resulting thermal stress-induced distortion, presuming all solid materials are linearly elastic.The stress can be articulated as follows [12]: where σ represents the stress tensor, D stands for the elastic matrix, ε el denotes the elastic strain, and σ 0 is the initial stress.The total strain comprises two components as follows: where ε th is the thermal strain.
The elastic strain can be acquired from [28]: where ∇u is the gradient of the displacement field u, representing the spatial variation in displacements within the material.(∇u) T is the transpose of the gradient of the displace- ment field.The thermal strain is derived from Energies 2024, 17, 2025 where α is CTE, and ∇T denotes the temperature change.
where T ref represents the state of the zero thermal stress adopted at 800 • C [29].
Commonly specified for the elastic matrix of an isotropic material is the following [12]: where Y is Young's modulus, and v is the Poisson ratio.
The PEN structure of SOFC stacks is composed of brittle ceramic materials, exhibiting very apparent brittle characteristics.The typical cermet electrodes with potential defects may appear; hence Weibull statistics are commonly used to evaluate SOFC failure probabilities.It is worth noting that this method does not consider the variation in the defects over time.In the Weibull method, the failure probability of the brittle materials under multi-axial stress is determined as follows [30]: where V is the volume of the brittle material, σ w is the Weibull strength of the brittle material, m is the Weibull modulus, and V 0 is the reference volume obtained by the following formula [31]: where σ max represents the maximum principal stress of the brittle materials.
The detailed parameters needed for the current study are outlined in Table 4 [32].

Effective Properties of Porous Composite Materials
According to Chin-Lung-Hsieh's model [33], the maximum and minimum limits of the material properties of Ni and YSZ after mixing are determined by both the material properties and respective volume fractions, while the average of the maximum and minimum limits assumed to be the mechanical properties is applied to determine Young's modulus, the Poisson ratio, and the CTE of the electrode materials under different mixing ratios.
The equivalent properties of the porous composite electrodes are related to not only the constituent materials but also the porosity, which is a key factor affecting their properties.It is revealed that the influence of the porosity on the CTE can be ignored [34], and the relationship of Young's modulus and the Poisson ratio with the porosity can be expressed as follows [35]: When the porosity is 0.3, the estimated mechanical parameters of the electrode materials and other components are shown in Table 5 [36].

Boundary Conditions
This section provides the boundary conditions and physical parameters.The cathode, electrolyte layer, and anode are made of LSM, YSZ, and Ni/YSZ materials, respectively.LSM is renowned for its high electronic conductivity and excellent catalytic properties, which are crucial for the oxygen reduction reaction on the cathode; YSZ is known for its high oxygen ion conductivity at high temperatures, which is essential for the operation of SOFCs; and Ni/YSZ exhibits good conductivity and high porosity, allowing for the effective diffusion and reaction of fuel gases, making it a suitable material for SOFCs.
The connecting plate and sealing material are made of stainless steel and Flexitallic 866, respectively.Stainless steel is suitable for use as the interconnect material in SOFCs due to its excellent high-temperature resistance, corrosion resistance, and mechanical strength.Flexitallic 866 is suitable for use as the sealing material in SOFCs due to its excellent sealing performance and high-temperature resistance.
As shown in Table 6, the fuel consists of 97% hydrogen and 3% water, with an intake flow rate of 1500 sccm.The gas inflow temperature is 800 • C, and the air intake flow rate is 4500 sccm.The initial SOFC stack temperature after preheating is also set at 800 • C. The stack's upper and lower surfaces are assumed to be adiabatic, while the remaining surfaces experience natural convection and radiative heat transfer, with a convective heat transfer coefficient of 2 W m −2 K −1 [20] and surface emissivity of 0.3 [37], exchanging heat with the surroundings.For the electron and ion potential field, the boundary conditions are grounded (0 V) at the bottom of the stack and N × 0.7 V at the top, where N represents the quantity of cell layers.

Grid Arrangement and Model Validation
The SOFC stack consists of multiple repeating SOFC unit cells, so the grid detail is exampled by one of the unit cells.In this study, the PEN is divided into three layers, as shown in Figure 2. As shown in Table 7, a coarser mesh is used for the interconnectors; the metal frame and the seals are solely engaged in heat transfer, so the meshing is also coarser, while the PEN, the ribs, and the fluid channels involve the multi-physics parameters, e.g., heat and mass transfer with a large temperature gradient, and a dense meshing is required.
dicted results.The root mean square error (i.e., RMSE) is around 1%, as defined in Equation (40).As depicted in Figure 3b, the modeling outcomes closely match the empirical data.The grid independence is assessed in the x, y, and z directions of each cell PEN in terms of the highest temperature at the top and bottom cells and at the midpoint of the centerline.When the number of the grids is doubled (i.e., increased to 6.01 × 10 5 ), the resulting variation in the highest temperatures of the three points is below 1%, as shown in Figure 3a, suggesting that this grid number (6.01 × 10 5 ) is adequate for achieving a solution independent of the grid.
Energies 2024, 17, x FOR PEER REVIEW 11 of 26 metal frame and the seals are solely engaged in heat transfer, so the meshing is also coarser, while the PEN, the ribs, and the fluid channels involve the multi-physics parameters, e.g., heat and mass transfer with a large temperature gradient, and a dense meshing is required.The grid independence is assessed in the x, y, and z directions of each cell PEN in terms of the highest temperature at the top and bottom cells and at the midpoint of the centerline.When the number of the grids is doubled (i.e., increased to 6.01 × 10 5 ), the resulting variation in the highest temperatures of the three points is below 1%, as shown in Figure 3a, suggesting that this grid number (6.01 × 10 5 ) is adequate for achieving a solution independent of the grid.
Under the same conditions, the experimental data [27] are compared with the predicted results.The root mean square error (i.e., RMSE) is around 1%, as defined in Equation (40).As depicted in Figure 3b, the modeling outcomes closely match the empirical data.Under the same conditions, the experimental data [27] are compared with the predicted results.The root mean square error (i.e., RMSE) is around 1%, as defined in Equation (40).As depicted in Figure 3b, the modeling outcomes closely match the empirical data.

Results and Discussion
The distribution of the gas flow, temperature, and thermal stress is predicted and analyzed for the studied 10-cell stack under different gas flow configurations (i.e., the co-/counter-/cross-flow stack design), while the difference in the temperature and thermal stress between the cell layers is outlined and discussed.The stack power generation, temperature distribution, and stress distribution of the 3-cell and 10-cell stacks arranged in co-flow mode are compared, aiming to identify the stack size effect.Finally, the failure probabilities of the electrodes and the stack are calculated and presented for all the above operating conditions.

Effect of Flow Arrangement and Configurations
The prediction of the stack parameters is intricately linked to the electrochemical reaction rates and temperature, subsequently influencing the distribution of thermal stress.Understanding these relationships is crucial for optimizing stack performance and durability.
Figure 4 presents the gas distribution within the stack across various flow arrangements.The observations indicate that the gas concentration distribution in each channel of the stack is highly consistent for both H 2 and O 2 under the co-and counter-flow cases, while it becomes uneven due to the flow pattern particularity in the cross-flow case.It is also true that the consumption of hydrogen is greater in the counter-flow stack.In the co-flow stack, as a result of electrochemical reactions, hydrogen diminishes in the direction of gas flow, experiencing a sharp decline in the initial portion of the H 2 /O 2 flow channels, followed by a more gradual decline in the latter portion of the channels.
For the counter-flow stack, the hydrogen distribution matches that of the co-flow case, but due to oxygen flowing in the opposite direction, the rate of its flow is three times greater than that of hydrogen.Oxygen experiences rapid depletion in the region corresponding to the hydrogen inlet area.In the cross-flow mode, the oxygen in the air is rapidly consumed near the hydrogen inlet region due to electrochemical reactions, hydrogen gradually decreases in the direction of flow, but its distribution is relatively uneven, while the oxygen concentration is higher due to the depletion of hydrogen in the channels on the opposite side near the hydrogen outlet region.
The distribution of the gases between different cells or within different flow channels in the same cell is commonly assessed by the normalized mass flow, which is calculated as shown below [38]: The mass flow rate at cell i is denoted as m L,i .Firstly, the gas distribution uniformity in the flow channels is investigated by selecting cells 1, 5, and 10 in the 10-cell stack (representing the bottom, the middle, and the top cell) in the co-flow mode, as shown in Figure 5.It was found that the gas distribution in the flow channels between different cells was similar with a minor difference.Based on this finding, a comparison of the hydrogen distribution within the flow channels in the three flow configurated stacks can be carried out by arbitrarily selecting a specific cell.In this paper, the first cell (i.e., the bottom one) is taken as an example, and it is found that the hydrogen distribution within each channel of the initial cell is analogous for the co-and counter-flow modes, owing to the comparable stack configuration, while a more uniform hydrogen distribution is predicted within the channels for the cross-flow mode due to the longer distance between the air inlet and the hydrogen inlet manifolds.In this case, less influence on each other can be expected, while the longer transfer distance makes the flow rate of hydrogen in both side channels (i.e., channel 1 and channel 13) smaller, which may further affect the stack performance as highlighted in the following sections.Firstly, the gas distribution uniformity in the flow channels is investigated by selecting cells 1, 5, and 10 in the 10-cell stack (representing the bottom, the middle, and the top cell) in the co-flow mode, as shown in Figure 5.It was found that the gas distribution in the flow channels between different cells was similar with a minor difference.Based on this finding, a comparison of the hydrogen distribution within the flow channels in the three flow configurated stacks can be carried out by arbitrarily selecting a specific cell.In this paper, the first cell (i.e., the bottom one) is taken as an example, and it is found that the hydrogen distribution within each channel of the initial cell is analogous for the coand counter-flow modes, owing to the comparable stack configuration, while a more uniform hydrogen distribution is predicted within the channels for the cross-flow mode due to the longer distance between the air inlet and the hydrogen inlet manifolds.In this case, less influence on each other can be expected, while the longer transfer distance makes the The generated current density for the three flow stacks is 0.639 Acm −2 , 0.656 Acm −2 , and 0.560 Acm −2 , respectively.This indicates that for the same operating conditions, the counter-current flow stack exhibits superior power generation capabilities, which aligns with the forecasted hydrogen utilization, while the cross-flow mode has the worst power generation performance.
The cloud plot of the hydrogen concentration distribution for all three flow stacks shows that the hydrogen concentration is lower in the flow channels close to the outlet region.If one wishes to use the assumptions (e.g., the average current density or uniform gas flow) to calculate the stack gas distribution and temperature field distribution, a certain pre-assessment procedure is required and should be conducted because the distribution of current density within the stack is not consistently uniform.Overall, only the coupled multiphysics field modeling may reflect or predict the realistic conditions of the multi-physics parameters distributed within the long stacks, such as in this study.
Energies 2024, 17, x FOR PEER REVIEW 14 of 26 flow rate of hydrogen in both side channels (i.e., channel 1 and channel 13) smaller, which may further affect the stack performance as highlighted in the following sections.The generated current density for the three flow stacks is 0.639 Acm −2 , 0.656 Acm −2 , and 0.560 Acm −2 , respectively.This indicates that for the same operating conditions, the counter-current flow stack exhibits superior power generation capabilities, which aligns with the forecasted hydrogen utilization, while the cross-flow mode has the worst power generation performance.
The cloud plot of the hydrogen concentration distribution for all three flow stacks shows that the hydrogen concentration is lower in the flow channels close to the outlet region.If one wishes to use the assumptions (e.g., the average current density or uniform gas flow) to calculate the stack gas distribution and temperature field distribution, a certain pre-assessment procedure is required and should be conducted because the distribution of current density within the stack is not consistently uniform.Overall, only the coupled multi-physics field modeling may reflect or predict the realistic conditions of the multi-physics parameters distributed within the long stacks, such as in this study.
Figure 6 illustrates the temperature distribution of various flow stack configurations.There are conspicuous differences in temperature distribution under different flow conditions.Temperatures are focused at the central region of the stack and decrease progressively towards the periphery.In the co-flow stack configuration, the maximum temperature reaches 895 °C, with the high-temperature region situated at the center of the stack.This is because the electrochemical reactions in the stack mainly occur against the inlet region, and this heat is further carried to the center of the stack through the flowing of the generated products.The highest temperature of the counter-flow stack is 905 °C, which is the highest within the three flow modes studied.The predicted high-temperature region primarily resides within the hydrogen inlet area, attributable to the air flow rate being three times greater than that of hydrogen, and the specific heat capacity of the products on the air side is also greater than that on the hydrogen side, so most of the heat is carried forward to the hydrogen inlet, forming a high-temperature zone; the highest temperature of the cross-flow mode stack is 896 °C, and the high-temperature zone is primarily situated in the central region near the air outlet, which is also caused by the gas flow effect.The flow direction of the gases, the electrochemical exotherm, and the thermal coupling between the fluid and the solid ultimately determine the temperature field distribution of the SOFC stack.There are conspicuous differences in temperature distribution under different flow conditions.Temperatures are focused at the central region of the stack and decrease progressively towards the periphery.In the co-flow stack configuration, the maximum temperature reaches 895 • C, with the high-temperature region situated at the center of the stack.This is because the electrochemical reactions in the stack mainly occur against the inlet region, and this heat is further carried to the center of the stack through the flowing of the generated products.The highest temperature of the counter-flow stack is 905 • C, which is the highest within the three flow modes studied.The predicted high-temperature region primarily resides within the hydrogen inlet area, attributable to the air flow rate being three times greater than that of hydrogen, and the specific heat capacity of the products on the air side is also greater than that on the hydrogen side, so most of the heat is carried forward to the hydrogen inlet, forming a high-temperature zone; the highest temperature of the cross-flow mode stack is 896 • C, and the high-temperature zone is primarily situated in the central region near the air outlet, which is also caused by the gas flow effect.The flow direction of the gases, the electrochemical exotherm, and the thermal coupling between the fluid and the solid ultimately determine the temperature field distribution of the SOFC stack.Due to the difference in the temperature field predicted for various flow arrangements, the thermal stress distributed in the stacks also differs considerably.The distribution of the thermal stresses at the electrodes for the three flow cases (i.e., co-/counter-/cross-flow stack configurations) is shown in Figure 7(a1-a3).The predicted findings, regardless of the flow arrangements, reveal that the electrolyte layer experiences a significantly higher thermal stress than the anode and cathode, while the thermal stress distribution trend is basically similar in all three layers.In the co-flow mode, the thermal stresses are mainly concentrated at the edges in the width direction, with the inlet and outlet regions being slightly higher than the intermediate regions, while in the counterflow mode, the stresses are not only concentrated at the edges in the width direction but also in the air outlet (i.e., the hydrogen inlet) region.On the other hand, in the cross-flow mode, the thermal stresses are primarily focused in the hydrogen inlet area, with the thermal stresses below the ribs of the flow channels being greater than those below the gas channels, which are similar to those found in the co-and counter-flow modes.The primary reasons for this are the thermal expansion of the ribs and the pre-compression of the cell functional layers of the stack, resulting in a concentration of the mechanical stresses in this region.The maximum thermal stresses in the electrolyte layer are 138 MPa, 142 MPa, and 146 MPa, respectively.In the anode, they are 32 MPa, 33 MPa, and 34 MPa, respectively.In the cathode, they are 22 MPa, 24 MPa, and 25 MPa, respectively.So the maximum thermal stress in all three flow configurations occurs in the electrolyte layer.The temperature distribution cloud is shown in Figure 6.A comprehensive analysis shows that the highest thermal stress distribution in the electrodes does not coincide exactly with that of the highest temperature.For instance, the highest temperature is predicted in the center region of the stack close to the hydrogen inlet, while the maximum thermal stress is mainly concentrated close to the edge of the electrode width close to the Due to the difference in the temperature field predicted for various flow arrangements, the thermal stress distributed in the stacks also differs considerably.The distribution of the thermal stresses at the electrodes for the three flow cases (i.e., co-/counter-/cross-flow stack configurations) is shown in Figure 7(a1-a3).The predicted findings, regardless of the flow arrangements, reveal that the electrolyte layer experiences a significantly higher thermal stress than the anode and cathode, while the thermal stress distribution trend is basically similar in all three layers.In the co-flow mode, the thermal stresses are mainly concentrated at the edges in the width direction, with the inlet and outlet regions being slightly higher than the intermediate regions, while in the counter-flow mode, the stresses are not only concentrated at the edges in the width direction but also in the air outlet (i.e., the hydrogen inlet) region.On the other hand, in the cross-flow mode, the thermal stresses are primarily focused in the hydrogen inlet area, with the thermal stresses below the ribs of the flow channels being greater than those below the gas channels, which are similar to those found in the co-and counter-flow modes.The primary reasons for this are the thermal expansion of the ribs and the pre-compression of the cell functional layers of the stack, resulting in a concentration of the mechanical stresses in this region.The maximum thermal stresses in the electrolyte layer are 138 MPa, 142 MPa, and 146 MPa, respectively.In the anode, they are 32 MPa, 33 MPa, and 34 MPa, respectively.In the cathode, they are 22 MPa, 24 MPa, and 25 MPa, respectively.So the maximum thermal stress in all three flow configurations occurs in the electrolyte layer.
Energies 2024, 17, x FOR PEER REVIEW 16 of 26 discussed in the next section (i.e., the probability of failure).It is found that the distribution of the predicted HTSR is completely different for the three flow cases.Although the maximum thermal stresses are similar for the three flow stacks, the counter-flow mode is significantly subject to a bigger region with high thermal stress than both the co-flow and cross-flow modes, while in the cross-flow stack, the thermal stress is maximum in the air outlet region, owing to the increased flow rate of the air, which may carry more heat out through the air flowing.
The total volume of the electrolyte layers in the stack is 10 × 91 × 91 × 0.02 mm 3 .It is found that the volume ratio with the thermal stress on the electrolyte layers exceeding 100 MPa in the three flow stacks is 19.7%, 37.0%, and 7.1%, respectively, with the biggest ratio predicted for the counter-flow stack.In other words, the biggest and smallest ratio is about five times differed in the counter-and cross-flow modes.Clearly, as far as the flow arrangement is concerned, the thermal stress distribution is predominantly affected by the stack air flow configuration.The SOFC stack in this study consists of 10 identical repeating unit cells, but the difference in the gas distribution and the electrochemical reaction rates between the cells will lead to different temperature and thermal stress distributions.The analysis of these differences between the cell layers will help to identify the weak stack components or regions with big thermal stress generated, which could potentially result in the component failures of the operating stacks.
The thermal stress on the electrolyte layers in the SOFC stack is much greater than The temperature distribution cloud is shown in Figure 6.A comprehensive analysis shows that the highest thermal stress distribution in the electrodes does not coincide exactly with that of the highest temperature.For instance, the highest temperature is predicted in the center region of the stack close to the hydrogen inlet, while the maximum thermal stress is mainly concentrated close to the edge of the electrode width close to the hydrogen inlet and air outlet regions.In other words, the formation of the thermal stresses is not entirely determined by the highest temperature.Further analysis should be carried out to correlate the highest thermal stress with the highest temperature and the temperature gradient under various stacks.
The concept of high stress regions is introduced in this study, as discussed below, aiming to provide both qualitative and quantitative analyses of the thermal stress-affecting factors.Figure 7(a2-c2) presents the regions where the thermal stress in the electrolyte layer exceeds 100 MPa (hereafter, the high thermal stress region or HTSR).The reason for choosing 100 MPa as the basis for analysis is mainly based on two points: first, as shown in Table 8, if 110 MPa is chosen as the analysis basis, the volume proportions of the thermal stress regions within the electrolyte layer in the three flow patterns are 3.5%, 11.1%, and 1.5%, respectively; and if 120 MPa is chosen as the analysis basis, then the proportions are 0.3%, 2.4%, and 0.2%, respectively.These proportions are essentially consistent with the characteristics obtained for the high-stress region when selecting 100 MPa as the analysis basis.Second, considering the data for failure probability analysis, which will be discussed in the next section (i.e., the probability of failure).It is found that the distribution of the predicted HTSR is completely different for the three flow cases.Although the maximum thermal stresses are similar for the three flow stacks, the counter-flow mode is significantly subject to a bigger region with high thermal stress than both the co-flow and cross-flow modes, while in the cross-flow stack, the thermal stress is maximum in the air outlet region, owing to the increased flow rate of the air, which may carry more heat out through the air flowing.The total volume of the electrolyte layers in the stack is 10 × 91 × 91 × 0.02 mm 3 .It is found that the volume ratio with the thermal stress on the electrolyte layers exceeding 100 MPa in the three flow stacks is 19.7%, 37.0%, and 7.1%, respectively, with the biggest ratio predicted for the counter-flow stack.In other words, the biggest and smallest ratio is about five times differed in the counter-and cross-flow modes.Clearly, as far as the flow arrangement is concerned, the thermal stress distribution is predominantly affected by the stack air flow configuration.
The SOFC stack in this study consists of 10 identical repeating unit cells, but the difference in the gas distribution and the electrochemical reaction rates between the cells will lead to different temperature and thermal stress distributions.The analysis of these differences between the cell layers will help to identify the weak stack components or regions with big thermal stress generated, which could potentially result in the component failures of the operating stacks.
The thermal stress on the electrolyte layers in the SOFC stack is much greater than that on the anodes and cathodes, so in this section, the electrolyte layer is taken as the main object for further evaluation.It should be noted that the first cell is specified near the inlet side and the 10th layer specified away from the inlet side.
In Figure 8a, it is found that the difference in the distribution of hydrogen uniformity across the cell layers of the stack is small for all three flow forms, with the flow rate of hydrogen monotonically decreasing from cell layer 1 to 10. Figure 8b shows the difference between the maximum and minimum temperatures predicted for the electrolyte layers in the cells.It is evident that the maximum temperature for the co-, counter-and cross-flow stacks is 895 • C, 905 • C and 896 • C, respectively, which is observed in cell 7; while the lowest temperatures is 847 • C, 844 • C and 835 • C, respectively, which is predicted in cell 1.The predicted maximum and minimum temperatures look to be similar in all cells, but the maximum temperature in the counter-flow stack is notably greater compared to that in the co-and cross-flow stacks, while the minimum temperature reaches its peak in the co-flow stack, then in the counter-flow stack, and is lowest in the cross-flow stack.

Failure Probability Analysis
The Weibull distribution is widely used in various fields, including electrical engineering, mechanical engineering, etc., for predicting the probability distribution of the potential faults, failures, and lifespans [39].This indicates the broad applicability of the Weibull distribution, which can be used for different types of equipment and systems, including SOFCs.There are various parameter estimation methods available, which are particularly important for SOFCs operated for a long period with only limited failure data [40].Moreover, even with a small sample size, the Weibull distribution can provide reliable parameter estimations and reliability assessments.However, there are also limitations, such as the uncertainty of the parameter estimation and a high dependence on data quality and initial settings [41].Therefore, when applying the Weibull distribution for assessing SOFC failure probabilities, it is necessary to carefully select suitable parameter estimation methods and ensure the quality and integrity of the data.
The Weibull analysis is employed for predicting the probability of failure of the electrodes within the stack, as plotted in Figure 9a-d.The electrolyte layer in the counter-flow configuration is investigated in this section, followed by the volume ratio of the electrolyte layer at different thermal stresses.The corresponding contribution to the failure probability is defined as the likelihood of failure in the region above the specified thermal stress as a proportion of the total failure likelihood, as shown in Figure 9a.It is found that the contribution of the failure probability increases rapidly at the beginning and then grows slowly with an increase in the proportion of high thermal stress areas, which is different from the linear growth pattern of volume proportion of the high thermal stress areas.When the thermal stress is above 100 MPa, the volume share of the electrolyte layer is Because the bottom cell has a high mass flow rate, the heat emitted from the electric stack is rapidly dissipated, while some of the heat carried away by the fluid is transferred to the upper cell layers, resulting in the highest temperature in the top cell.The gases in the cross-flow mode flow in the opposite direction, and the high-temperature gas in the exhaust manifold can transfer some of the heat to the gas in the inlet region to warm it up and speed up the electrochemical reaction.The electrochemical reaction rate in the co-and cross-flow modes is lower than that in the counter-flow, where the heat production is also lower.As depicted in Figure 6, the high-temperature region in the co-flow mode is at the center of the stack and radiates evenly around the stack, with heat gradually transferred upwards with the gas flow, resulting in a gradual rise in the minimum temperature; the high-temperature region in the cross-flow mode is biased towards the hydrogen inlet, while the minimum temperature is lower in the air inlet region due to the low hydrogen content and the slower electrochemical reaction.
Using 100 MPa as the cut-off value of the predicted thermal stress, the volume ratio of each electrolyte layer subjected to thermal stress beyond this value is calculated and compared with that in different stacks, as shown in Figure 8c.In all three flow cases, the volume ratio is larger in the 1st cell electrolyte layer and becomes gradually smaller and then increases until the maximum value is reached in the 10th cell layer.The volume ratios are 44.9%,48.5%, and 19.5%, respectively.The counter-flow stack has the largest volume ratio, while the cross-flow stack has the smallest.Not only is the counter-flow stack subjected to the higher thermal stress, but the volume where the thermal stress in each electrolyte layer exceeds 100 MPa is also significantly larger than the other two modes.The cross-flow stack is subjected to the maximum thermal stress more than the co-and counterflow modes, but the volume ratio where the stress in each cell layer exceeds 100 MPa is the smallest, indicating that most of the electrolyte layers are subjected to low thermal stress.

Failure Probability Analysis
The Weibull distribution is widely used in various fields, including electrical engineering, mechanical engineering, etc., for predicting the probability distribution of the potential faults, failures, and lifespans [39].This indicates the broad applicability of the Weibull distribution, which can be used for different types of equipment and systems, including SOFCs.There are various parameter estimation methods available, which are particularly important for SOFCs operated for a long period with only limited failure data [40].Moreover, even with a small sample size, the Weibull distribution can provide reliable parameter estimations and reliability assessments.However, there are also limitations, such as the uncertainty of the parameter estimation and a high dependence on data quality and initial settings [41].Therefore, when applying the Weibull distribution for assessing SOFC failure probabilities, it is necessary to carefully select suitable parameter estimation methods and ensure the quality and integrity of the data.
The Weibull analysis is employed for predicting the probability of failure of the electrodes within the stack, as plotted in Figure 9a-d.The electrolyte layer in the counterflow configuration is investigated in this section, followed by the volume ratio of the electrolyte layer at different thermal stresses.The corresponding contribution to the failure probability is defined as the likelihood of failure in the region above the specified thermal stress as a proportion of the total failure likelihood, as shown in Figure 9a.It is found that the contribution of the failure probability increases rapidly at the beginning and then grows slowly with an increase in the proportion of high thermal stress areas, which is different from the linear growth pattern of volume proportion of the high thermal stress areas.When the thermal stress is above 100 MPa, the volume share of the electrolyte layer is relatively small, accounting for 37%, but it contributes significantly to the failure probability, accounting for 78%.When the pressure is below 100 MPa, the volume share of the electrolyte layer is relatively large, accounting for 63%, but it contributes minimally to the failure probability, accounting for 22%.The probability of the component failure increases slowly.Therefore, 100 MPa is considered as the basis value of the thermal stress, as mentioned previously.
Energies 2024, 17, x FOR PEER REVIEW 19 of 26 relatively small, accounting for 37%, but it contributes significantly to the failure probability, accounting for 78%.When the pressure is below 100 MPa, the volume share of the electrolyte layer is relatively large, accounting for 63%, but it contributes minimally to the failure probability, accounting for 22%.The probability of the component failure increases slowly.Therefore, 100 MPa is considered as the basis value of the thermal stress, as mentioned previously.
From the cell component perspective, the failure probability is greatest in the porous cathode, followed by the electrolyte layer and the porous anode, and the Weibull strength and modulus presented in Table 4 can account for this.A higher Weibull strength means a more stable material and vice versa.The dense electrolyte layers have the biggest Weibull strength and the second biggest Weibull modulus, which means that the electrolyte layer is extremely stable.While the porous cathode has the lowest Weibull strength and Weibull modulus, it is most probable to result in cracks.This is due to the fact that the electrodes in the counter-flow stack mode are subjected to a higher thermal stress than those in the other flow stacks, not only in terms of the maximum thermal stress but also in terms of the volume ratio with higher thermal stress.The cross-flow stack has the lowest probability of failure, and although the thermal stress is the highest, the volume ratio of the higher thermal stress is the smallest, which is consistent with the previous findings.

Comparative Study with Shorter and Longer Stacks
The co-flow mode is a conventional and commonly employed flow mode in SOFC stacks.For the distribution of thermal stress within a shorter stack, a comparative investigation is conducted under identical configurations and operating conditions, aiming to identify any unique features only appearing in the longer stacks as well as the stack size From the cell component perspective, the failure probability is greatest in the porous cathode, followed by the electrolyte layer and the porous anode, and the Weibull strength and modulus presented in Table 4 can account for this.A higher Weibull strength means a more stable material and vice versa.The dense electrolyte layers have the biggest Weibull strength and the second biggest Weibull modulus, which means that the electrolyte layer is extremely stable.While the porous cathode has the lowest Weibull strength and Weibull modulus, it is most probable to result in cracks.This is due to the fact that the electrodes in the counter-flow stack mode are subjected to a higher thermal stress than those in the other flow stacks, not only in terms of the maximum thermal stress but also in terms of the volume ratio with higher thermal stress.The cross-flow stack has the lowest probability of failure, and although the thermal stress is the highest, the volume ratio of the higher thermal stress is the smallest, which is consistent with the previous findings.

Comparative Study with Shorter and Longer Stacks
The co-flow mode is a conventional and commonly employed flow mode in SOFC stacks.For the distribution of thermal stress within a shorter stack, a comparative investigation is conducted under identical configurations and operating conditions, aiming to identify any unique features only appearing in the longer stacks as well as the stack size effects.In addition, the short stack consisting of three cells may represent the long stack's geometrical features, represented by the top-and bottom-end cells and the middle cell, as shown in Figure 10.The comparison revealed that the 3-cell stack has an average current density of 0.641 Acm −2 and reaches its highest temperature of 887 °C, whereas the 20-cell stack has an average current density of 0.634 Acm −2 and reaches its highest temperature of 903 °C.As shown in Figure 11, with the temperature scales adjusted to the same range, the area of the elevated temperatures of the 20-cell stack is mainly located at the upper cells compared to the 10-cell stack due to the bigger gas flow rate of the pile and the fact that the main electrochemical reactions occur in the upper half of the stack.The comparison revealed that the 3-cell stack has an average current density of 0.641 Acm −2 and reaches its highest temperature of 887 • C, whereas the 20-cell stack has an average current density of 0.634 Acm −2 and reaches its highest temperature of 903 • C. As shown in Figure 11, with the temperature scales adjusted to the same range, the area of the elevated temperatures of the 20-cell stack is mainly located at the upper cells compared to the 10-cell stack due to the bigger gas flow rate of the pile and the fact that the main electrochemical reactions occur in the upper half of the stack.The comparison revealed that the 3-cell stack has an average current density of 0.641 Acm −2 and reaches its highest temperature of 887 °C, whereas the 20-cell stack has an average current density of 0.634 Acm −2 and reaches its highest temperature of 903 °C.As shown in Figure 11, with the temperature scales adjusted to the same range, the area of the elevated temperatures of the 20-cell stack is mainly located at the upper cells compared to the 10-cell stack due to the bigger gas flow rate of the pile and the fact that the main electrochemical reactions occur in the upper half of the stack.C respectively.The maximum temperature, the minimum temperature, and average temperature of all cells except the bottom one gradually increase along the upwind direction of the stack.This may result from the fact that, in the same inlet manifold size and arrangement, the longest 20-cell stack has the biggest gas inlet flow velocity, resulting in a slower electrochemical reaction in the bottom cell, and the lowest temperature as heat is quickly removed.Figure 13 shows a comparison of the thermal stress distribution among three stacks with different numbers of cells under the same conditions.The thermal stress distribution between the 3-cell and 10-cell stacks shows similarities, while the thermal stress in the electrolyte layer at the bottom of the 20-cell stack is significantly smaller; the higher the cell location, the greater the volume ratio of high thermal stress.Comparing the bottom-end, middle, and top-end cells of the three stacks, it can be found that there is no considerable distinction between the representing cell layers; the thermal stress in the top-end layer of the 10-cell stack is significantly greater compared to the middle and bottom layers, with a difference of 9.5% in the top and middle layers.This As mentioned above, the longer stacks with a big number of single cells have higher temperatures and larger temperature differences between the cells due to their more complex structure with poor heat dissipation.Therefore, the modelling and analysis of the heat transfer and temperature distribution in kW-scaled stacks is necessary, such as the one in this study.
The temperature distribution of the electrolyte layers is compared along the central line of the symmetry plane, following the path of the primary hydrogen flow, as shown in Figure 12d-f.The predicted temperature gradually increases with the increased L (here, L refers to the length of electrolyte along the direction of the hydrogen flow) and then gradually decreases.This feature of the temperature distribution is consistent with that observed in the cloud-style plot of the temperature distribution in the co-flow mode of the 10-cell stack presented in Figure 6a.
Figure 13 shows a comparison of the thermal stress distribution among three stacks with different numbers of cells under the same conditions.The thermal stress distribution between the 3-cell and 10-cell stacks shows similarities, while the thermal stress in the electrolyte layer at the bottom of the 20-cell stack is significantly smaller; the higher the cell location, the greater the volume ratio of high thermal stress.
Figure 13 shows a comparison of the thermal stress distribution among three stacks with different numbers of cells under the same conditions.The thermal stress distribution between the 3-cell and 10-cell stacks shows similarities, while the thermal stress in the electrolyte layer at the bottom of the 20-cell stack is significantly smaller; the higher the cell location, the greater the volume ratio of high thermal stress.Comparing the bottom-end, middle, and top-end cells of the three stacks, it can be found that there is no considerable distinction between the representing cell layers; the thermal stress in the top-end layer of the 10-cell stack is significantly greater compared to the middle and bottom layers, with a difference of 9.5% in the top and middle layers.This Comparing the bottom-end, middle, and top-end cells of the three stacks, it can be found that there is no considerable distinction between the representing cell layers; the thermal stress in the top-end layer of the 10-cell stack is significantly greater compared to the middle and bottom layers, with a difference of 9.5% in the top and middle layers.This difference is more clearly predicted in the 20-cell stack, with the highest thermal stress of 156 MPa in the top-end cell and 126 MPa in the bottom-end cell, i.e., a difference of 23.8%, as shown in Figure 14.According to Figure 15, the average thermal stress in the electrolyte layer of the 3-cell stack is lower in the middle layer and slightly greater in the bottom-end and top-end cells, but the variance is minimal.The distribution in the 10-cell stack resembles that of the 3cell stack, but the variation is larger, with the lowest in the middle layer being 92 MPa and the highest in the top-end cell being 100 MPa, i.e., a difference of 8.7%, and the lowest is in the bottom-end cell of the 20-cell stack (93 MPa) and the highest in the top-end layer (112 MPa), i.e., a difference of 20.4%.In addition, the volume ratio with high thermal stress (above 100 MPa) in the top-end electrolyte layer is 5.3% for the 3-cell stack, 19.7% for the According to Figure 15, the average thermal stress in the electrolyte layer of the 3-cell stack is lower in the middle layer and slightly greater in the bottom-end and top-end cells, but the variance is minimal.The distribution in the 10-cell stack resembles that of the 3-cell stack, but the variation is larger, with the lowest in the middle layer being 92 MPa and the highest in the top-end cell being 100 MPa, i.e., a difference of 8.7%, and the lowest is in the bottom-end cell of the 20-cell stack (93 MPa) and the highest in the top-end layer (112 MPa), i.e., a difference of 20.4%.In addition, the volume ratio with high thermal stress (above 100 MPa) in the top-end electrolyte layer is 5.3% for the 3-cell stack, 19.7% for the 10-cell stack, and 92.4% for the 20-cell stack.In comparison to the 3-cell stack, the 20-cell stack has an about 17 times increased volume ratio of the high thermal stress in the top-end electrolyte layer.This finding is mainly because the longer stacks have a higher gas flow velocity in the entrance, resulting in slower electrochemical reactions and less heat produced.Furthermore, due to the complexity of the longer stack structure, the heat dissipation is much poorer, which results in higher overall temperatures and a bigger temperature difference between the cells.It is also a fact that, due to the structural constraints, the variation in the thermal stress distributed in the longer stacks is significantly greater than that in the shorter stacks.In view of the above findings and discussion, it is reasonable to point out that as the quantity of the cells rises, the average current density gradually decreases, while the maximum temperature and thermal stress increase gradually, as does the variation between cells.It would be inaccurate or even unreasonable if the findings and conclusions drawn for small-sized and shorter stacks are directly applied or extrapolated to large-sized longer stacks.

Conclusions
This study presents the CFD model for the planar designed full-size SOFC stack.The effects of different stack flow configurations on the stack performance and the thermal stress distributions are evaluated, together with the corresponding failure probabilities.The forecasted outcomes are also contrasted with the representing cells in the shorter and longer stacks to identify the stack size effect.The principal discoveries include the following: (1) In identical situations, the counter-flow stack has the highest power output together with the greater forecasted thermal stress; the cross-flow stack has the lowest power output accompanied by the highest thermal stress.The electrolyte layer experiences greater thermal stress, and the volume ratio with the thermal stress beyond 100 MPa on the electrolyte layer is the highest in the counter-flow stack and the lowest in the cross-flow stack (the reduction in the volume ratio is about 30% in the cross-flow stack).The thermal stress distribution is primarily influenced by the stack air flow direction under the same conditions.(2) The highest and lowest temperatures vary similarly between the unit cells, and the highest temperature appears in the counter-flow stack, and the lowest temperature is slightly higher in the co-flow stack than in the other stacks.The HTSR on the electrolyte layer in each cells decreases first and then increases, reaching a maximum Furthermore, due to the complexity of the longer stack structure, the heat dissipation is much poorer, which results in higher overall temperatures and a bigger temperature difference between the cells.It is also a fact that, due to the structural constraints, the variation in the thermal stress distributed in the longer stacks is significantly greater than that in the shorter stacks.In view of the above findings and discussion, it is reasonable to point out that as the quantity of the cells rises, the average current density gradually decreases, while the maximum temperature and thermal stress increase gradually, as does the variation between cells.It would be inaccurate or even unreasonable if the findings and conclusions drawn for small-sized and shorter stacks are directly applied or extrapolated to large-sized longer stacks.

Conclusions
This study presents the CFD model for the planar designed full-size SOFC stack.The effects of different stack flow configurations on the stack performance and the thermal stress distributions are evaluated, together with the corresponding failure probabilities.The forecasted outcomes are also contrasted with the representing cells in the shorter and longer stacks to identify the stack size effect.The principal discoveries include the following: (1) In identical situations, the counter-flow stack has the highest power output together with the greater forecasted thermal stress; the cross-flow stack has the lowest power output accompanied by the highest thermal stress.The electrolyte layer experiences greater thermal stress, and the volume ratio with the thermal stress beyond 100 MPa on the electrolyte layer is the highest in the counter-flow stack and the lowest in the cross-flow stack (the reduction in the volume ratio is about 30% in the cross-flow stack).The thermal stress distribution is primarily influenced by the stack air flow direction under the same conditions.
(2) The highest and lowest temperatures vary similarly between the unit cells, and the highest temperature appears in the counter-flow stack, and the lowest temperature is slightly higher in the co-flow stack than in the other stacks.The HTSR on the electrolyte layer in each cells decreases first and then increases, reaching a maximum value in the 10th cell (the top-end cell).(3) The longer and shorter stacks have similar temperature and thermal stress distributions, but their difference between the unit cells of the longer stack is bigger than that predicted for the shorter stack.It would be inaccurate or even unreasonable if the findings and conclusions drawn for small-sized and shorter stacks are directly applied or extrapolated to large-sized longer stacks.(4) The probability of failure of the porous cathodes significantly exceeds that of the porous anodes and electrolyte layers; the probability of failure of the counter-flow stacks is the highest, particularly in the longer stack.
This study is based on a coupled CFD model of the multi-physics fields in SOFC stacks, aiming to simulate and analyze the distribution of temperature and thermal stress as well as the influencing factors.Due to the difficulty in measuring the thermal stress of the stack under complex operating conditions, the simulation method proposed in this paper may serve as a good reference for predicting the electrochemical performance and thermal stress of the stack.However, the model does have some limitations and shortcomings, which are worthy of further investigations, e.g., (1) The accuracy of the model has been validated by comparing with experimental I-V curves, which needs further improvement.Future research could incorporate temperature distribution within the stack as the modeling validation.(2) The focus of this paper is on the distribution of thermal stress in the stack under steady-state operation.Future research could be extended to capture the transient operating performance of the stack with different stack configuration parameters and start-stop operations.

Figure 1 .
Figure 1.Geometry of modeled SOFC stack with (a) co-and count-flow mode; (b) cross-flow mode; (c) its unit cell with enlarged PEN and active region.

Figure 1 .
Figure 1.Geometry of modeled SOFC stack with (a) co-and count-flow mode; (b) cross-flow mode; (c) its unit cell with enlarged PEN and active region.

Figure 2 .
Figure 2. Schematic diagram of the grid and meshing in the exampled SOFC unit cell.

Figure 3 .
Figure 3. (a) Grid independence testing; (b) predicted I-V characteristics compared with the experimental one.

Figure 2 .
Figure 2. Schematic diagram of the grid and meshing in the exampled SOFC unit cell.

Figure 2 .
Figure 2. Schematic diagram of the grid and meshing in the exampled SOFC unit cell.

Figure 3 .
Figure 3. (a) Grid independence testing; (b) predicted I-V characteristics compared with the experimental one.

Figure 3 .
Figure 3. (a) Grid independence testing; (b) predicted I-V characteristics compared with the experimental one.

Figure 5 .
Figure 5. Normalized mass flow rates predicted for (a) different cells in the co-flow stack; (b) first layer in co-/counter-/cross-flow stacks.

Figure 5 .
Figure 5. Normalized mass flow rates predicted for (a) different cells in the co-flow stack; (b) first layer in co-/counter-/cross-flow stacks.

Figure 6
Figure 6 illustrates the temperature distribution of various flow stack configurations.There are conspicuous differences in temperature distribution under different flow conditions.Temperatures are focused at the central region of the stack and decrease progressively towards the periphery.In the co-flow stack configuration, the maximum temperature reaches 895 • C, with the high-temperature region situated at the center of the stack.This is because the electrochemical reactions in the stack mainly occur against the inlet region, and this heat is further carried to the center of the stack through the flowing of the generated products.The highest temperature of the counter-flow stack is 905 • C, which is the highest within the three flow modes studied.The predicted high-temperature region primarily resides within the hydrogen inlet area, attributable to the air flow rate being three times greater than that of hydrogen, and the specific heat capacity of the products on the air side is also greater than that on the hydrogen side, so most of the heat is carried forward to the hydrogen inlet, forming a high-temperature zone; the highest temperature of the cross-flow mode stack is 896 • C, and the high-temperature zone is primarily situated in the central region near the air outlet, which is also caused by the gas flow effect.The flow direction of the gases, the electrochemical exotherm, and the thermal coupling between the fluid and the solid ultimately determine the temperature field distribution of the SOFC stack.

Figure 8 .
Figure 8.(a) Normalized mass flow rates; (b) highest and the lowest temperature; and (c) volume ratio with thermal stress beyond 100 MPa predicted for individual cell electrolyte layers in different flow stacks.

Figure 8 .
Figure 8.(a) Normalized mass flow rates; (b) highest and the lowest temperature; and (c) volume ratio with thermal stress beyond 100 MPa predicted for individual cell electrolyte layers in different flow stacks.

Figure 9 .
Figure 9. Failure probability predicted for (a) volume ratio and probability of failure contribution; (b) anode; (c) electrolyte layer; (d) cathode.The stack designed with different flow configurations.

Figure 9 .
Figure 9. Failure probability predicted for (a) volume ratio and probability of failure contribution; (b) anode; (c) electrolyte layer; (d) cathode.The stack designed with different flow configurations.
Energies 2024, 17, x FOR PEER REVIEW 20 of 26 effects.In addition, the short stack consisting of three cells may represent the long stack's geometrical features, represented by the top-and bottom-end cells and the middle cell, as shown in Figure 10.

Figure 10 .
Figure 10.Bottom-end, middle, and top-end cells in three-cell stack represent the typical ones in a multi-cell long stack.

Figure 11 .
Figure 11.Temperature distribution predicted for (a) 3-cell; (b) 10-cell; and (c) 20-cell stack.As shown in Figure12a-c, the highest, lowest, and average temperatures of the 3-cell stack are basically the same (about 886 °C), while those in the 10-cell stack differ by 5 °C, 10 °C, and 4 °C, respectively, and those in the 20-cell stack by 17 °C, 24 °C, and 15 °C respectively.The maximum temperature, the minimum temperature, and average tem-

Figure 10 .
Figure 10.Bottom-end, middle, and top-end cells in three-cell stack represent the typical ones in a multi-cell long stack.

Energies 2024 ,
17, x FOR PEER REVIEW 20 of 26effects.In addition, the short stack consisting of three cells may represent the long stack's geometrical features, represented by the top-and bottom-end cells and the middle cell, as shown in Figure10.

Figure 10 .
Figure 10.Bottom-end, middle, and top-end cells in three-cell stack represent the typical ones in a multi-cell long stack.

Figure 11 .
Figure 11.Temperature distribution predicted for (a) 3-cell; (b) 10-cell; and (c) 20-cell stack.As shown in Figure12a-c, the highest, lowest, and average temperatures of the 3-cell stack are basically the same (about 886 °C), while those in the 10-cell stack differ by 5 °C, 10 °C, and 4 °C, respectively, and those in the 20-cell stack by 17 °C, 24 °C, and 15 °C respectively.The maximum temperature, the minimum temperature, and average tem-

Energies 2024 ,
17, x FOR PEER REVIEW 21 of 26L refers to the length of electrolyte along the direction of the hydrogen flow) and then gradually decreases.This feature of the temperature distribution is consistent with that observed in the cloud-style plot of the temperature distribution in the co-flow mode of the 10-cell stack presented in Figure6a.

Figure 12 .
Figure 12.Temperature distribution in 3-cell and 10-cell stacks: (a) highest temperatures; (b) lowest temper-atures; (c) average temperature; and (d) bottom-end cell; (e) middle cell; (f) top-end cell along the hydrogen main flow direction predicted for representing cells (in the figure, L refers to the length of electrolyte along the flowing direction of hydrogen gas).

Figure 12 .
Figure 12.Temperature distribution in 3-cell and 10-cell stacks: (a) highest temperatures; (b) lowest temper-atures; (c) average temperature; and (d) bottom-end cell; (e) middle cell; (f) top-end cell along the hydrogen main flow direction predicted for representing cells (in the figure, L refers to the length of electrolyte along the flowing direction of hydrogen gas).

Energies 2024 ,
17, x FOR PEER REVIEW 22 of 26difference is more clearly predicted in the 20-cell stack, with the highest thermal stress of 156 MPa in the top-end cell and 126 MPa in the bottom-end cell, i.e., a difference of 23.8%, as shown in Figure14.

Table 1 .
Geometric characteristics of the modeled cell.

Table 2 .
Parameters related to electrochemical reactions.

Table 3 .
Physical parameters of SOFC stack components.

Table 4 .
Weibull parameters for the SOFC materials utilized in the current study.

Table 5 .
The parameters estimated in this study.

Table 6 .
Boundary conditions of SOFC stack model.

Table 7 .
Grid type of SOFC stack.

Table 8 .
Proportion of different flow patterns exceeding corresponding thermal stress regions.

Table 8 .
Proportion of different flow patterns exceeding corresponding thermal stress regions.