2D Simulation for CH₄ Internal Reforming-SOFCs: An Approach to Study Performance Degradation and Optimization

Abstract

Solid oxide fuel cells (SOFCs) are a well-developed technology, mainly used for combined heat and power production. High operating temperatures and anodic Ni-based materials allow for direct reforming reactions of CH₄ and other light hydrocarbons inside the cell. This feature favors a wider use of SOFCs that otherwise would be limited by the absence of a proper H₂ distribution network. This also permits the simplification of plant design avoiding additional units for upstream syngas production. In this context, control and knowledge of how variables such as temperature and gas composition are distributed on the cell surface are important to ensure good long-lasting performance. The aim of this work is to present a 2D modeling tool able to simulate SOFC performance working with direct internal CH₄ reforming. Initially thermodynamic and kinetic approaches are compared in order to tune the model assuming a biogas as feed. Thanks to the introduction of a matrix of coefficients to represent the local distribution of reforming active sites, the model considers degradation/poisoning phenomena. The same approach is also used to identify an optimized catalyst distribution that allows reducing critical working conditions in terms of temperature gradient, thus facilitating long-term applications.

1. Introduction

Solid oxide fuel cells (SOFCs) are well-known electrochemical devices suitable for a sustainable energy production in both high capacity power plant and residential use [1,2]. Current research on SOFCs focuses on (i) performance and durability improvement through characterization of new innovative materials [3], (ii) experimental detection of degradation mechanisms [4], and (iii) modeling tools for system simulation and control [5]. Nevertheless, a relevant limit to their wider spread on the energy market is the lack of a proper network for distribution of hydrogen, the main cell reactant, at low price [6]. Several existing processes are repeatedly improved in order to obtain a more efficient and sustainable H₂ production route. For instance, steam reforming of hydrocarbons, the main approach used at industrial level, is enhanced through the introduction of new reactor designs to obtain a higher quality of outlet syngas, such as catalytic membranes or sorption enhanced steam reforming [7,8]. Partial oxidation reaction has been introduced in different applications thanks to the need for lower working temperature and high fuel conversion obtained through the introduction of self-sustained electrochemical promoted catalysts [9]. Gasification of biomass/waste has also gained literature attention in the last decades, considering the wide diversity and the easy availability of feed [10]. Water electrolysis is also becoming a competitive application through the integration with renewable sources in order to reduce requested external power and operating cost [11].

An alternative way to overcome the issue of H₂ production and distribution is the possibility to directly use light hydrocarbons such as CH₄, CH₃OH, gasoline, diesel, syngas, or biogas as SOFC fuel [12,13,14]. Differently from low temperature fuel cells, SOFCs do not require a pure H₂ feeding stream since carbon compounds such as CO and CO₂ are not poisonous for the materials used. Coupled with high operating temperatures, this allows for the in-situ production of H₂ necessary in the electrochemical mechanisms through reforming reactions. The use of alternative fuel makes SOFC safer application due to highly flammable and volatile of pure H₂. There are three possible configurations based on the integration of cells with reforming processes [15]:

• External reforming (ER): SOFCs and reforming units are two distinct blocks. In the reforming reactor, H₂ is produced from light hydrocarbons and consequently fed to the fuel cell anode inlet;
• Indirect internal reforming (IIR): SOFCs and reforming units are two neighboring blocks in order to favor thermal exchanges;
• Direct internal reforming (DIR): the reforming takes places inside SOFCs using the Ni-based anode material as catalyst.

While the first one is the easiest to operate thanks to the lack of interactions between two units, the other cases are more efficient since the requested heat for the endothermic reforming reaction is provided by the exothermic electrochemical process. At the same time, since achieved usually by flowing excess air through the cathode, the cell cooling decreases because the reforming prevents an excessive temperature rise. Moreover, capital and operating costs are reduced without separated external units [16]. Compared to IIR, DIR configuration also allows further process optimization as the continuous consumption of H₂ due to electrochemical reactions enhances the hydrocarbon conversion and results in a more uniform H₂ distribution [17]. However, this design is characterized by larger temperature gradients and relevant carbon deposition [18]. The high anodic Ni content favors a faster reforming process compared to the electrochemical one. Thus, hydrocarbon conversion is usually complete within a small distance from the inlet, leading to initial severe cooling effects. Even if proper cell materials are used in order to reduce the mismatch among thermal expansion coefficients, the subsequent temperature gradient on the fuel cell plane may induce system failure [19]. Indeed, creep formation is favored using brittle ceramics. However, it also occurs in the metallic interconnects that are characterized by a thermal behavior strongly dependent on temperature [20]. Appropriate modifications of the anodic structure, such as the partial poisoning of active sites or doping processes [21], are under investigation to decrease the reforming rate. Another common problem is the carbon deposition that is catalyzed by the presence of Ni and results in decreasing the active sites. This phenomenon is reversible and can be minimized by increasing the inlet steam to carbon ratio (S/C ratio) without a high fuel dilution to avoid the reduction of electrical efficiency [22].

Another relevant issue is the presence of different pollutant compounds in the fuel stream that can degrade cell materials and highly decrease both reforming and electrochemical performance. As for their amount, the type of such pollutants highly depends on the fuel source and usually consists of sulfur-based and chlorine-based compounds [23]. If biogas is used as fuel, siloxanes may also represent relevant issues [24]. Among these possible compounds, H₂S is the most common poison for catalyst activity. The adsorption of S and the subsequent formation of secondary Ni-S phases can cause serious but reversible degradation of the cell performance [25], making the fuel pre-treatment fundamental. To reduce the possible damage, alternative sulfur tolerant materials, such as metal sulfides [26], are under investigation.

In a view of process improvement, all three configurations of SOFC-reforming integration are analyzed in literature using different levels of details. Such experimental and theoretical studies focus on both new industrial power generation systems based on SOFC technology and on its integration into existing plants. In different works, the simulation is usually performed considering ER or IIR units or, when DIR configurations are presented, simplifying the cell through 0D [27,28] and 1D [22,29] approaches. Although effective in the first feasibility analysis, such models disregard local effects, especially in terms of temperature and current density. This may lead one to contemplate possible solutions that could bring the cell to failure. A more detailed analysis is performed through 2D simulation. Considering a planar geometry, cross-section is commonly assumed as system domain to evaluate the main changes in chemical-physical features in flow direction and along cell thickness [18,30,31]. These studies guarantee a quite good overview of cell behavior in co- and counter-flow configurations, since flow channels can be approximated as plug flow model, whereas in the case of cross-flow design, the analysis should consider occurring gradients on cell plane section that is a less common assumption [32]. The complete knowledge of the system is reached only through a 3D approach [33,34]. This permits a more detailed modeling but requires long computational times to reach a solution penalizing its use. On the other hand, when a tubular cell is simulated, 2D modeling is sufficient to describe completely system behavior [35]. A further step consists of the analysis of occurring degradation and poisoning phenomena. In literature, they have been mostly investigated with experimental tests on both single cells [36] and stacks [37] as well as in terms of regeneration [38]. However, there are few modeling efforts that evaluate coking and H₂S poisoning influence on electrochemical performance based on empirical [39] or theoretical formulation [40].

This work presents a 2D simulation tool developed for DIR-SOFC in industrial applications using a biogas type fuel (mixture of CH₄, CO, CO₂, H_2, and H₂O). To present its features, the authors consider an anode-supported planar cell design in cross-flow configuration. The chemical-physical property changes are described through local balance resolution over the cell plane, whereas material transport phenomena along cell thickness are considered only for the electrochemical kinetics following a previously validated model for H₂-feeding SOFC [41,42]. For methane DIR process, the code considers two different approaches presented in literature: one based on the thermodynamic equilibrium of the reforming reaction and one based on its reaction kinetics. Subsequently, the authors introduce a simplified way to model cell degradation phenomena affecting mainly the reforming reaction. This approach is possible thanks to the introduction of a corrective term that allows reducing the effective active catalysts surface. This coefficient is implemented in matrix form to enable a local distribution and thus identify possible different levels of degradation in each single point of the cell. Following a similar approach, this matrix is also tuned to identify an optimized catalyst distribution in order to guarantee more uniform working temperatures and reactant utilizations. Despite the modeling complexity, the results are obtained with calculations lasting only few minutes.

2. CH₄ Internal Reforming Modeling

2.1. The Cell Reactions

At the cathodic side of an anionic conductive-electrolyte SOFC, O₂ reacts electrochemically to produce O⁼ ions (Equation (1)). These migrate through the ceramic electrolyte to the anodic side where the oxidation of H₂ (Equation (2)) or CO (Equation (3)) completes the circuit with electrical energy production, as shown by the following reaction paths.

$$\frac{1}{2}{\mathrm{O}}_2+2{\mathrm{e}}^{-}\to {\mathrm{O}}^{=} \mathrm{at} \mathrm{the} \mathrm{cathode}$$

(1)

$${\mathrm{H}}_2+{ \mathrm{O}}^{=}\to {\mathrm{H}}_2\mathrm{O} +2{\mathrm{e}}^{-} \mathrm{at} \mathrm{the} \mathrm{anode}$$

(2)

$$\mathrm{CO} +{ \mathrm{O}}^{=}\to {\mathrm{CO}}_2+2{\mathrm{e}}^{-} \mathrm{at} \mathrm{the} \mathrm{anode}$$

(3)

Coupled with Ni catalysts usually employed, the high operating temperature allows for the direct feeding of light hydrocarbons at the SOFC anode. This is possible thanks to internal reforming reactions that convert these hydrocarbons to the H₂ necessary for cell operation.

CH₄ reforming can be carried out using H₂O (steam reforming, SR) or CO₂ (dry reforming, DR) as co-reactants. In literature, both cases are largely investigated, and different reaction paths resulting in the following formulations (Equations (4) and (5) for SR, Equations (6) and (7) for DR) have been proposed.

$${\mathrm{CH}}_4+{ \mathrm{H}}_2\mathrm{O} \leftrightarrow \mathrm{CO} +3{\mathrm{H}}_2$$

(4)

$${\mathrm{CH}}_4+2{\mathrm{H}}_2\mathrm{O}\leftrightarrow {\mathrm{CO}}_2+4{\mathrm{H}}_2$$

(5)

$${\mathrm{CH}}_4+{ \mathrm{CO}}_2\leftrightarrow 2\mathrm{CO} +2{\mathrm{H}}_2$$

(6)

$${\mathrm{CH}}_4+3{\mathrm{CO}}_2\leftrightarrow 4\mathrm{CO} +2{\mathrm{H}}_2\mathrm{O}$$

(7)

Experimental remarks prove that dry reforming reactions (Equations (6) and (7)) develop when the feed has a minimum moisture level and is mainly composed of CH₄ and CO₂. Thus, in SOFC operating conditions, DR can occur only near the anode inlet due to the water generated by electrochemical reactions (Equation (2)). For this reason, in kinetic modeling related to DIR-SOFCs, only the steam reforming route is usually assumed [43].

Due to the applied electric load, parasite electrochemical reactions of CH₄ partial or total oxidation can occur (Equations (8) and (9)). However, reforming is usually faster and thus these last ones are assumed negligible [30].

$${\mathrm{CH}}_4+{ \mathrm{O}}^{=}\to \mathrm{CO} +2{\mathrm{H}}_2+2{\mathrm{e}}^{-}$$

(8)

$${\mathrm{CH}}_4+4{\mathrm{O}}^{=}\to {\mathrm{CO}}_2+2{\mathrm{H}}_2\mathrm{O} +8{\mathrm{e}}^{-}$$

(9)

The presence of all requested reactants and of a suitable catalyst permits the water gas shift (WGS) reaction (Equation (10)) to occur altering the gas phase concentration.

$$\mathrm{CO} +{ \mathrm{H}}_2\mathrm{O}\leftrightarrow {\mathrm{CO}}_2+{ \mathrm{H}}_2$$

(10)

Since Equation (10) is the preferential path for CO consumption [43], only H₂ oxidation is usually considered as anodic electrochemical reaction (Equation (2)). If both the steam reforming (Equation (4)) and the WGS (Equation (10)) are considered reversible, Equation (5) is usually neglected, resulting in the linear combination of these two reactions [44].

In addition to these, carbon deposition mechanisms such as the Boudouard reaction (Equation (11)) can occur due to the high operating temperature (Equations (12)–(16)). Carbon deposition reduces the active area of the catalyst, penalizing reforming reactions and consequently global cell performance.

$$2\mathrm{CO}\leftrightarrow \mathrm{C} +{ \mathrm{CO}}_2$$

(11)

$${\mathrm{CH}}_4\leftrightarrow \mathrm{C} +2{\mathrm{H}}_2$$

(12)

$$\mathrm{CO} +{ \mathrm{H}}_2\leftrightarrow \mathrm{C} +{ \mathrm{H}}_2\mathrm{O}$$

(13)

$${\mathrm{CO}}_2+2{\mathrm{H}}_2\leftrightarrow \mathrm{C} +2{\mathrm{H}}_2\mathrm{O}$$

(14)

$${\mathrm{CH}}_4+2\mathrm{CO} \leftrightarrow 3\mathrm{C} +2{\mathrm{H}}_2\mathrm{O}$$

(15)

$${\mathrm{CH}}_4+{ \mathrm{CO}}_2\leftrightarrow 2\mathrm{C} +2{\mathrm{H}}_2\mathrm{O}$$

(16)

However, working with excess steam hinders Equations (11)–(16) that consequently are usually neglected in the modeling [30]. In addition to carbon deposition, also the presence of poisonous compounds can induce other degradation phenomena. We overlook their effects in the basic model formulation, introducing them in the catalyst deactivation simplified simulation that is discussed afterwards.

Thus, to define a kinetic model that properly describes the SOFC performance, we have to consider interactions between electrochemical (Equations (1) and (2)), SR (Equation (4)), and WGS (Equation (10)) reactions. While electrochemical and WGS processes have been already analyzed in previous studies for the simulation of high temperature fuel cells [45,46], in the following, we focus on possible SR modeling solutions.

2.2. A Model for Steam Reforming Reaction

Literature mainly follows three different modeling approaches [17]:

• equilibrium;
• power law kinetic formulation;
• surface reaction kinetic model.

The equilibrium approach assumes that the reforming reaches the thermodynamic equilibrium (Equation (17), where K_eq,SR is the equilibrium constant, and p_i are the partial pressures of reactants [47].

$$K_{eq,SR}=\frac{p_{CO,eq}{p}_{H_2,eq}^3}{p_{C{H}_4,eq}{p}_{H_{2}O,eq}}$$

(17)

The equilibrium constant can be expressed according to either rigorous Van’t Hoff formulation (Equation (18)) [48] or simplified semi-empirical approach (Equation (19)) [49].

$$\Delta G=-RT\ln K_{eq,SR}$$

(18)

$$K_{eq,SR}={\exp }^{\left(30.114-\frac{26830}{T}\right)}$$

(19)

where ΔG is the Gibbs free energy variation, R the ideal gas constant, and T the temperature. Experimental results show that the CH₄ conversion is lower than the equilibrium one [44]. Thus, the use of this approach may require opportune corrections to properly reproduce experimental data. Since the kinetics are independent from catalyst amount and distribution, they can be used to describe a 0D system.

A power law kinetic formulation approach is based on semi-empirical equations, as described in Equation (20). The SR reaction rate (r_SR) is usually expressed as the product between a kinetic constant (k_SR) and the partial pressure of reactants elevated to different exponents. Both kinetic constant and exponents are fitted through analysis of experimental data.

$$r_{SR}=k_{SR}S{p}_{C{H}_4}^{\alpha }{p}_{H_{2}O}^{\beta }{p}_{H_2}^{\gamma }{p}_{C{O}_2}^{\delta }{p}_{CO}^{\epsilon }$$

(20)

The dependence on H₂, CO_2, and CO is usually negligible (γ, δ, and ε are close to zero), thus the reaction rate is usually assumed as a function of CH₄ and H₂O. Consequently, α and β values may vary among different studies [21,43]. Literature generally agrees that the reforming has a first order dependence on CH₄ partial pressure, while the order of water seems highly influenced by the steam to carbon ratio. It can be positive for low S/C ratio, zero for S/C close to two, and negative for higher values of S/C [50]. This is explained by negative effects that great amount of water has on the CH₄ adsorption on the catalyst surface [44]. This approach does not require the knowledge of the mechanisms involved and can be easily applied when a large number of experimental data is available for model tuning. However, results are specific to the analyzed case and cannot be applied to different systems.

The surface reaction kinetic model approach describes occurring mechanisms as a sequence of intermediate phenomena, consisting of adsorption, surface reaction, and desorption of all present gases. The rate of the total kinetics is determined by the slowest phenomenon that changes at different temperatures and reactant-product compositions. Such kinetics are usually modeled following the Langmuir–Hinshelwood (LH) or the Hougen–Watson (HW) approaches. The first one assumes a bimolecular reaction between two reactants adsorbed on neighboring sites as the rate-limiting step and the water dissociation into atomic hydrogen H and hydroxyl groups OH (Equation (21)) [51].

$$r_{SR}=\frac{k_{SR}S\prod p_i^{{\phi }_i}}{{\left(1+\sum K_i{p}_i^{{\phi }_i}\right)}^{\chi }}\left(1-\frac{Q_{SR}}{K_{eq,SR}}\right)$$

(21)

The second approach also takes into account the sorption and the reaction of intermediates (Equation (22)) [51].

$$r_{SR}=\frac{k_{SR}S\prod \frac{p_i^{{\phi }_i}}{p_j^{{\lambda }_j}}}{{\left(1+\sum K_i\frac{p_i^{{\phi }_i}}{p_j^{{\lambda }_j}}\right)}^{\chi }}\left(1-\frac{Q_{SR}}{K_{eq,SR}}\right)$$

(22)

In both Equations (21) and (22), the numerator shows the kinetics dependency on involved gases, while the denominator considers the availability of active sites through adsorption isotherm. The last term, expressed as ratio between the reaction quotient Q_SR (Equation (23)) and the equilibrium constant K_eq,SR, represents the driving force of the overall process.

$$Q_{SR}=\frac{p_{CO}{p}_{H_2}^3}{p_{C{H}_4}{p}_{H_{2}O}}$$

(23)

Both kinetic k_SR and adsorption K_i coefficients can be described by an Arrhenius type dependency on the operating temperature (Equations (24) and (25)).

$$k_{SR}=k_0{\exp }^{-\frac{E_{act, SR}}{RT}}$$

(24)

$$K_i=K_{0,i}{\exp }^{- \frac{\Delta H_{ads, i}}{RT}}$$

(25)

where k₀ and K_0,i are the pre-exponential coefficients, E_act the activation energy of SR reaction, and ΔH_ads the adsorption enthalpy variation. In DIR-FC modeling, Equation (22) is commonly used, considering the values detected for Ni-MgAl₂O₃-spinel catalysts as reference of kinetics parameters [49]. Since the process rate is strongly influenced on catalyst features, such as used support, Ni percentage, and particle size, this approach is not always effective for SR occurring inside a fuel cell due to higher Ni content compared to the traditional SR catalyst needed to guarantee a good conductivity [21].

Experimental data suggest that, for Ni/YSZ, the material usually employed as the SOFC anode, the rate-limiting step is the CH₄ dissociative adsorption. Thus, a first order expression function of only CH₄ partial pressure is formulated in accordance with power law models. Under the assumption that the surface cannot be covered by other components, the adsorption dependency is neglected, and the kinetic rate is expressed through Equation (26) as reported in different works [21,44,51].

$$r_{SR}=k_{SR}S{p}_{C{H}_4}\left(1-\frac{Q_{SR}}{K_{eq,SR}}\right)$$

(26)

The kinetics constant k_SR usually depends on available catalyst active area. This kinetics has been validated over a wide range of temperatures and S/C ratios.

For the DIR-SOFC simulation, we implement in our code both equilibrium (Equation (17)) and surface reaction kinetics approach (Equation (26)). It is important to underline that the first is independent of the specific used catalyst, while the second one is expressed in function of its distribution on cell plane. This study takes as reference the kinetics proposed by [52] for a nickel cermet electrode (

Table 1. Reference kinetics coefficients [52].

Kinetics Parameter	Value
kSR [mol m−2 s−1 bar−1]	4274
Eatt,SR [J mol−1]	82,000

).

3. The DIR-SOFC simulation

To simulate a DIR-SOFC, the authors integrate the presented reforming kinetics into a home-made Fortran written code called SIMFC (SIMulation of Fuel Cell) that has been developed and validated in previous works. It is a 2D model that simulates the operation of high temperature molten carbonate [46,53] and solid oxide [41,42] fuel cells on the basis of local mass, energy, charge and momentum balances.

The SOFC electrochemical performance is evaluated through a semi-empirical formulation, derived by detailed experimental campaigns on H₂-feeding cell [42,45], where the main dependences of different polarization terms have been investigated. The equations proposed to describe cell voltage and used parameters are presented in

Table 2. Solid oxide fuel cell (SOFC) electrochemical kinetics (for the complete description refer to [45]).

Cell voltage	$\begin{array}{cc}\hfill V=E-{\eta }_{leakage}& -P_{1}T{\exp }^{\frac{P_2}{T}}J-\frac{RT}{F}{\sin \mathrm{h}}^{-1}\left(\frac{J}{2J_{0,an}}\right)\hfill \\ -& \frac{RT}{2F}{\sin \mathrm{h}}^{-1}\left(\frac{J}{2J_{0,cat}}\right)\hfill \\ -& \frac{RT}{2F}\ln \left(\frac{{\left(1+\frac{RT{d}_{an}J}{6F{D}_{H_{2}O}^{eff}{p}_{H_{2}O,an}}\right)}^{2B}}{{\left(1-\frac{RT{d}_{an}J}{6F{D}_{H_2}^{eff}{p}_{H_{2,an}}}\right)}^{2A}}\right)\hfill \end{array}$
OCV potential	$E=E^0+\frac{RT}{2F}\ln \left(\frac{p_{H_2,an}{p}_{O_2,cat}^{0.5}}{p_{H_{2}O,an}}\right)$
Anodic exchange current density	$J_{0,an}=P_3{\left(y_{H_{2,an}}\right)}^A{\left(y_{H_{2}O,an}\right)}^B{\exp }^{- \frac{E_{act,an}}{RT}}$
Cathodic exchange current density	$J_{0,cat}=P_4{\left(y_{O_2,cat}\right)}^C{\exp }^{- \frac{E_{act,cat}}{RT}}$
Electrochemical parameters	ηleakage = 0.03, P1 = 2 × 10−9 [Ω cm2 K−1], P2 = 10,986 [K], P3 = 2.8 × 105 [A cm−2], P4 = 4 × 106 [A cm−2], Eact,an = 110 [kJ mol−1], Eact,cat = 120 [kJ mol−1], A = 0.5, B = 0.55, C = 0.25

. The same previous coefficients identified for H₂-feeding systems are valid also for DIR configuration, since the reactive sites for two reaction paths do not coincide [21].

The proposed model can be used not only to obtain overall cell information, such as cell voltage or produced power, but also to evaluate local variables, such as compositions, polarizations, and reactant utilization factors. This feature is extremely important for both theoretical and feasibility studies. Local data are fundamental for the understanding of the phenomena occurring inside SOFC as well as to avoid critical points and guarantee a more uniform utilization of the cell surface. As many properties cannot be directly measured, a detailed modeling is necessary to assess their values. Moreover, since the code needs only few minutes of computations to produce the results, it could be integrated successfully in more complex software for plant simulation (e.g., Aspen Plus).

4. Results and Discussion

The simulation tool is presented considering a single cell anode-supported SOFC in cross-flow configuration. The cell has a planar geometry with a total active area equal to 1 m² (anode inlet length = 71 cm, cathode inlet length = 142 cm). This design is commonly used in industrial applications due to both lower fabrication costs and a relatively more compact structure compared to tubular design [54]. In view of common used materials [55], it is assumed a Ni/YSZ porous anode provides the base for the thin dense YSZ electrolyte and the porous LSC cathode. The main cell features are presented in

Table 3. Physical and microstructure properties of different cell layers.

Property	Anode Ni/YSZ	Electrolyte YSZ	Cathode LSC
Density [g cm−3]	7.7	6.0	5.3
Heat Capacity [cal mol−1 K−1]	50	29	34
Porosity [-]	0.40	0.01	0.35
Thickness [μm]	350	5	30
Tortuosity [-]	4	-	n.a.

.

The SOFC performance has been studied setting up a current density load of 0.15 A cm⁻² coupled with an H₂ utilization of about 74%, common target value according to numerous SOFC producers [56]. The proposed feed conditions, with inlet S/C ratio of about two (safety condition to avoid carbon deposition at the cell operating temperature [57]), are represented in

Table 4. Inlet operating conditions.

Inlet Condition	Anode	Cathode
T [K]	1023	1023
N [Nm3 h−1]	0.75	8.60
yCH4 [-]	0.25	-
yH2 [-]	0.04	-
yH2O [-]	0.51	-
yCO2 [-]	0.13	-
yCO [-]	0.03	-
yN2 [-]	0.04	0.79
yO2 [-]	-	0.21

.

4.1. Analysis of DIR-SOFC Operation at Local Level

The comparison between DIR-SOFC simulation assuming reforming reaction at the equilibrium (Figure 1A,C,E) and DIR-SOFC simulation assuming reforming reaction as kinetics driven (Figure 1B,D,F) is presented in term of local maps for the H₂ molar fraction (Figure 1A,B), the temperature of the solid structure (Figure 1C,D), and the applied current density (Figure 1E,F). In both cases, CH₄ is almost completely consumed in proximity of the anode inlet (about 10 cm), thus its mapping on cell plane is not reported here.

As verified before, the high operating temperature allows for a rapid conversion of CH₄. This is confirmed in literature, where CH₄ is observed consuming in the so called “reforming zone” [18]. The equilibrium conversion represents the maximum value that can be theoretically achieved. This is also evident in the kinetics reaction rate expression (Equation (26)), where the imbalance between actual and equilibrium composition represents the driving force of the process (Q_SR/K_eq,SR). The consequence is that the equilibrium case foresees a slightly faster SR to form H₂ and CO, reducing local temperature. However, in both approaches, the newly formed H₂ reacts electrochemically to produce water at the anode (Equation (2)), while CO mainly produces additional H₂ via WGS (Equation (10)). These two reactions are exothermic and balance the local temperature decrease due to the reforming. As result, CH₄ is almost all converted in proximity of the anode inlet. H₂ is immediately formed thanks to the reforming and rapidly depletes along the cell plane to sustain electrochemical reactions (Figure 1A,B). As expected, the temperature (Figure 1C,D) decreases at the inlet due to the reforming endothermicity and increases moving towards the anode and, to a lesser degree, the cathode outlet with a peak (Figure 1E,F) where both exothermic electrochemical and WGS reactions occur [58]. Due to the cross-flow configuration, inlet cathodic gas (bottom left corner of maps) reduces the anode inlet temperature proceeding toward the cathode outlet. This in turn penalizes the reforming and the subsequent anodic reactions, forbidding the temperature increase. The assumed kinetics have a relevant influence on electrochemical processes. If the fastest equilibrium reforming kinetics induce an initial peak of H₂ production followed by a slow electrochemical conversion (Figure 1A), in the surface kinetics formulation a wider H₂ conversion zone (Figure 1B) is detected, causing lower peaks of temperature and local current density (Figure 1D,F). However, as shown in

Table 5. Main results of different simulations basing on equilibrium and surface reaction mechanism approaches.

Reforming Operating Condition	Equilibrium	Kinetics
V [V]	0.747	0.856
Jmax [A cm−2]	0.88	0.42
Tmax [K]	1110	1084
Tmin [K]	873	828
Taverage [K]	1037	1034
yCH4,max [-]	0.048	0.196
yH2,max [-]	0.402	0.525
Cell Power [W]	1121	1284

, the macroscopic results of the simulation are not particularly dissimilar, thus both approaches can be used as preliminary analysis.

However, these obtained results are not completely satisfactory. As previously mentioned, the rapid conversion of CH₄ near the anode inlet induces a relevant temperature gradient with a difference between maximum and minimum temperature of more than 200 K in the equilibrium case and more than 250 K in the kinetics case. This, coupled with the high peak current density, greatly speeds up degradation processes limiting the cell operability.

For this reason, the catalyst active area is usually reduced by a partial poisoning or by changing the material microstructure [21].

4.2. Possible Applications of Local Simulation Tool

As already underlined, the local simulation has several advantages: it allows for not only a detailed knowledge of the main chemical-physical features on the cell plane but also a local description of system structure. Thus, using a local modeling approach, we can study the degradation of the reforming catalyst and how it affects the cell performance.

To simulate the catalyst deactivation that happens due to several phenomena such as sintering, poisoning by sulfur and other pollutants, or carbon deposition, the authors adjust the reforming kinetics introducing a coefficient (σ) as shown in Equation (27). This coefficient represents a corrective parameter that permits one to consider a reduced active area or an unevenly distributed catalyst. Since the equilibrium kinetics are independent of the catalysts active area; this point is added only to the surface reaction approach.

$$r_{SR}=k_{SR}\sigma S{p}_{C{H}_4}\left(1-\frac{Q_{SR}}{K_{eq,SR}}\right)$$

(27)

An uneven distribution of SR active sites due to advanced degradation processes can be introduced considering local value of σ to differentiate the cell area. This is specifically relevant when there is the need to simulate long-term applications when such deactivation phenomena cannot be neglected any more. To demonstrate this possibility, the authors introduce into the code a matrix of σ, as presented in Figure 2. The cell surface is divided into a 20 × 20 mesh with each sub-cell having an area equal to 25 cm². The start of degradation processes is assumed at the anode inlet where the sigma value is the smallest (i.e., fewer catalyst active areas). Moving towards the anode outlet, σ increases, reflecting a milder degradation effect (i.e., more catalyst active areas). Such distribution may derive from carbon deposition (Equations (11)–(16)) that usually begins at the anode inlet, where reactions initially take place and the temperature is lower, and then it sequentially shifts to the outlet as more and more active sites are covered. Another possible cause could be the long exposition to poisoning of sulfur or other poisonous compounds present in biogas feed. Indeed, it is well known that sulfur can react with the reforming catalyst inhibiting its properties [25]. Additionally, in this case, the poisoning starts at the anode inlet and, as more sites are poisoned, it shifts towards the outlet. The local results of the simulation are reported in Figure 3.

Close to the anode inlet, the CH₄ conversion is negligible due to the highly reduced catalytic surface area that corresponds to a value of 0.00001 in the σ matrix. This can be interpreted as an extremely deactivated or poisoned catalyst. As in Figure 3A, then moving towards the anode outlet, the conversion begins to increase when the fuel flow encounters the more active catalyst (σ from 0.01 to 0.7 progressively). The maps of H₂, current density, and temperature on the cell plane are also influenced, resulting in a more uniform distribution (Figure 3B–D). As

Table 6. Main results of the simulation using a local map of σ to describe catalyst deactivation.

Reforming Operating Condition	Degradation
V [V]	0.823
Jmax [A cm−2]	0.47
Tmax [K]	1088
Tmin [K]	962
Taverage [K]	1034
yCH4,max [-]	0.250
yH2,max [-]	0.373
Cell Power [W]	1234

shows, there is not an excessive global performance variation compared to the kinetics case previously described in Table 5; still, an inefficient use of the catalyst occurs.

It is important to underline that the decreased temperature gradient (less than 150 K between maximum and minimum values) should not be interpreted as positive outcomes. It is the result of a highly deactivated catalyst that penalizes the performance in terms of power output decreasing by about 50 W in comparison with the kinetics case. In this simulation, the authors do not consider the degradation effects on electrochemical reactions, which could be added with available experimental data and would more highly affect power losses.

Through a similar local approach, it is also possible to identify an optimized catalyst distribution that should yield a reduced solid temperature variation compared to the kinetics base one shown in Figure 1D. Having a lower solid temperature gradient allows for the improvement of cell stability and consequently for a longer cell durability by reducing the mechanical stress on materials. An uneven catalyst distribution can be obtained in the manufacturing process through several ways, such as different Ni particle size, addition of alkali metals in lattice, or Cu doping [21]. To simplify the optimization process, the authors decide to divide the cell surface in rectangular shaped areas of similar dimension and catalyst amount. Then, using iterative calculations, the σ matrix is tuned in order to reduce the temperature gradient, resulting in the configuration of Figure 4.

Using this σ matrix, the cell performance is reevaluated, as presented in Figure 5 and Figure 6. A more uniform temperature distribution is obtained on the cell plane, which maintains good global performance (

Table 7. Main results of the simulation using local maps of σ to optimize cell performance.

Reforming Operating Condition	Optimized Configuration
V [V]	0.847
Jmax [A cm−2]	0.29
Tmax [K]	1059
Tmin [K]	963
Taverage [K]	1034
yCH4,max [-]	0.233
yH2,max [-]	0.322
Cell Power [W]	1270

). For an easy comparison, the temperature profiles of the solid structure from anode inlet to anode outlet are presented for kinetics base case (Figure 6A,C,E) and optimized approach (Figure 6B,D,F) in different cell positions in correspondence with the cathode inlet (Figure 6A,B), in the middle between cathode inlet and outlet (Figure 6C,D), and in correspondence with the cathode outlet (Figure 6E,F). As can be seen, the temperature change in optimized catalyst distribution is much less sharpened in all analyzed points. Moreover, the difference between the lower and the higher temperature is less than 100 K, whereas it is more than 250 K in the kinetics base case. Since it is estimated that about 30% of occurring stresses are due to thermal gradients, this configuration offers a relevant improvement to guarantee long-term operations [59]. This result is obtained by expanding the CH₄ reforming zone in order to have a more homogeneous conversion and consequently a minor temperature drop caused by the reforming endothermicity. Compared with previous solution, the current density and the H₂ molar fraction also have lower maximum values (J_max = 0.29 A cm⁻² vs. 0.42 A cm² in kinetics base case, y_H2,max = 0.322 vs. 0.525 in kinetics base case), signifying less stressed working conditions. However, the measured voltage and the total output power are slightly decreased in this configuration; reductions of 9 mV and 14 W are obtained.

5. Conclusions

This work aims at presenting a suitable tool for the simulation of DIR-SOFCs for industrial applications, with regard to local effects on the main chemical-physical variables (e.g., composition of reactants, current density, and solid temperature) that are fundamental for both feasibility studies and to favor long-term operations.

The SIMFC code, a home-made tool previously validated for H₂-feeding high temperature fuel cells, is improved, introducing the possibility to consider internal CH₄ reforming. Based on experimental observations and literature modeling, both thermodynamic equilibrium and surface reaction kinetic approaches are implemented to simulate the direct internal reforming using a biogas as feeding. In both cases, similar results are obtained; the CH₄ conversion occurs in the initial 10 cm, and the obtained voltage is around 0.8 V.

Being independent of the catalysts used, the equilibrium approach does not permit to differentiate catalyst operation on the cell plane. However, this is possible using the kinetic approach that is consequently more suitable for a detailed local modeling. As examples of possible applications, the authors analyze (i) how to simulate SOFC performance affected by the deactivation of reforming catalysts and (ii) how to use also local control in order to detect an optimized catalyst distribution that would provide minimized thermal stresses on the cell. This is achieved with the introduction of a specific coefficient (σ) in matrix form to describe the catalyst amount on the cell surface.

When catalyst deactivation is considered, the term σ indicates the degree of degradation or poisoning of reforming active sites, with lower values corresponding to reduced catalytic activity. At such conditions, the reforming reaction is deeply penalized, reaching the complete conversion over half along cathode inlet direction in the presented example.

When used to determine an optimized catalysts distribution, the value of σ indicates the quantity of catalyst or, in other terms, the active surface available for the reactions. In order to reduce the solid temperature gradient on the cell plane, in the example, the authors tune the matrix of σ coefficients to obtain an optimized distribution that can be easily realized in manufacturing process. A steep temperature profile induces higher mechanical stresses that can result in creep formation and ultimately global anode degradation. Thanks to the identified simple distribution, it is possible to sensitively reduce the temperature gradient by about 60% compared to the kinetics base example. The proposed configuration also permits one to halve the peak of current density. At the same time, the complete CH₄ conversion is guaranteed without a relevant power penalization that is reduced by only 14 W.