Numerical Simulation of Processes in an Electrochemical Cell Using COMSOL Multiphysics

: Fuel cells are a promising source of clean energy. To ﬁnd optimal parameters for their operation, modeling is necessary, which is quite difﬁcult to implement taking into account all the signiﬁcant effects occurring in them. We aim to develop a previously unrealized model in COMSOL Multiphysics that, on one hand, will consider the inﬂuence of electrochemical heating and non-isothermal ﬂuid ﬂow on the temperature ﬁeld 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 ﬁndings 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 2 with a current density of 0.25 A/cm 2 . When fueled by methane, it produces 1200 W/m 2 at the same current density. Using synthesis gas, it reaches its peak power of 1340 W/m 2 at a current density of 0.3 A/cm 2 .


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].
Energies 2023, 16, 7265 2 of 20 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.

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 2 , with a designated active region measuring 10 × 10 cm 2 .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.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: 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 COM-SOL 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  = 333330 [S/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  of 7937 [S/m].

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.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: 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 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 σ c of 7937 [S/m].

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.
Anode : Energies 2023, 16, 7265 4 of 20 The electromotive force of the aforementioned electrochemical system is depicted by the Nernst equation [19]: E eq -the reversible voltage or Equilibrium potential in terms of COMSOL; x v pt and x v rt -the concentration of the reactants and products at the reaction sites, and the superscript v stands for the stoichiometric coefficient where E eq,re f 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: COMSOL Multiphysics employs the same approach, but in this case, E eq,re f is determined based on changes in enthalpy and entropy, using data from NIST [21].

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: where i-the volumetric exchange current density, A/m 3 ; i 0 -is the equilibrium exchange current density per unit volume associated with the TPB reaction in the anode (A/m 3 ) when E = E eq and it is estimated based on the fitting to the experimental data [28]; η act is the activation overpotential, V; n is the number of electrons involved in the electrode reaction; α a and α 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: φ s -denotes the electric potential at the surface; φ l -is the electrolyte potential; E eq -denotes the equilibrium electrode potential; E th = φ s − φ l To utilize the Butler-Volmer equation, one needs to determine the equilibrium exchange current density, charge transfer coefficients, and TPB length.

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.
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 2 ), representing the exchange current density when Eeq matches E eq = E eq,ref .
So equilibrium exchange current density: However, we still need to specify i o,re f for anode and cathode.Based on the empirical formula [30]: According to (10) i o,re f _c = 0.11A/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.

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 3 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 2 /µm 3 , 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 2 /µm 2 .

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

•
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 Thermodynamic formulas for gaseous species of H i (T) at a standard partial pressure of 1 atm are derived from NIST datasets [21].

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: 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.

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.

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.

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: where ρ-denotes the mixture density, kg/m 3 , u the mass averaged velocity of the mixture, m/s; ω i -the mass fraction; j i -the mass flux relative to the mass averaged velocity, kg/(m 2 •s)]; R i -the rate expression describing its production or consumption, kg/(m 3 •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]: where ∼ D ik are the multicomponent Fick diffusivities, m 2 /s; D T i 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: In the case of ideal gas mixtures, the force driving diffusion is: [36]: 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 2 .For ionic species, the external force is generated by the presence of an electric field.

Binary Diffusion Coefficients
The determination of the inherent binary diffusion coefficients is carried out utilizing the Fuller-Schettler-Giddings (FGS) approach [36]: 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 3 .
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: In this equation, u is the Darcy's velocity or specific discharge vector, m/s; k is the permeability of the porous medium m 2 ; η is the fluid's dynamic viscosity Pa•s; The mixture viscosity of the gas phase is based on the kinetic theory by Brokaw [36]: where 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. with where η i,v is the vapor viscosity of each individual species, and M i are the molecular weights.

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]: where ε is the porosity; e f f 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: 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,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.

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 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 3 •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.
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 3 , mol/(s•m 2 •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 3 [mol/(s•m 2 •Pa)]; and E is the activation energy of the reaction.E = 82 × 10 3 [J/mol] [43] By quantifying the area density of Ni, we can represent the rate of methane consumption in the subsequent manner: S Ni−pore represents the cumulative area of Ni oriented towards the pores in the sample.S Ni−pore is obtained as follows: W-SOFC volume, obtained from geometrical parameters.So, according to the reactions and stoichiometric coefficients, we can obtain the consumption/production The heat of the reaction is calculated via enthalpy change as 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].

Model Multiphysics
To simultaneously solve the aforementioned equations while accounting for nonisothermal flow, we utilize additional modules found in the multiphysics section of COM-SOL 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 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

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. = 0.36;  = 0.39.

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.

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.

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.
Along the length, the mesh step was 1 mm, and in width, it was 0.15 mm.The final mesh consists of 112,752 elements, 26,928 boundary elements, and 2188 edge elements.Its appearance is shown in Figure 3.

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.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.

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: 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.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.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].

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.

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.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].

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.

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.
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.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.
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.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).
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.

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 2 ] 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 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.
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.

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 2 ] 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

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 2 ] 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.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 multichannel 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 Figures 11 and 12   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 multichannel 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 Figures 11 and 12.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 multichannel 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 Figures 11 and 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.

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 2 at a current density of 0.25 A/cm 2 .For methane, the power output is 1200 W/m 2 at a current density of 0.25 A/cm 2 .For synthesis gas, the maximum power of 1340 W/m 2 is achieved at a current density of 0.3 A/cm 2 .
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.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.

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 2 at a current density of 0.25 A/cm 2 .For methane, the power output is 1200 W/m 2 at a current density of 0.25 A/cm 2 .For synthesis gas, the maximum power of 1340 W/m 2 is achieved at a current density of 0.3 A/cm 2 .
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.

Figure 1 .
Figure 1.Schematic model of the fuel cell indicating processes and boundary conditions.

Figure 1 .
Figure 1.Schematic model of the fuel cell indicating processes and boundary conditions.

Figure 2 .
Figure 2. Schematic model of the fuel cell indicating multiphysics conditions.

Figure 2 .
Figure 2. Schematic model of the fuel cell multiphysics conditions.

Figure 3 .
Figure 3. Computational mesh of the model.

Figure 3 .
Figure 3. Computational mesh of the model.
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-49].However, simulations indicate different current densities at lower voltages.themajority 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.

Energies 2023 ,
16, x FOR PEER REVIEW 16 of 20 the majority of the operating range, the deviation does not exceed 6 percent.It also aligns with the results obtained by other researchers [47-49].However, simulations indicate different current densities at lower voltages.