Numerical Simulation of Processes in an Electrochemical Cell Using COMSOL Multiphysics

Abstract

Fuel cells are a promising source of clean energy. To find optimal parameters for their operation, modeling is necessary, which is quite difficult to implement taking into account all the significant effects occurring in them. We aim to develop a previously unrealized model in COMSOL Multiphysics that, on one hand, will consider the influence of electrochemical heating and non-isothermal fluid flow on the temperature field and reaction rates, and on the other hand, will demonstrate the operating mode of the Solid Oxide Fuel Cell (SOFC) on carbonaceous fuel. This model incorporates a range of physical phenomena, including electron and ion transport, gas species diffusion, electrochemical reactions, and heat transfer, to simulate the performance of the SOFC. The findings provide a detailed view of reactant concentration, temperature, and current distribution, enabling the calculation of power output. The developed model was compared with a 1-kW industrial prototype operating on hydrogen and showed good agreement in the volt-ampere characteristic with a deviation not exceeding 5% for the majority of the operating range. The fuel cell exhibits enhanced performance on hydrogen, generating 1340 W/m² with a current density of 0.25 A/cm². When fueled by methane, it produces 1200 W/m² at the same current density. Using synthesis gas, it reaches its peak power of 1340 W/m² at a current density of 0.3 A/cm².

1. Introduction

Fuel cells are fundamentally made up of an electrolyte layer flanked by two porous electrodes. The central role of the electrolyte layer, often considered the core of the fuel cell, is to determine the operational temperature range of the cell, the nature of electrochemical reactions occurring at the electrodes, and the kind of catalyst materials within the electrodes.

In SOFCs, electrochemical reactions occur at the triple phase boundary (TPB) where the electrolyte, electrode, and gas phase meet [1,2].

At the anode, hydrogen gas is oxidized, releasing electrons into the circuit and protons into the electrolyte. At the cathode, oxygen from the air accepts these electrons and reacts with the protons to form water. This electrochemical process is accompanied by transport phenomena, including electronic and ionic conduction, gas diffusion, and heat transfer.

Various configurations of planar SOFC designs exist, including electrode-backed SOFCs, SOFCs with electrolyte support, and symmetrical SOFCs [3,4]. Symmetrical SOFCs enhance the harmony between cell parts and streamline the cell design. Anode-backed SOFCs can notably diminish resistive losses, leading to a decrease in the SOFC’s operational temperature, and a more robust anode can facilitate internal reforming [5].

Our goal is to create a model for different types of fuel, providing data on species composition, temperature, and power, while considering the heat transfer due to the movement of the working gas. Modeling data is crucial for understanding the operation of the fuel cell and optimizing its performance.

COMSOL Multiphysics is a versatile tool that offers advanced capabilities for modeling various physical processes. It is efficient and effective in numerical simulation in several fields such as: porous media, desalination, renewable energy, chemical engineering, and heat exchangers [6,7,8,9,10,11] making it particularly suitable for simulating fuel cell dynamics.

COMSOL is widely used for SOFC simulations. However, in most of these, either the temperature distribution was not accounted for, the effects of non-isothermal flow and electrochemical heating were not considered, or the selection of the exchange current was not disclosed [3,4,12,13,14].

In our study, we aim to address these issues. The developed model is compared with a real, industrial prototype running on hydrogen. To the best of our knowledge, there have been no similar models implemented in COMSOL Multiphysics according to the literature. However, there are similar models for methane fuel cells, such as the research conducted by Cai, W [5]. Nevertheless, in this study, a different fuel composition was considered, and the temperature regime, power, and configuration of the fuel cell also differed.

2. Materials and Methods

The 3D geometric representation was crafted based on the actual stack design’s specifications. This model takes inspiration from a 1 kW planar SOFC with anode support, originating from China. The SOFC’s dimensions are 16 × 16 cm², with a designated active region measuring 10 × 10 cm². Gas pathways for entry and exit are facilitated through specific inlets and outlets in the stack. Every module of the fuel cell encompasses a membrane-electrode assembly, which includes the cathode (or the positive electrode), the electrolyte, the anode (or the negative electrode), channels for both air and fuel and interconnecting structures. Consequently, the electrochemically active region for a stack of 30 SOFC elements is composed of 900 uniform module blocks. Given their identical nature, a computational simulation was executed for a single SOFC cell block. The block’s specific measurements can be found in

Table 1. SOFC channel geometric parameters.

SOFC Geometric Parameter	Value
Channel width	1.4 mm
Rib width	2 mm
Anode thickness	550 μm
Cathode thickness	30 μm
Electrolyte thickness	15 μm
Gas flow channel height	1.6 μm

.

The model being developed is schematically shown in Figure 1. Here, the computational domains are correlated with the set of physics being solved, and materials and boundary conditions are shown. The sets of equations corresponding to each physics will be discussed further.

The conductivity of the electrolyte is primarily represented by ionic conductivity. Given that the electronic conductivity is extremely minimal, we will not consider it in this simulation. Typically, ionic conductivity falls within the range of 1 to 10 S/m and is temperature-dependent. Additionally, there exists a model for estimating conductivity based on activation energy:

$${\sigma }_{el}={\sigma }_o{T}^{-1}\mathrm{e}\mathrm{x}\mathrm{p}(-\frac{E_{el}}{RT})$$

(1)

COMSOL already incorporates a temperature-dependent conductivity model for the 8YSZ electrolyte material. This dependency is established through the interpolation of tabular data. For instance, the conductivity is 5.1 S/m at 800 °C, 7.44 S/m at 850 °C, and 10.66 S/m at 900 °C.

Simultaneously, an experiment cited in [15] found a value quite close to this, at 4.2 S/m, at a temperature of 800 °C. Moving forward, we will use the data embedded in COMSOL for our model.

High conductivity is crucial for an anode. In the case of Ni/YSZ material, Ni provides an effective electronic conductivity of about 100,000 S/m, [16,17], while YSZ contributes an ionic conductivity of about 1 S/m [18].

As per the manufacturer’s data, our anode exhibits an electronic conductivity ${\sigma }_a=333330[\mathrm{S}/\mathrm{m}]$.

LSCF + GDC cathode enables both kinds of conduction. The ionic conductivity is approximately 5 S/m, and the electronic conductivity is also significantly higher.

According to the manufacturer of our sample, our cathode displays an electronic conductivity ${\sigma }_c$ of 7937 [S/m].

2.1. Electrochemistry

2.1.1. Electromotive Force

The following reactions occur in a fuel cell when hydrogen is used as a fuel. The inclusion of additional components such as methane and syngas is elaborated upon in Section 2.5.

$$\begin{gathered} Anode:H_2\to 2H^{+}+2e^{-} \\ Cathode:\frac{1}{2}{O}_2+2e^{-}\to O^{2-} \\ Anode:O^{2-}+2H^{+}\to H_{2}O \\ Overall:H_2+\frac{1}{2}{O}_2\to H_{2}O \end{gathered}$$

(2)

The electromotive force of the aforementioned electrochemical system is depicted by the Nernst equation [19]:

$$E_{eq}=E_{eq,ref}-\frac{RT}{nF}ln\frac{{\mathsf{\Pi }(x}_{pt}^v)}{{\mathsf{\Pi }(x}_{rt}^v)}$$

(3)

$E_{eq}$—the reversible voltage or Equilibrium potential in terms of COMSOL; $x_{pt}^v$ and $x_{rt}^v$—the concentration of the reactants and products at the reaction sites, and the superscript v stands for the stoichiometric coefficient

$$E_{eq,ref}=\frac{-∆G_r^0}{nF}$$

(4)

where $E_{eq,ref}$ represents the standard electrochemical cell voltage [20] (also known as the reference equilibrium potential in COMSOL terms). It is tabulated for most common reactions under standard conditions (∆G0r).

For instance, [19] provides the following equation for the reference equilibrium potential:

$$E_{eq,ref}=1.271-{2.7311}^{-4}T$$

(5)

COMSOL Multiphysics employs the same approach, but in this case, $E_{eq,ref}$ is determined based on changes in enthalpy and entropy, using data from NIST [21].

2.1.2. Electrode Kinetics

The Butler-Volmer equation, in the form proposed by Kawada et al. [22], is extensively used in the numerical simulation of SOFCs, as evidenced by references [23,24,25,26,27]. The equation accounts for the rate of an electrochemical reaction by considering both forward and reverse reaction rates.

This is frequently employed in a form that associates current density per unit volume with the TPB (termed as ‘Specific Surface Area’ in COMSOL).

The equation is then given by:

$$i=l_{TPB}{i}_0[\exp \left(\frac{{\alpha }_{a}nF}{RT}\right){\eta }_{act}-\exp \left(-\frac{{\alpha }_{c}nF}{RT}\right){\eta }_{act}]$$

(6)

where $i$—the volumetric exchange current density, A/m³; $i_0$—is the equilibrium exchange current density per unit volume associated with the TPB reaction in the anode (A/m³) when $E_{}=E_{eq}$ and it is estimated based on the fitting to the experimental data [28]; ${\eta }_{act}$ is the activation overpotential, V; n is the number of electrons involved in the electrode reaction; ${\alpha }_a$ and ${\alpha }_c$ are the charge transfer coefficients that are obtained by fitting to the experimental data [29], and F is the Faraday constant.

The activation overpotential is defined as:

$${\eta }_{act}={{\phi }_s-{\phi }_l-E}_{eq}$$

(7)

${\varphi }_s$—denotes the electric potential at the surface; ${\varphi }_l$—is the electrolyte potential; $E_{eq}$—denotes the equilibrium electrode potential; $E_{th}={\varphi }_s-{\varphi }_l$

To utilize the Butler-Volmer equation, one needs to determine the equilibrium exchange current density, charge transfer coefficients, and TPB length.

2.1.3. Equilibrium Exchange Current Density

In the anodes of SOFCs, electrochemical oxidation of hydrogen occurs at the TPB. Currently, the estimation of the charge-transfer rate in SOFC electrodes remains a subject of active research, with numerous studies devoted to investigating the charge-transfer rate in SOFC electrodes. T. Yonekura [21] consolidated the calculations of many other authors, presenting empirical formulas for the exchange current density at the anode and cathode, while also reviewing the existing coefficients in the specified empirical formula. Kota Mioshi [30] demonstrates a similar approach in his studies of porous LSM cathodes.

$$\begin{gathered} i_{0,c}={\gamma }_c{\left(\frac{P_{O_2}}{P_{O_{2,ref}}}\right)}^C\exp \left(\frac{-E_{act,c}}{RT}\right) \\ {i}_{0,a}={\gamma }_a{\left(\frac{P_{H_2}}{P_{H_{2,ref}}}\right)}^A{\left(\frac{P_{H_{2}O}}{P_{H_{2}O,ref}}\right)}^B\exp \left(\frac{-E_{act,a}}{RT}\right) \end{gathered}$$

(8)

While determining the equilibrium potential using the Nernst Equation, the equilibrium exchange current density i0’s concentration reliance can be consistently established thermodynamically. This is in line with the Nernst equation, combined with a reference exchange current density i_0,ref (A/m²), representing the exchange current density when Eeq matches E_eq = E_eq,ref.

So equilibrium exchange current density:

$$i_o=i_{o,ref}(T)\prod _{i:v_i>0}{\left(\frac{p_i}{p_{i,ref}}\right)}^{{\alpha }_c{v}_i/n}\prod _{i:v_i<0}{\left(\frac{p_i}{p_{i,ref}}\right)}^{{-\alpha }_a{v}_i/n}$$

(9)

However, we still need to specify $i_{o,ref}$ for anode and cathode. Based on the empirical formula [30]:

$$i_{o,ref}=\gamma \mathrm{e}\mathrm{x}\mathrm{p}(-\frac{E_{act}}{RT})$$

(10)

According to (10) $i_{o,ref\_c}=0.11\mathrm{A}/{\mathrm{m}}^2$

Unfortunately, a literature review has shown that the values can vary greatly, so we will choose this value based on our own experimental data.

The most commonly used form of the Butler-Volmer formula, given by Equation (6), assumes charge transfer coefficients ${\alpha }_a=0.5$ and ${\alpha }_c={1-\alpha }_a=0.5$ [31].

2.1.4. Specific Surface Area

This parameter designates the surface area where the catalytic reaction will occur. It can be derived from the outcomes of 3D-FIB modeling. However, it is important to highlight that not all of the surface area will necessarily partake in the reaction.

As per Vivet et al. [32], who conducted a FIB-SEM procedure to obtain a high-quality 3D structure of the electrode, a sample volume of 8.66 × 9.79 × 11.41 μm³ was reconstructed. This unveiled the presence of 99.8%, 99.1%, and 87.4% percolation of the pore, YSZ, and Nickel in the structure, respectively. The specific surface areas of the Pore, YSZ, and Nickel phases were observed as 4.27, 4.24, and 2.33 μm²/μm³, respectively. Approximately 50% of the Nickel surface area was assigned as a pore region which could potentially be utilized for surface catalytic reactions. So, for the anode $l_{TBP,a}$ = 1.165 μm²/μm².

2.2. Heating

In an electrochemical cell, certain factors can lead to irreversible voltage losses [33], including:

• Electrochemical heating
• Charge transport in the electrolyte (Joule heating)
• Charge transport in the solid conductor materials (Joule heating)
• Heat of mixing

2.2.1. Electrochemical Heating

Heat is produced as a result of the electrochemical reactions taking place at the active sites where the electrolyte meets the electrodes [34].

In order to calculate the reaction heat source, we can use thermoneutral voltage

$$E_{therm}=\frac{-∆H_r^0}{nF}$$

(11)

Thermodynamic formulas for gaseous species of $H_i\left(T\right)$ at a standard partial pressure of 1 atm are derived from NIST datasets [21].

$$q=i(E_{eq}-E_{th})$$

(12)

2.2.2. Joule Heating Due to Charge Transport

When charged particles move within an electric field, electrical energy transforms into thermal energy. The terms representing Joule heating in both the electrode and electrolyte phases are:

$$q=\frac{i^2}{\sigma }$$

(13)

Obviously, it is necessary to take into account both the ionic and electronic components of the current in the calculation. The electronic and ionic conductivities will also vary for different sections of the fuel cell.

2.2.3. Heat Generation from Mixing Effects

When the enthalpy is influenced by the local concentration of the species involved in the reaction, it is essential to consider the heat sources related to concentration gradients. These gradients lead to the molecular flux of the reacting species moving from the bulk to the surface, ensuring the cell’s thermal equilibrium. While the effects of heat from mixing are usually minimal (non-existent for ideal gases), they are often omitted in Electrochemistry interfaces.

2.3. Fluid Dynamic Model

This section describes mass transport in porous media.

We consider the fluid flow in the channel to be laminar, as for the given geometry and flow velocity, the Reynolds number will be significantly less than 2000. For instance, for methane flow in the channel at a temperature of 900 °C, the Reynolds number will be 3.6.

In the gas phase nodes, calculations are made for a series of mass fraction variables, denoted as ω_i, with i being the species index. The primary equations used for these calculations stem from the Maxwell–Stefan equations combined with Darcy’s Law.

2.3.1. Maxwell–Stefan Description

Consider a reactive flow composed of a blend with species i ranging from 1 to Q and reactions j from 1 to N. The subsequent section delineates the mass transport for a specific species:

$$\frac{\partial }{\partial t}\left(\rho {\omega }_i\right)+\nabla ·\left(\rho {\omega }_{i}u\right)=-\nabla ·{\mathrm{j}}_i+R_i$$

(14)

where ρ—denotes the mixture density, kg/m³, u the mass averaged velocity of the mixture, m/s; ${\omega }_i$—the mass fraction; ${\mathrm{j}}_i$—the mass flux relative to the mass averaged velocity, kg/(m²·s)]; $R_i$—the rate expression describing its production or consumption, kg/(m³·s).

For a mixture with multiple components, the mass flux, when compared to the average mass velocity, j_i, can be characterized using the generalized Fick’s equations. [35]:

$${\mathrm{j}}_i=-\rho {\omega }_i\sum _{k=1}^Q{\tilde{D}}_{ik}{d}_k-D_i^T\nabla lnT$$

(15)

where ${\tilde{D}}_{ik}$ are the multicomponent Fick diffusivities, m²/s; $D_i^T$ are the thermal diffusion coefficients, kg/m·s; $d_k$ is the diffusional driving force acting on species k, 1/m.

Utilizing the Maxwell–Stefan diffusion approach, the transport equations related to the mass of the species are as follows:

$$\rho \frac{\partial }{\partial t}\left({\omega }_i\right)+\rho \left(u·\nabla \right){\omega }_i=\nabla \left(\rho {\omega }_i\sum _{k=1}^Q{\tilde{D}}_{ik}{d}_k+D_i^T\frac{\nabla T}{T}\right)+R_i$$

(16)

In the case of ideal gas mixtures, the force driving diffusion is: [36]:

$$d_k=\nabla x_k+\frac{1}{p}\left[{(x}_k-{\omega }_k)\nabla \mathrm{p}-\rho {\omega }_k{g}_k+{\omega }_k\sum _{l=1}^Q\rho {\omega }_l{g}_l\right]$$

(17)

where $x_k$ mole fraction; p—is the total pressure, Pa; $g_k$—an external force (per unit mass) acting on species k, m/s².

For ionic species, the external force is generated by the presence of an electric field.

2.3.2. Binary Diffusion Coefficients

The determination of the inherent binary diffusion coefficients is carried out utilizing the Fuller–Schettler–Giddings (FGS) approach [36]:

$$D_{i,j}^0=\frac{1.01325·{10}^{-2}{T}^{1.75}{\left(\frac{1}{M_i}+\frac{1}{M_j}\right)}^{1/2}}{P\left({\left(\sum _{i}v\right)}^{\frac{1}{3}}+{\left(\sum _{j}v\right)}^{\frac{1}{3}}\right)}$$

(18)

where T denotes the temperature, K; M_i the molecular weight of species i, g/mol; and v_i are the atomic diffusion volumes (Fuller diffusion volume), cm³.

According to Darcy’s law, the velocity field is influenced by factors such as the pressure gradient, the viscosity of the fluid, and the architecture of the porous material:

$$u=-\frac{k}{\eta }\nabla p$$

(19)

In this equation, u is the Darcy’s velocity or specific discharge vector, m/s; $k$ is the permeability of the porous medium m²; $\eta$ is the fluid’s dynamic viscosity Pa·s;

The mixture viscosity of the gas phase is based on the kinetic theory by Brokaw [36]:

$${\eta }_v=\sum _i\left(\frac{x_i{\eta }_{i,v}}{\sum _j{x}_j{\psi }_{i,j}}\right)$$

(20)

where

$${\psi }_{i,j}=A_{i,j}\sqrt{\frac{{\eta }_{i,v}}{{\eta }_{j,v}}}$$

(21)

It can be obtained by Herning and Zipperer approximation as molar mass ratio [37] or directly via interaction parameter $A_{i,j}$. We use the last-mentioned approach.

$$A_{i,j}=\frac{{\beta }_{i,j}}{\sqrt{R{M}_{i,j}}}\left(1+\frac{R{M}_{i,j}-R{M}_{i,j}^{0.45}}{2\left(1+R{M}_{i,j}\right)+\frac{1+{RM}_{i,j}^{0.45}}{\sqrt{{\beta }_{i,j}}(1+R{M}_{i,j})}}\right)$$

(22)

with

$${\beta }_{i,j}={\left(\frac{4M_i{M}_j}{{\left(M_i+M_j\right)}^2}\right)}^{1/4},R{M}_{i,j}=\frac{M_i}{M_j}$$

(23)

where ${\eta }_{i,v}$ is the vapor viscosity of each individual species, and $M_i$ are the molecular weights.

2.4. Heat Transfer

A prevalent assumption in SOFC modeling is the premise of local thermal equilibrium (LTE) [38,39,40]. In the porous catalyst layer, this approach enables the use of effective heat conductivity and heat capacity, which can be calculated as follows [41,42]:

$$k_{eff}=\epsilon ·k_f+(1-\epsilon )k_f$$

(24)

$$c_{p,eff}=\epsilon ·c_{p,f}+(1-\epsilon )c_{p,f}$$

(25)

where ε is the porosity; $eff$ means effective; $s$ means solid and $f$ means fluid (gas) phase.

Yet, certain common conditions observed in the porous SOFC electrodes challenge this assumption: (1) extremely low flow with a minimal Reynolds number, (2) occurrence of heat generation in volume, and (3) a significant disparity in thermal conductivity between the gaseous and solid phases.

The gas phase mixture thermal conductivity, which may be used in heat transfer simulations, is calculated according to:

$${\lambda }_v=\sum _i\left(\frac{x_j{\lambda }_{i,v}}{\sum _j{x}_j{\phi }_{i,j}}\right)$$

(26)

$${\phi }_{i,j}=\frac{1}{4}{\left(1+{\left(\frac{{\eta }_{i,v}}{{\eta }_{j,v}}{\left(\frac{M_j}{M_i}\right)}^{\frac{3}{4}}\frac{T+\frac{3}{2}{T}_{b,i}}{T+\frac{3}{2}{T}_{b,j}}\right)}^{1/2}\right)}^2\left(\frac{T+S_{ij}\sqrt{\frac{9}{4}{T}_{b,i}{T}_{b,j}}}{T+\frac{3}{2}{T}_{b,i}}\right)$$

(27)

In the aforementioned equation, the summations account for all active species. Here, λ_i,_v and η_i,v represent the individual vapor thermal conductivity and dynamic viscosity correlations based on temperature for the pure gases, respectively. The formula also takes into account the normal boiling points T_b,i, and the molecular weights of the species. The coefficient S_ij is set at 0.773 when either i or j corresponds to the index of water (or steam). In other scenarios, S_ij is set to 1.

Given that this equation is dependent on the viscosities of individual species, thermal conductivity is determined solely in gas phase domains using the built-in dynamic viscosities, without any user-defined modifications.

The heat capacity is calculated as:

$$c_{p,eff}=\sum _i\frac{{{\omega }_{i}c}_{p,i}}{M_i}$$

(28)

$c_{p,i}$ is the species heat capacity in J/(mol·K).

According to the manufacturer’s specifications, the specific heat for the YSZ + NiO anode and the LSCF + GDC cathode is 450 J/kgK each.

2.5. Model for Additional Reactions

In order to model SOFC on methane several additional reactions are required. The first one is methane reforming. Following this is the Water Gas Shift (WGS) reaction.

The reforming process of methane is characterized by the following two equations

$$\mathrm{C}{\mathrm{H}}_4+{\mathrm{H}}_2\mathrm{O}=3{\mathrm{H}}_2+\mathrm{C}\mathrm{O}$$

$$\mathrm{C}\mathrm{O}+{\mathrm{H}}_2\mathrm{O}={\mathrm{H}}_2+\mathrm{C}{\mathrm{O}}_2$$

Steam reforming is widely recognized as a highly endothermic reaction, whereas the water-gas shift reaction exhibits mildly exothermic characteristics. In our analysis, the focus on the reaction rate is mainly on steam-methane reforming. This is because, at elevated temperatures, the water-gas shift reaction is swift and achieves equilibrium nearly immediately [43,44]. Therefore, for simplicity, we assume the water-gas shift reaction rate as 10,000 [mol/(m³·s)].

As highlighted earlier, three formulas are utilized to outline reaction kinetics. To sidestep intricacy and maintain computational efficiency, we employ the subsequent expression in our study. This expression incorporates the partial pressures of both methane and steam to depict the kinetics.

$$R_{smr}=k_{smr}{\left(P_{CH4}\right)}^a{\left(P_{H2O}\right)}^b$$

(29)

$$k_{smr}=Aexp(-\frac{E}{RT})$$

(30)

In Equation (29),

R_smr is the methane consumption rate per second and per unit exposed Ni area; k_smr is the reaction coefficient k_smr = 1.42 × 10³, mol/(s·m²·Pa) [45] P_CH4 and P_H2O are the partial pressures of methane and steam at the reaction site, respectively.

In Equation (30),

A is the frequency factor A = 1.42 × 10³ [mol/(s·m²·Pa)]; and E is the activation energy of the reaction. E = 82 × 10³ [J/mol] [43]

By quantifying the area density of Ni, we can represent the rate of methane consumption in the subsequent manner:

$$R_{chvol}=R_{chvol}\ast S_{Ni-pore}$$

(31)

$S_{Ni-pore}$ represents the cumulative area of Ni oriented towards the pores in the sample. S_Ni−pore is obtained as follows:

$$S_{Ni-pore}=\frac{W}{\rho A_{Ni-pore}}$$

(32)

W—SOFC volume, obtained from geometrical parameters.

So, according to the reactions and stoichiometric coefficients, we can obtain the consumption/production

$$R_i=v\ast Rchvol\ast P_{CH4}\ast M_i$$

(33)

The heat of the reaction is calculated via enthalpy change as

$$Q=R\ast ∆H$$

(34)

The enthalpy change for methane reforming is calculated according to the method depicted detailed in the reference [45]. Water Gas-shift reaction enthalpy change is ∆H = −41.6 kJ/mol [46].

2.6. Model Multiphysics

To simultaneously solve the aforementioned equations while accounting for non-isothermal flow, we utilize additional modules found in the multiphysics section of COMSOL Multiphysics.

In Multiphysics, the Reacting Flow module facilitates a two-directional coupling mechanism. This allows the velocity field and pressure layout, as determined by the fluid flow interface, to be integrated into the Fuel Cell module. Simultaneously, parameters like density and dynamic viscosity, as deduced by the Fuel Cell module, are incorporated into the fluid flow interface.

The Electrochemical Heating module within Multiphysics is designed to specify domain and boundary heat sources in a heat transfer interface. This is based on the aggregate of irreversible heat components (which encompasses Joule heating and activation losses) and reversible heat within an electrochemical interface. Additionally, this module aligns the temperature in the electrochemical interface with that of the heat transfer interface.

The Nonisothermal Flow module in Multiphysics is employed to address the conservation principles of energy, mass, and momentum in fluids and porous structures, along with energy conservation in solid materials. The described conditions are shown schematically in Figure 2

2.7. Model Inputs

A comprehensive model of a solid oxide fuel cell requires various inputs that describe the materials, operating conditions, physical geometry, and other factors.

Porosity: The porosity of the electrodes and electrolytes affects the diffusion of gas species within these materials. The porosity must be specified for each material in the cell. ${por}_a=0.36;{por}_c=0.39$.

Permeability: This is a measure of the ability of a porous material to allow fluids to pass through it. The permeability of the anode and cathode must be defined. ${perm}_a=2·{10}^{-10}\left[{\mathrm{m}}^2\right];{perm}_c=2·{10}^{-10}\left[{\mathrm{m}}^2\right]$.

Operating Temperature: The temperature at which the fuel cell operates. In this case, the operating temperature was about 900 degrees Celsius.

Inlet parameters:

Air: 0.23 O₂, 0.765 N₂, 0.005 H₂O; 2/15 L/min; Tin = 650 °C.

Methane fuel: 0.995 CH₄, 0.005 H₂O; 1/15 L/min; Tin = 900 °C.

Syngas fuel: 0.08 N₂, 0.067 H₂, 0.225 CO, 0.25 CO₂, 0.245 CH₄, 0.133 H₂O; 1/15 L/min; Tin = 900 °C.

2.8. Numerical Model Data

2.8.1. Mesh

When constructing the mesh, we set a condition that its step should be 12 times smaller than the thickness of each of the modeled elements. Thus, since the electrolyte is the thinnest element with a thickness of 15 μm, the maximum mesh step across the thickness of this element was 1.25 μm.

Additionally, the mesh step was set to decrease closer to the boundaries of the elements.

Along the length, the mesh step was 1 mm, and in width, it was 0.15 mm.

The final mesh consists of 112,752 domain elements, 26,928 boundary elements, and 2188 edge elements. Its appearance is shown in Figure 3.

2.8.2. Mesh Sensitivity Analysis

To ensure the accuracy and reliability of our numerical simulations for the SOFC model, a grid sensitivity analysis was conducted. This step is pivotal in confirming that our results are not influenced by the mesh size and that the chosen mesh captures the intricate physical phenomena within the SOFC effectively.

Methodology:

Three distinct mesh configurations were considered:

Coarser Configuration:

34,560 domain elements, 12,192 boundary elements, and 1464 edge elements.

Minimum Element Size: 1.875 μm

Maximum Element Size: 1.5 mm

Standard Configuration:

112,752 domain elements, 26,928 boundary elements, and 2188 edge elements.

Minimum Element Size: 1.25 μm

Maximum Element Size: 1 mm

2× Configuration:

866,880 domain elements, 105,244 boundary elements, and 4344 edge elements.

Minimum Element Size: 0.625 μm

Maximum Element Size: 0.5 mm

The primary metrics for comparison across these grids were temperature distribution, current density, average cell power, and reactant concentrations. The aim was to ascertain the deviation in these metrics across different mesh sizes and determine the optimal configuration.

Results:

Comparing the Standard and Coarser configurations, the difference in maximum temperature was about 3% for all fuels. For other parameters, the difference was less than 2%. Between the Standard and 2× configurations, the difference was less than 0.001%, which is close to the specified tolerance.

Computational Time:

Standard Configuration: 21 min

1.5× Configuration: 5 h 19 min

3× Configuration: 32 h

Conclusion:

Based on the grid sensitivity analysis, the 1.5× mesh was selected for further simulations presented in this paper. This choice strikes a balance between computational efficiency and the precision required for SOFC simulations, ensuring the robustness of our findings.

2.8.3. Solver Configuration and Convergence

The precision and efficiency of our numerical simulations are contingent upon the solver configuration. This section elucidates the solver settings and parameters employed for our stationary SOFC simulations in COMSOL Multiphysics.

Linear Solver

Type: Direct.

Method: MUMPS (Multifrontal Massively Parallel Sparse Direct Solver).

Preconditioner: ILU (Incomplete LU factorization).

Non-linear Solver:

Method: Newton.

Convergence Criterion: Residuals of all equations dropping below a tolerance of 1 × 10⁻⁶.

Maximum Iterations: 50.

Parameter Variation:

For our simulations, the cell voltage was varied from 0.05 to 0.95 volts in increments of 0.1 volts.

2.8.4. Hardware Specifications

Processor: AMD Ryzen 7 3700X

RAM: 64 GB DDR4

Storage: 2TB NVMe SSD

Computational Time:

The total computational time for the simulations to converge was about 5 h for hydrogen fuel and about 6 h for methane and syngas fuels.

3. Results and Discussion

Upon running the simulations with the defined parameters and boundary conditions, we obtained the distributions of temperature, current, and reactant concentrations across the fuel cell. From there, we can also calculate the overall power output. Let us delve into the details.

3.1. Temperature Distribution

From the temperature distribution illustrated in Figure 4, we can observe that as cold air travels through the channel, it warms up from an initial temperature of 923 K to a maximum of 1556 K near the channel’s midpoint. This heating is a result of the electrochemical processes occurring within the fuel cell, as previously described. Moreover, it can be seen that towards the end of the channel, the air temperature drops to approximately 1250 K. This temperature decrease is attributable to two factors: the diminishing intensity of electrochemical processes due to reduced air concentration, and the cooling effect caused by heat conduction from the proximity to the hydrogen inlet.

This illustration also distinctly depicts the forward shift of the temperature gradient along the channel relative to the electrode temperature, an effect resulting from heat dissipation due to the movement of the medium.

When comparing different types of fuels, it is evident that carbonaceous fuels lead to approximately 70 K less heating, which is attributed to endothermic reactions. This also impacts the performance of the fuel cell, as will be demonstrated later.

A temperature profile similar to ours has also been observed by other researchers. The temperature peak is typically seen closer to the center of the fuel cell [9].

3.2. Current Distribution

Regions with high current density indicate areas where fuel utilization is efficient (Figure 5). Essentially, the current distribution serves as a map showing the zones where electrochemical processes are actively occurring. As the fuel becomes depleted, there is a corresponding decrease in current density.

3.3. Reactants Concentration

In this section, we examine the component fractions while operating on methane and synthesis gas, as these represent the most intriguing scenarios.

In the case of methane (Figure 6), we observe a progressive decrease in its concentration due to the ongoing reforming reaction. This reaction produces hydrogen and carbon monoxide, causing an increase in their concentrations along the channel length.

As the concentration of carbon monoxide increases, the water-gas shift reaction becomes more active, leading to a rise in carbon dioxide concentration, as shown in Figure 5.

Towards the end of the channel, with a scarcity of methane, the rate of carbon monoxide production becomes slower than its consumption due to the water-gas shift reaction. This results in a decreased carbon monoxide concentration, as illustrated in Figure 6. We notice an intriguing situation regarding hydrogen distribution. As a result of the ongoing reforming and water-gas shift reactions, its concentration increases, peaking towards the middle of the channel. However, with a depletion in reactants on one hand, and an acceleration of the electrochemical reaction on the other, the hydrogen proportion starts to decline. This process is shown in Figure 7.

When using a synthesis gas mixture, we witness similar reactions. Nevertheless, due to the varying initial concentrations of components, their distribution along the channel length differs compared to methane fuel.

For instance, the high flow rate of the water-gas shift reaction and the initial presence of its components result in a concentration peak of hydrogen at the beginning of the channel. Subsequently, as a result of electrochemistry, the concentration of hydrogen falls (Figure 8).

The high rate of the water-gas shift reaction also impacts the distribution of carbon monoxide and carbon dioxide, as evidenced in Figure 9. Initially, carbon monoxide rapidly transforms into carbon dioxide. However, with a lack of water vapor, this reaction concludes swiftly. As the concentration of water vapor and carbon monoxide increases, the concentration of carbon dioxide begins to rise along the length of the channel, reaching its peak at the end. This is due to the absence of reactions that would consume carbon dioxide.

3.4. Power Output, Model Verification

The average current density of the cell is calculated using COMSOL tools. In this instance, the z-component of the computed current density vector in the electrolyte [A/cm²] is averaged over the volume of the electrolyte. Averaging is achieved through fourth-order volume integration. The power output of the cell (power density, to be precise) is calculated by multiplying the average current density of the cell by the cell voltage.

As shown in Figure 10, the predicted current density versus voltage from the numerical model demonstrates acceptable accuracy compared to the manufacturer’s data. For the majority of the operating range, the deviation does not exceed 6 percent. It also aligns with the results obtained by other researchers [47,48,49]. However, simulations indicate different current densities at lower voltages.

We assume that this deviation is a result of the fact that at high current density, the heating of the fuel cell should occur more intensely than it does in our model. This heating should be greater because, when modeling a single channel, the influence of adjacent channels in the assembly is not taken into account. In our case, the initial temperature of the fuel has a stronger influence on the temperature distribution than it would in a multi-channel model. The elevated temperature in the multi-channel model would inevitably lead to the acceleration of electrochemical reactions and increased current density.

We plan to test this hypothesis in our next study, where the entire multi-channel assembly will be modeled.

Finally, the overall power output of the SOFC can be calculated from the voltage and total current. This gives us a measure of the cell’s energy conversion efficiency, enabling comparison with experimental data or other models. The results are displayed in Figure 11 and Figure 12.

It is noteworthy that the fuel cell demonstrates the highest power output when operating on hydrogen. A slight decrease in power for methane is attributed to the endothermic reaction of steam reforming, which leads to a reduction in the temperature of the fuel cell and a slowdown in electrochemical reactions. At the same time, the power reduction is minimal, indicating that the rate of hydrogen generation from steam reforming and water-gas shift reaction is not less than the consumption of hydrogen due to electrochemical reactions. Thus, since there is no shortage of hydrogen, the overall power remains close.

Interestingly, syngas proves to be a more efficient fuel than methane, as it facilitates a faster production of hydrogen through the water-gas shift reaction.

4. Conclusions

This article presents a comprehensive model of a solid oxide fuel cell (SOFC) using the COMSOL Multiphysics 6.1 software. The model considers key physical phenomena such as electron and ion transport, gas diffusion, electrochemical reactions, and heat transfer, and it integrates these processes to accurately simulate the operation of an SOFC.

The model inputs include material properties, operating conditions, and geometric parameters which are grounded on experimental data and literature. We have incorporated detailed considerations for the calculation domain, boundary conditions, and multiphysics interactions, providing an extensive simulation platform.

Carbonaceous Fuels vs. Hydrogen: The fuel cell, when operating on carbonaceous fuels, exhibited a temperature reduction due to endothermic reactions, contrasting the behavior observed in hydrogen-fueled cells. The fuel cell delivers higher power when operating on hydrogen, producing 1340 W/m² at a current density of 0.25 A/cm². For methane, the power output is 1200 W/m² at a current density of 0.25 A/cm². For synthesis gas, the maximum power of 1340 W/m² is achieved at a current density of 0.3 A/cm².

Temperature Distribution: Our simulations revealed distinct temperature profiles within the cell. The central part of the channel in hydrogen-fueled SOFCs was identified as the high-temperature zone, with a peak temperature that was about 70 K higher than that of cells fueled by carbon.

Reactants Concentration: The concentration of components in the developed model displays the anticipated physical behavior of the fuel cell operation. When operating on carbonaceous fuels, a change in substance concentrations is observed due to reforming reactions and electrochemical reactions.

The developed model was compared with a 1-kW industrial prototype operating on hydrogen and showed good agreement in the volt-ampere characteristic with a deviation not exceeding 5% for the majority of the operating range.

The work presented here is an important step towards more realistic modeling of SOFCs. It provides a foundation for future studies to further investigate and optimize SOFC performance under various conditions. Although the model currently employs certain assumptions, these can be refined in future studies to provide even more precise simulations. A clear enhancement to the model would be the incorporation of multi-channel geometry and accounting for thermal radiation, which we plan to implement in subsequent studies.