Temperature Gradient Analyses of a Tubular Solid Oxide Fuel Cell Fueled by Methanol

Thermal management in solid oxide fuel cells (SOFC) is a critical issue due to non-uniform electrochemical reactions and convective flows within the cells. Therefore, a 2D mathematical model is established herein to investigate the thermal responses of a tubular methanol-fueled SOFC. Results show that unlike the low-temperature condition of 873 K, where the peak temperature gradient occurs at the cell center, it appears near the fuel inlet at 1073 K because of the rapid temperature rise induced by the elevated current density. Despite the large heat convection capacity, excessive air could not effectively eliminate the harmful temperature gradient caused by the large current density. Thus, optimal control of the current density by properly selecting the operating potential could generate a local thermal neutral state. Interestingly, the maximum axial temperature gradient could be reduced by about 18% at 973 K and 20% at 1073 K when the air with a 5 K higher temperature is supplied. Additionally, despite the higher electrochemical performance observed, the cell with a counter-flow arrangement featured by a larger hot area and higher maximum temperature gradients is not preferable for a ceramic SOFC system considering thermal durability. Overall, this study could provide insightful thermal information for the operating condition selection, structure design, and stability assessment of realistic SOFCs combined with their internal reforming process.


Introduction
The solid oxide fuel cell (SOFC) is a promising future energy device since electricity can be electrochemically generated from the chemical energy of fuels through SOFCs with higher efficiency and lower greenhouse gas emission in comparison to traditional thermal power plants [1]. SOFCs are typically operated at high temperatures (873-1273 K) to maintain moderate oxygen ion conductivity of ceramic electrolytes [2], thus offering various beneficial aspects to the SOFC systems. For example, cheaper nickel (Ni) metal could be utilized as the chemical/electrochemical catalyst in SOFCs, rendering a more economical prospect than protonexchange membrane fuel cells (PEMFCs) [3][4][5]. Besides, the utilization of hydrocarbon [6,7], alcohols [8], ammonia [9,10], biomass [11], or even solid carbon [12,13] in SOFCs is possible for power generation because these fuel sources can be internally reformed into more effective fuels (H 2 and CO) for faradic processes under the anodic environment [14]. The internal reforming process provides SOFCs with an efficient and sustainable operation since the thermal coupling of endothermic reforming and exothermic processes can improve the cell efficiency [15], and the integration of renewable fuel productions with direct hydrocarbon or alcohol SOFCs is promising in achieving carbon-neutral power generation. Moreover, given the considerable amount of waste heat produced due to cell inefficiency, accounting for 20-50% LHV (Lower heating value) of fuels [16], exhausted gases could be further used for bottoming cycles to generate extra power, maximizing overall system efficiency [17,18]. However, despite the numerous benefits involved in the SOFC system, thermal management has been a major challenge [19]. Owing to the uneven electrochemical reactions and heat transfer of gas convective flows, non-uniform temperature distributions were frequently observed within the SOFC assembly or stack, especially for the internal reforming SOFCs, posing a threat to material deterioration and structural integrity. Local temperatures higher than the design value could lead to nickel coarsening because of the relatively low melting temperature of nickel and thermally favorable sintering process [20], degrading anode functions and consequently lowering the cell performance. For example, electrode sintering was observed at the bottom part of a tubular SOFC operated on the methane flame mode as the temperature of the corresponding position reached 1300 K [21]. Another thermal issue of temperature gradient arising from the non-uniform temperature distribution negatively affects the cell structure, thus cracking or delaminating brittle ceramic components that are more vulnerable under tensile stress [22]. For example, delaminations between the electrolyte and the electrode as well as cracks in the electrolytes were found in experimental tests due to the temperature gradient [21,23], likely resulting in the complete failure of SOFCs as their components cannot be replaced or repaired [24]. Therefore, under the complicated thermal environment of SOFCs, the cell components are subjected to the accelerated rate of property evolution and degradation, eventually causing the fuel cell failure.
To solve the aforementioned problems, various efforts have been devoted to developing effective temperature control strategies in recent years. Zeng et al. [21] experimentally fabricated a tubular SOFC with a liquid-sodium heat pipe to equalize the temperature distribution. It was found that with the help of high thermal conductivity induced by the evaporation-transport-condensation process in liquid sodium, the axial temperature gradient was significantly decreased, enhancing the cell electrochemical performance and life span. Dillig et al. [16] and Marocco et al. [25] have drawn similar conclusions in their simulation investigations. Recently, Promsen et al. [26] proposed a novel concept of SOFC cooled by saturated water. The calculation results indicated that the water-cooled stack demonstrated a better performance due to a more uniform temperature profile. Besides, due to the efficient cooling of the saturated water, the air flow was greatly reduced since conventionally heat dissipation in an SOFC is conducted by oversupplying cathode gas, decreasing the parasitic loads related to air preheating and blower energies. Combining an SOFC with a heat pipe or a water cooling tube enhances the temperature uniformity and electrochemical performance of a cell and reduces the costs related to auxiliary devices. However, the manufacturing cost and system complexity remain a challenge, especially for the SOFC stacks, as numerous heat pipes must be fabricated to achieve sufficient cooling effect [25]. Besides, possible leakage of liquid metal from the heat pipe might catastrophically damage the cell system given the complicated gas environment and chemically reactive nature of the liquid metal. Based on the same approach of increasing the thermal conductivity of solid cell components, utilization of metal support [27,28], metallic interconnect [29], or metal cathode flow distributor [30] in SOFC fabrication also exhibited a more uniform temperature distribution, leading to improved cell stability in terms of thermal stress [29]. However, as indicated from the numerical study conducted by Park et al. [31], the bonding layer between the ceramic cell and the metallic support resisted hydrogen mass transfer, resulting in a 17% decrease in the average current density compared to the cermet anode-supported SOFC. Besides, possible metal oxidation risk [32], accelerated poisoning of cathode catalyst by Cr from ferritic stainless steel [33], as well as relatively poor electrochemical performance [34] have limited the commercialization of this technology.
However, despite numerous theoretical and experimental works conducted with a focus on thermal control and management, to the best of our knowledge, a detailed thermo-electrochemical simulation has not been performed to identify the thermal characteristics, especially the temperature gradient profile of SOFC operated on methanol fuel. Importantly, unlike the battery, SOFC is an open system that thermally interacts with reacting gases; thus, the cell thermal state can be considerably influenced by gas conditions [35], especially for cathode air, but the cooling effects and thermal/electrical responses of air conditions have not been fully understood yet. Therefore, to fill this research gap, a validated mathematic model extended from our previous simulation works [36,37] is developed to further investigate the temperature gradient profiles in cell components, especially ceramic electrolytes. This work will be critically essential for the operating condition selection, structure design, and failure assessment of realistic SOFCs.

Working Principles and Modeling Assumptions
A fuel mixture of methanol and steam (molar ratio 1:1) is supplied into the fuel channel to electrochemically generate electricity. As shown in Fig. 1, the computation domain thermally consists of the anode and cathode channels, porous Ni-YSZ and YSZ-LSM composite electrodes, as well as dense YSZ electrolyte [38]. Macro-and micro-structural parameters, as well as the physical properties of materials adopted in the present study, are summarized in Tables 1 and 2, respectively.
Methanol is used as the primary fuel in this study due to its several promising features, such as storage/transportation  [39], relatively high volumetric energy density [40], and significantly low coking threat [8]. To further eliminate carbon formation in cell operation, an equimolar amount of steam is added to the main methanol stream [6,14]. Importantly, under the active catalytic effect of nickel, methanol could be readily converted into a H 2 and CO mixture through the methanol decomposition reaction (MDR, Eq. (1)), followed by water gas shift reaction (WGSR, Eq. (2)) to further generate H 2 in the presence of water [39]. At the anode functional layer, the generated H 2 and CO could be oxidized into H 2 O and CO 2 by oxygen ions (Eqs. (3) and (4)), produced at the cathode side (Eq. (5)) and subsequently transported to the anode through an electron-insulated electrolyte. Electrons generated from electrochemical oxidation reactions can be conducted to the cathode layer via an external circuit producing electrical power.
Based on the working mechanism, assumptions adopted in the simulation are listed below: (1) Due to kinetically fast reaction rates, only hydrogen and carbon monoxide could be electrochemically oxidized. (2) Local temperature equilibrium hypothesis is applied for porous electrodes because of the negligible temperature difference in the gas/solid phases [41].

Governing and Constitutive Equations
Based on the chemical/physical processes within the cell, five sub-models will be calculated simultaneously to elucidate the thermo-electrochemical environment of methanolfueled tubular SOFC.

Chemical Reaction Model
Chemical reaction rates and heat changes could be calculated by the chemical reaction model. As a result of the active catalytic activity of Ni, MDR, and WGSR will principally occur at the porous anode [42]. The corresponding reaction rates (R r ) could be expressed by [43,44]:  where T and p l i represent temperature and species local partial pressure, respectively; R denotes the universal gas constant, and tuning parameter k D is used to validate the model.
Chemical reactions could lead to heat generation or consumption according to thermodynamics; thus, total chemical heat changes (Q che ) are determined by both reaction enthalpy changes of formation (ΔH) and rates, as shown below: where h i and C p,i denote the species enthalpy of formation and heat capacity, respectively.

Electrochemical Reaction Model
Electrical output and reversible/irreversible heat generations related to the electrochemical processes are determined by the electrochemical reaction model. The current density is described by the general Butler-Volmer (BV) equation (Eq. (17)), wherein the activation overpotential (η act ) reflects the electrochemical activation barrier [45].
Here, α, n, and F are the charge transfer coefficient, the number of electrons for one single electrochemical reaction, and the Faraday constant, respectively; i 0 represents the exchange current density. In the present model, both H 2 and CO are oxidized, and the electrochemical rate of hydrogen oxidation is faster than that of carbon monoxide [42]; thus, i 0 is set as 2000, 3000, and 5300 A/m 2 for O 2 , CO, and H 2 , respectively. Because of various overpotential losses, the output potential (V, Eq. (18)) could be much lower than the thermodynamic equilibrium potential (E, Eqs. (19) and (20)), especially under high current density conditions.
In addition to the activation overpotential, ohmic overpotentials for ionic/electronic transport resistances can be calculated by Eqs. (23) and (24) according to Ohm's law [46].
Here, Ø and σ eff represent the phase potential and effective conductivity, respectively. During the electricity output, heat is released simultaneously, involving reversible electrochemical heat generation (Q ele ) due to change of entropy, irreversible joule heating (Q ohmic ) because of transport of charged species, and irreversible heat generation (Q act ) upon activation polarization [19]. Therefore, corresponding equations can be written as:

Computational Fluid Dynamics Model
Fuel and gas flows described by the velocity (U) and pressure (p) fields in electrodes and gas channels will be simulated by the computational fluid dynamics (CFD) model. In CFD, gas behaviors in the channels can be described by general continuity and Navier-Stokes (NS) equations (Eqs. (32) and (33)) [47]. However, in fuel and air electrodes, the modified continuity equation (Eq. (34)) with a mass source term and modified NS equation (Eq. (35)) with Darcy's term account for the mass exchange and porous microstructure, respectively [48].
Here, ε and κ represent the porosity and permeability of electrodes, respectively; Q m is the net mass change; the dynamic viscosity (μ, Eq. (36)) of the gas mixture is estimated based on the species molecular weight (M i ) and molar fraction (y i ).

Mass Transfer Model
Owing to various physical/chemical processes in an SOFC, the mass transfer model is used to calculate the species concentration distributions via steady-state species molar conservation equations with and without a source term applied in electrodes (Eq. (38)) and channels (Eq. (37)), respectively.
Here, R i and J i are the molar change rate and molar diffusion flux of species i, respectively; u and c represent the molar average velocity and molar concentration of the gas mixture, respectively. Molar diffusion flux in the species molar conservation equations could be estimated by the Stefan-Maxwell approach (Eqs. (39) and (40)) [49]. Notably, the effective binary diffusion coefficient ( D eff ij , Eq. (42)) corrected by molecular and Knudsen diffusions, as well as microstructure [50] is used in porous electrodes [51], while only binary diffusion coefficient (D ij ), which can be determined by Eq. (41) [52,53] is considered in gas channels. As claimed by the kinetic theory [54], the Knudsen diffusion coefficient (D ik ) could be expressed as Eq. (43).

Heat Transfer Model
Thermal characteristics within the entire computational domain can be described by the heat transfer model. To demonstrate the complicated thermal processes in an SOFC, general energy conservation equations with and without a heat source term are applied in electrodes (Eq. (45)) and gas (44)), respectively. Besides, the mole fraction averaging method is used to estimate both the heat capacity (C p , Eq. (46)) and thermal conductivity (λ g , Eq. (47)) of a gas mixture [52,55]. The gas species thermodynamic properties could be found in our previous publication [36].

Mathematical Method and Operating Conditions
The developed model is simulated by solving fully-coupled nonlinear equations through the finite element method with a relative tolerance of 10 -3 . To ensure grid independence, the model is composed of 30, 150, 10, 40, and 30 elements for anode channel, anode, electrolyte, cathode, and cathode channel widths, respectively, as well as 200 for cell length, generating 52,000 elements. Besides, the validation is conducted by comparing the calculated polarization curves with the experimental data [34], which can be readily found in our published studies [32,33]. Working conditions for the following parametric calculation are given in Table 3.

Effects of Air Flow Rate
Calculations are employed to examine the effects of the air flow rate on tubular cell thermal/electrical performance under several typical operating temperatures. Here, the operating potential is set as 0.5 V because although a large current density or a low operating potential is not desired in practical operations, negative thermal features can become quite obvious under a considerably high current density.
Current densities (Fig. 2a) can be highly influenced by the air flow rate. At a relatively low flow rate, current densities are nearly the same under different temperature conditions and considerably increase as the air flow rate increases. It is mainly because the cell performance is restricted by an insufficient oxygen supply for electrochemical reactions, evident from the complete consumption of oxygen (Fig. 2b), especially near the cathode outlet (Fig. 2c). With a further increase to 1200 SCCM, the current density at 1073 K continues to enhance slightly, while the power outputs for three relatively low-temperature conditions experience a gradual  [59]. The possible reason for these two opposite phenomena is because of different cell sensitivities to the temperature change. At low operating temperatures, the cell performance is more sensitive to temperature variation, as indicated in Fig. 2a, where a temperature increase by 50 K from 873 K leads to a substantial current density improvement (about 5800 A/m 2 ), while the cell electrochemical performance at 1073 K is just 1400 A/m 2 higher compared to that at 973 K. Therefore, despite the increased oxygen concentration and thus decreased overpotential, the decrease in the cell average temperature (Fig. 2d) caused by the cooling air dominates the power output at 873 K, 923 K, and 973 K, reducing the cell performance. Conversely, a high operating temperature weakens the negative effect of the average temperature decrease, resulting in an insignificant but continuous increase in the power output when the air flow rate is faster than 500 SCCM at 1073 K.
In addition, the air flow rate can strongly impact the cell thermal characteristics involving axial temperature gradient and temperature distributions in the ceramic electrolyte (Fig. 3a, b). At 873 K, a peak temperature gradient approximately occurs at the center of the electrolyte because of the relatively moderate current density (Fig. 2a) and a smooth increase in temperature along the cell length near the gas inlet (Fig. 3a). Although limited by the inadequate oxygen and low chemical/electrochemical kinetics due to a slow air flow and a low operating temperature, respectively, the current density is only about 11,000 A/m 2 at the air flow rate of 200 SCCM, the temperature gradient still reaches 43.5 K/ cm because of poor heat convection and easy heat accumulation. Fortunately, with the rising air flow rate and enhanced flow convective heat transfer capacity, the temperature gradient can be effectively suppressed within 20 K/cm, thus reducing the possibility of structural failure during practical SOFC operation. Besides, the reduction in current density also plays a role in controlling the high-temperature gradient (Fig. 2a), especially at a high air flow rate, since various heat sources strongly depend upon the current density. Unlike the distribution at 873 K, the largest axial temperature gradient likely exists near the fuel inlet when the fuel cell is operating at 1073 K (Fig. 3b), primarily because the high operating temperature benefits the generation of direct fuels (hydrogen and carbon monoxide), electrochemical reactions, and YSZ electrolyte ion transport, resulting in an extremely elevated current density and thus a huge temperature rise near the fuel inlet (Fig. 3b). Although the increased air flow rate ensures sufficient oxygen supplement (Fig. 2c) and enhanced power generation (Fig. 2a), the induced strong heat convection still promotes heat dissipation under 1073 K, slightly improving the temperature uniformity and reducing the peak temperature gradient.
As the air flow rate continuously increases, mitigations of high temperature and peak temperature gradient become negligible for both operating temperature conditions ( Fig. 3a  and b), indicating that heat convective capacity arising from the air flow can show a limited cooling effect on the exothermic cell when the air flow rate exceeds a certain value. Therefore, to better evaluate the effectiveness of this active cooling strategy, an excess air ratio is introduced to represent the actual air flow rate [26]. As observed in Fig. 4, all peak temperature gradients occur in two stages with the increase in excess air ratio (air flow rate: 200-2400 SCCM), involving sudden and then gradual declines. At 1073 K, tripling the air ratio from 2 to 6 merely results in about an 18% reduction in the maximum temperature gradient, indicating that the local high-temperature gradient cannot be efficiently eliminated to a safe value (less than 10 K/cm [15,19]) by simply increasing the airflow rate at such a high current density (0.5 V), especially for high-temperature conditions, which is consistent with Promsen's statement [26]. The underlying reason for the airflow's poor cooling effect could be the restricted heat convection between the cell assembly and the convective gas flow in the air channel, as evidenced by the heat transfer mechanism at the cathode [35]. Unlike the anode, where the heat generated from the thin interface is principally transported by mass transfer. Since the convective flow resulting from the mass exchange moves away from the anode functional layer [60], heat conduction serves as the predominant form of heat transfer in the cathode as the air flows toward the cathode functional layer. Besides, the reported dimensionless value of π z , a parameter comparing heat convection with conduction, is 1.075 for the anode, three orders higher as that for the cathode. Therefore, simply accelerating the convective flow in the air channel is not the most effective strategy to control temperature and temperature gradient within the safe region. Besides, an additional air flow normally increases the cost [59] since, reportedly, over 10% of the total cell electricity generation will be consumed by the air blower devices [16,21,61]. Therefore, the air flow rate must be carefully controlled considering material deterioration, structural integrity, and cell efficiency.

Effects of Operating Potential
As the key factor, the operating potential can have a dramatic effect on the cell's thermal characteristics. The effects of operating potential at different temperatures are illustrated in Fig. 5. Due to the endothermic MDR, the electrolyte temperature initially decreases and subsequently increases because of the synergistic promotional effects between local current density and actual cell temperature downstream. A local cooling spot near the gas inlet results in an abrupt fall of the axial temperature gradient, especially at a relatively high operating potential (Fig. 5a, b). The low localized temperature could negatively influence the current density (Fig. 5c, d), contributing little to the cell's electrical performance, and the induced temperature gradient (− 60 K/cm) does not fall within the safe region for a long-term stable operation based on the recommendation value (10 K/cm) of maximum temperature gradient [15,19], thus increasing the failure possibility related to the structural defects, such as cracks and delamination due to the applied thermal stress on the brittle ceramic structure. Fortunately, the peak temperature gradient caused by the internal reforming operation could be effectively inhibited by properly adjusting the operating potential as a free condition of extremely high positive or negative peak temperature gradient could be achieved between 0.5 and 0.6 V at 873 K, and 0.6 and 0.7 V at 1073 K (Fig. 5a, b), thereby achieving a localized thermal neutral state. Therefore, an abnormal temperature gradient could be avoided by selecting suitable operating conditions for the fuel cell without sacrificing significant power generation since excessive heat generated from the cell's inefficiency could be consumed by the endothermic decomposition reaction, showing an efficient and attractive operation of direct internal reforming in an SOFC.
In addition, temperature reduction near the cell inlet induced by the MDR could also be affected. A higher operating temperature and a lower potential favor the decomposition reaction (Fig. 6a), mainly arising from the higher reaction kinetics and promotional effects of hydrogen and carbon monoxide consumption by the faradic processes, respectively. Despite the high average current density at a high working temperature, a more favored endothermic reaction leads to a large temperature drop at the operating Fig. 4 Effect of excess air ratio on the maximum axial temperature gradient at the middle of the electrolyte under different temperatures potential of 0.9 V (Fig. 6b). Moreover, contrary to the published fact that cells operated on methane fuel generally experience large temperature drops, a relatively lower temperature reduction (below 30 K) is found in the present model, as demonstrated in Fig. 7. Dokmaingam et al. [62] have drawn a similar conclusion that methanol-fueled IIR-SOFC (indirect internal reforming) provided the smoothest temperature change among various promising fuels, including methane, biogas, methanol, and ethanol. The potential reason is related to different thermodynamical intensities of internal reforming reactions. Within the actual working temperature of SOFC (from 1000 to 1200 K), MDR is moderately endothermic as the associated reaction heat consumption is about 105 kJ/ mol, while methane steam (226 kJ/mol) and dry (260 kJ/ mol) reforming processes are highly endothermic [40,42]. Therefore, compared to methane reforming reactions, MDR   Fig. 5 a, b Effects of operating potential on axial temperature gradient and temperature; c, d local current density distributions at the middle of the electrolyte under different temperatures of a, c 873 K and b, d 1073 K Fig. 6 Effects of operating potential on a energy consumption via the endothermic methanol decomposition reaction; b maximum temperature reduction at the middle of the electrolyte under different temperatures is less endothermic and could lead to a lower temperature drop, which is beneficial to cell performance and durability, especially while operating under a start-up or partial load condition, showing a promising prospect for methanol as the fuel for SOFCs from a thermodynamic perspective.

Effects of Air Inlet Temperature
Besides the cathode gas flow rate, the air inlet temperature can also influence the cell performance and thermal features. Therefore, the effects of cathode gas temperature are studied at the fuel gas inlet temperatures of 973 and 1073 K.
With the rising air inlet temperature, the calculated current densities increase because of the increased average temperature of the cell assembly, as shown in Fig. 8. Besides, the current density at 973 K experiences a larger improvement than that at 1073 K because of higher cell sensitivity to temperature change, as demonstrated previously. Surprisingly, a small amount of air temperature increment can decrease the peak axial temperature gradient of the electrolyte near the inlet to some extent (Fig. 9a, b). For example, increasing the air temperature by 5 K could reduce the temperature gradients by about 18% for 973 K and 20% for 1073 K. However, with a further increase, a negative peak temperature gradient can appear because of the newly generated lowtemperature spot, especially for the large air temperature increment. The reduction in the peak temperature gradient is due to the warmer air flow serving as the heat source to heat the front end of the cell assembly, which could lead to a gradual temperature increase at the corresponding position (Fig. 9a, b), smoothing the temperature gradient distribution. In contrast, the excessively warmer air flow could cause a low-temperature spot because of the heat transfer between the air stream and the cell assembly when the temperature difference becomes larger than 30 K, creating a new negative peak temperature gradient. Meanwhile, a comparable positive peak temperature gradient could also be formed at the subsequent position at 973 K, and this is mainly because of the enhanced local current density upstream (Fig. 9c). Besides, low molar fractions of effective fuels at the anode (Fig. 9d) due to the improved cell performance cause a more rapid drop in the local current density downstream for the high increment of air temperature (Fig. 9c), thus lowering the heat generation and subsequently, the temperature gradient. Compared to the condition of 973 K, the cell with a fuel inlet temperature of 1073 K controls the temperature gradient more optimally since a minimal air temperature increase could lead to more temperature gradient reduction, and a further increment could not result in either huge negative or positive gradients because of a higher current density and reduced cell sensitivity to temperature change, respectively.
However, the increase in the air inlet temperature can have a complicated effect on the temperature gradients of other components (Fig. 10). Given the anode-supported structure used in this model, an increased air temperature increases the temperature gradient for the anode mainly because of a relatively long distance to the air channel and enhanced exothermic processes. Conversely, a larger gradient drop could be observed at the cathode since it is the component with the shortest distance to the air channel.
In general, simply increasing the temperature of the incoming air could help reduce the peak temperature gradient of the electrolyte and thus the structural failure possibility since it is frequently reported that the maximum tensile stress can exist in the electrolyte of SOFC, especially for the internal reforming operation [27-29, 63, 64], practically making the electrolyte the most vulnerable component. However, the negative effects of this approach on other important components must be paid attention to as there is a trade-off in the temperature gradient changes at different radial positions. Fig. 7 Comparison of local temperature reduction due to the internal reforming reactions between methanol and methane-fueled SOFCs. Data from Ref. [2,15,42,62,[68][69][70] Fig. 8 Effect of air temperature increment on the electrochemical performance and the average temperature of anode-electrolyte-cathode at 973 K and 1073 K

Effects of Fuel/Air Flow Arrangements
Flow arrangements of the tubular fuel cell, including coflow and counter-flow are critically important to its electrochemical characteristics and thermal management (Table 4); thus, cells with different air flow directions are simulated to assess the heat dissipation capacity at 973 K and 0.5 V. As indicated in Fig. 11, a distinctive discrepancy of temperature distributions is that a cell with a counter-flow setup exhibits a larger high-temperature area compared to the co-flow counterpart, especially in the fuel inlet region where the reforming reaction is primarily performed, thus producing more effective fuels for faradic processes, indicating a more efficient waste heat recovery for the endothermic decomposition reaction. Besides, a considerably higher average cell assembly temperature of 1102.4 K enables the counter-flow cell to perform better, leading to a higher current density  (Table 4). Therefore, the cell with the counter-flow arrangement is characterized by higher power output and thereby efficiency, which coincides with many research works [15,29,58,[65][66][67] from literature. However, a large hightemperature area is generated by the poor heat dissipation because opposite directions of air and fuel flows diminish the heat convection capacity, causing heat accumulation near the cell center. Quite different temperature distributions between these two configurations can create much disparate temperature gradient distributions, as demonstrated in Fig. 12. In addition to a higher positive peak gradient (136.66 K/cm) near the fuel inlet, a negative one (− 382.58 K/cm) is also found near the air inlet for the counter-flow setting because the axial temperature can decrease from the cell center (Fig. 12b) due to the combined effects of the cooler incoming air and the decreased local current density (Fig. 12c) evident from the  slight increase in hydrogen fuel near the air inlet (Fig. 12d).
Notably, the absolute value of the negative peak gradient is much higher than that of the positive one, mainly because the rate of temperature drop at the air inlet is higher compared to the temperature rising rate near the fuel inlet due to the anode-supported structure used in the current model and the more rapid drop of current density near the air inlet than the current density enhancement near fuel inlet (Fig. 12c). Therefore, the co-flow configuration has a safer temperature profile, while much higher temperature gradients are found in both fuel and air inlets for the counter-flow counterpart. Radial temperature and temperature gradient distributions of the cell assembly have not been discussed frequently due to negligible temperature variation along the thickness direction since relatively high thermal conductivities enable solid materials to conduct heat efficiently. For the co-flow arrangement, the radial temperature difference is less than 1 K (Fig. 13) because fuel and air flows follow the same temperature variation pattern, leading to relatively small temperature gradients at three typical positions (Fig. 14). However, the temperature along the cell thickness exhibits different tendencies at different axial positions for the counter-flow setting (Fig. 13), showing opposite temperature changes at the fuel and air inlets. Besides, a relatively high gradient intensity (− 60.5 K/cm) could be found at the interface of electrolyte and cathode near the air inlet (fuel outlet) (Fig. 14a), again because anode supportive structure enables the electrolyte to dissipate heat more rapidly through the thin cathode near the air inlet.

Conclusion
The two-dimensional axisymmetric model studies the thermal and electrochemical characteristics of tubular SOFC fed by equimolar amounts of methanol and steam. This validated numerical model fully considers the mass and heat transfer processes, hydrogen/carbon monoxide electrochemical oxidations, as well as methanol conversion and simulates the cell's thermo-fluid environment under various working conditions involving operating potential, air flow rate inlet temperature, and gas flow arrangements. The main conclusions are as follows: (1) The air flow rate can have a complicated effect on a cell's performance. Limited by the insufficient oxygen at relatively low air flow rates, current densities increase dramatically with the rising air flow rate. However, with a further increase, the power output experiences a gradual drop because of the induced temperature decrease at 873, 923, and 973 K, but a slight increase at 1073 K due to the low cell sensitivity to temperature variation. Besides, as an active cool- ing strategy, oversupplying the cathode air could not effectively suppress the detrimental thermal features to a safe level at a considerably high current density due to the intrinsic heat transfer in the cathode. (2) Owing to the thermal coupling of cell inefficiencies and the endothermic decomposition reaction, a localized thermal neutral state was achieved by suitably adjusting the operating potential and eliminating the peak axial temperature gradient of the electrolyte. Compared to methane-fueled SOFCs, an SOFC operating on methanol experiences a less temperature reduction (below 30 K) near the fuel inlet due to the less endothermic thermodynamic nature, causing fewer negative effects on the cell performance and durability. (3) Interestingly, increasing the air inlet temperature by 5 K could suppress the peak axial temperature gradients by about 18% at 973 K and 20% at 1073 K since warmer fresh air can be the heat source to heat the front end of the fuel cell, while a further increase could generate negative and positive temperature gradient peaks because of the new cooling spot and elevated current density, respectively. However, as the result of the anode's supportive structure, an increased air inlet temperature can have different effects on other components, thus increasing the corresponding failure possibility. (4) Unlike the cell with the co-flow arrangement, high temperature concentrates at the cell center for the counterflow counterpart, enhancing the electricity output and cell efficiency. However, because of poor heat dissipation capacity and enhanced heat generation, the cell with the counter-flow setting is characterized by higher local temperature gradients along the cell length and radial direction, making the counter-flow undesirable from the perspective of thermal stability.
Notably, in addition to the temperature gradient distribution, mechanical properties of cell components, including thermal expansion coefficient, elasticity, and Poisson's ratio, can affect the cell's thermal stress conditions. Therefore, the electro-thermo-mechanical modeling of SOFCs is considered in our subsequent simulation work to study the operation stability related to thermal stress.