The mathematical, multiphysic, multidimensional, and electrochemical modelation of a high temperature solid oxide fuel cell system (planar electrolyte-supported configuration) is discussed in the present paper. The mass transport within the cell is studied using the Stefan-Maxwel model, and the momentum balance is solved by means of Navier-Stokes and Brinkman equations, respectively. On the other hand, the energy balance includes the generation term coupled with the convection and conduction equations. It was demonstrated that the diffusion resistances play an important role in the cell performance, and the oxidant concentration is enough high to work at fuel utilization coefficient of 0.8. The current density suffers a reduction (10 A/m2 to 1.5103 A/m2) due to the variation of reactants concentration at the cell outlet and the diffusive flux resistances. The developed models can be used to further analyses and to study a solid oxide fuel cell working with other fuels but hydrogen.

1. Introduction

The use of fossil fuels is a one of the major concerns due to its effect on the environment pollution and global warming. The fuel cells (FCs) are electrochemical devices which can produce electric power and heat, by chemical combination among a fuel (generally H2) and an oxidizer (O2) [1]. The solid oxide fuel cells are promising candidates for future energy systems for their higher energy efficiency than heat engines and other fuel cell models. In addition the SOFCs, allows working at high temperatures (650–1000°C), supporting the internal conversion of light alcohols, hydrocarbons, and other substances that are considered poisons for the other types of fuel cells (CO) [2, 3].

A simple SOFC system is composed of two electrodes (anode and cathode), separated by an electrolyte through the one to which the ions O2− are transferred. The current is generated fundamentally in the anode by means of the direct conversion of the chemical energy of the fuel [4]. The SOFC are generally coupled in stacks (packages of cells) with different configurations, according to the operation conditions [5]. The planar design has received special attention recently, because it is simpler to fabricate and it exhibits higher current density than cylindrical designs [3].

As can be corroborated the low negative effect on the environment and higher energy efficiency than traditional generating technologies made the fuel cells one of the candidates for the near energy future [6].

Nowadays this devices are built of ceramic and cermet materials (metallic-ceramic) [7], generally Ni/Yttria-stabilized zirconia YSZ (anode), perovskite and lanthanum strontium manganese oxide with YSZ (cathode), while to conform the electrolyte it is used Yttria-stabilized zirconia (YSZ) which provides the appropriate ionic conduction (O2) [8]. The combination between charge, mass, and energy transport with the diffusive fluid flow within channels and electrodes is the main phenomena affecting the process efficiency. Moreover the overpotentials losses by activation, ohmic, and concentration are closely related with the total current density and polarization curves [911].

The aim of the present communication is to develop a comprehensive model which can be used to describe the processes mentioned above and to be used as design tool. The mathematical solution of the partial differential equations was performed using the finite element method to a microsolid oxide fuel cell with a planar geometry.

2. Mathematical Model

2.1. Problem Definition

Since a lot of complicated models are needed to simulate the whole cell performance in the present paper, a single unit was taken as reference case. The domain includes the electrodes, electrolyte, and the channels for the fuel and oxidizer, in an SOFC system working with H2 like fuel and O2 (air), as oxidizer. The following processes are described starting from the multiphysics combination of the phenomenological models as follows.(a)Species transport in channels and solid elements.(b)Momentum balance in channels and electrodes.(c)Dissociation and ionization of the oxidizer in the interface cathode/electrolyte.(d)Fuel reaction at interface anode/electrolyte.(e)Heat balance for overall domain.

2.2. Model Description

The conservation laws are applied to the geometry represented in Figure 1 to obtain a general model in two dimensions for the system. In the solid parts of the system (anode/electrolyte/cathode) the electrokinetics, charge, momentum, and mass transport phenomenon are mainly studied, while in the channels the transfer of heat, mass, and momentum is analyzed as well. The kinetic aspects of the electrochemical reactions and its contribution as source or sinks for transports equations are discussed in the active area of the electrodes. The electrochemical reactions that are considered are the following ones: H2+O2=H2O+2e1(Anodereaction),(1)2O2+2e=O2(CathodeReaction).(2)

The mathematical formulation and the solution procedure were simplified assuming isotropic and homogeneous electrodes, with uniform morphological properties (porosity and permeability), parallel flow between the fuel and the oxidizer, negligible ohmic drop in electronically conductive solid womb of porous electrodes and ideal gas mixtures.

Weight fractions of the fuel and oxidizer in the regions where the balances of mass are applied, potential in any point of the solid surfaces, the temperature in the whole cell and the falls of pressure in the channels and electrodes, respectively, are the main incognito in the modelation. The general operation conditions are summarized in Table 1.

2.3. Conservation Laws
2.3.1. Mass and Species Balance

The mass balance is based on the Stefan-Maxwell equation (3) [12], for the divergence of the mass flow for diffusion (porous elements) and convection (channels) taking into account the interactions among the dissolved species by means of the symmetrical diffusivity defined by (4) reported elsewhere [13]: 𝜌𝑤𝑖𝑁𝑗=1𝐷𝑖𝑗𝑀𝑀𝑖𝑗𝑤𝑗+𝑤𝑗𝑀𝑀+𝜌𝑤𝑖𝑢𝐷=𝑅𝑖,(3)𝑖𝑗𝑇=𝑘1,75𝑝𝑣𝑖1/3+𝑣𝑗1/321𝑀𝑖+1𝑀𝑗1/2.(4)

2.3.2. Momentum Balance [13, 14]

The momentum balance in the channels and the porous elements is described by the Navier-Stokes (5) [15, 16] and Brinkman (6) equations, respectively, in conjunction with the continuity law (7), conforming in this way a nonlinear system, of partial differential equations (PDE): 𝜌𝜕𝑢𝜕𝑡𝜇𝑢+(𝑢)𝑇𝜇𝑘𝜌𝑢+𝑝𝐹=0,(5)𝑑𝑢𝜇𝑑𝑡𝜇Δ𝑢+𝑘𝑢+𝑝=0,(6)𝑢=0.(7)

2.3.3. Potential and Charge Transfer [5, 17]

The potential variation in the solid elements of the cell is subject to fluctuations within the interfaces electrode/electrolyte, due to finite discontinuities that appear because of variations of species concentrations in those places [18]. The charge and potential balances applied to a planar solid oxide fuel cell are described by the Ohm law (8) adding a generation term (𝑅𝑝𝑖𝑗) in the active surfaces: 𝜎𝜙+𝑅𝑝=0,(8)𝑅𝑝=𝑖𝑗,𝑗=𝑎,𝑐.(9)

The exchange current densities at anode and cathode (10) and (11) are related with the irreversibility losses in the device mainly by concentration, activation, and ohmic. The total cell voltage is calculated considering the influence of the irreversibilities mentioned above: 𝑖𝑎=𝑖ref0𝐶H2𝐶H2,ref𝛾H2𝐶H2O𝐶H2Oref𝛾H2O𝛼exp𝐴𝐹𝜂,𝑅𝑇(10)𝑖𝑐=𝑖ref0𝐶O2𝐶O2,ref𝛾O2𝛼exp𝐶𝐹𝜂,𝑅𝑇(11)𝜂=𝜂𝑐+𝜂𝑎+𝜂Ohm.(12)

The activation overpotential is related to the electrode kinetics at the reaction site and the relationship between overpotential-current density can be expressed by the Butler-Volmer equation [9, 19, 20], which for a typical SOFC is expressed as 𝜂act=𝑅𝑇𝐹sinh1𝑖2𝑖𝑜,𝑗,𝑗=𝑎,𝑐.(13)

This expression is valid when 2 electrons are transferred in the electrochemical reaction, the symmetric factor of the SOFC (alpha) is 0.5 and could vary if the oxidations of CH4 and CO are considered [20].

On the other hand, the concentration overpotential is evaluated considering the limit current density, defined by Wang [21]. This parameter is closely related to the transport properties of the fuel and oxidant and the morphological characteristics of the cell electrodes: 𝜂conc=𝑅𝑇𝑖𝑛𝐹ln1𝑖𝑙,𝐷𝑖𝑙=0.68𝑛𝐹𝐶e2/3𝜇1/6𝐴𝑏1/2𝑈1/2.(14)

The effect of the Ohmic overpotential on the cell voltage is calculated using the equation presented by Ni et al. [9]; it is only affected by the electrolyte properties, the temperature, and the current density within the cell.

2.3.4. Heat Balance [21]

The general heat balance equation for convection and conduction is used to calculate the heat effects and temperatures profiles along the channels and electrodes, adding a source term defined as 𝑄, in the interfaces where the electrochemical reactions take place (15)–(17): 𝜌𝐶𝑝𝜕𝑇𝜕𝑡+𝑘ei𝑇+𝜌𝐶𝑝𝑇𝑢=𝑄,(15)𝑄=2𝐾e𝑇+𝑖𝜂+𝑖𝑑𝐸,𝐾𝑑𝑇(16)e=1𝜀𝑝/3𝐾𝑖+𝜀𝑝/(2𝐾𝑠+𝐾)2𝐾.(17)

2.3.5. Numerical Procedure

To solve the complex system of partial differential equations, the finite element method was implemented. The solution region was divided into several subregions, and the generated mesh in the workspace was refined to increase the effectiveness in all the calculations at the electrode/electrolyte boundary. The iterative procedure is chosen to solve the problem, establishing a minimum error between iterations of 10−6 and a total of 500 iterations.

3. Results and Discussion

The simulation of the system is carried out for a fuel utilization coefficient of 80% and oxygen around 30%. At these conditions and a fuel cell temperature of 923 K the species concentration, current density curve, and the convective flux influence on the mass balance within the device were studied. Figures 2 and 3 show how the oxidizer concentration (O2) in the air channel have a more marked decrease in comparison to the hydrogen (channel of the fuel); this phenomenon is closely related to the incidence of the flow field, the diffussional resistances, and the electrokinetics on the catalyst boundaries, those which on the whole affect the accessibility of H+ toward the reaction surface. The lowest concentration of both substances is located at the boundary between the electrodes and the electrolyte where the reaction is verified at higher extension than other zones within the cell.

All previously outlined depends on the species migration and the diffusive flux of hydrogen and oxygen in the channels which can be corroborated analyzing Figures 4 and 5, where it is appreciated as the decrease in the diffusive flux and causes an important fall in the accessibility of hydrogen to the active surface of the electrode; this problem is translated in a negative influence on the kinetics of the electrochemical process and consequently in a decrease of the current density at anode (Figure 5).

In Figure 5 the water concentration is increased proportionally to the loss in the hydrogen weight fraction. The arrows in the same figure allows identifying the influence of the concentration gradients in the process mass balances and the migration of the species toward the reaction zones placed at the electrodes surfaces. Even at these conditions the oxygen fraction is maintained enough high (~0.93%) to complete the hydrogen conversion.

The current density along the cell is represented in Figure 6. The reduction of the current density (10 A/m2 to 1.5 * 10−3 A/m2) is due to the lower reactants concentration at the cell outlet and the diffusive flux reduction while the process takes place. Also the concentration overpotential and the activation losses have a great influence on this behavior; similar results were reported in a previous paper [19]. The velocity profile in the channels is flat (Figure 7), demonstrating as the system of models describes the uniformity in the distribution fields and the absence of stagnant films within the cell. Similar results in this sense were reported by Beale, 2004 [17].

4. Conclusions

The developed model allows describing the phenomenon that take place inside an SOFC. The coupling of a multiphysics system of complex mathematical equations has been created taking into account the transport phenomena that describe the process and the kinetic aspects of the electrochemical reactions inside the device. Although the study is even preliminary, it is considered an advance step in the development of useful tools in the reliable and good design of fuel cells. If other fuel than hydrogen is used, the conversion reaction of methane, carbon monoxide, or other compounds should be included in the generation terms defined in (3) and (15).


𝜌:Density (Kg/m3)
𝑤𝑖𝑗:Weight fraction of component 𝑖
𝑣,𝑢:Velocity (m/s)
𝑝:Pressure (Pa)
𝑇:Temperature (K)
𝜇:Viscosity (Pa s)
𝑔:Constant of graveness (9.8 m/s2)
𝐶𝑝:Heat capacity P cte (J/KgK)
𝐾e:Gaseous mixture thermal conductivity (W/mK)
𝐷𝑖𝑗:Gaseous mixture binary diffusivity (m2/s)
𝑖0+,𝑖0:Reference exchange current density for anode and cathode (A/m2)
𝑖𝐿:Current limit density (A/m2)
𝜂:Over potential (V)
𝑖𝐴, 𝑗𝐶:Current density at anode and cathode (A/m2)
𝜎:Electric conductivity (Ω1m1)
𝑈:Fuel flow (m3/s)
𝐸:Electric potential from Nerst equation(V)
𝐹:Faraday constant (96487 C/mol)
𝛼:Electronic transfer coefficient
𝑅:Constant universal of the gases (8.32 J/Kmol)
𝛿:Film thickness (m)
𝐶𝑖:Concentration of species 𝑖 (mol/m3)
𝑀𝑖:Molecular weight species 𝑖 (Kg/mol)
𝐷𝑖𝑗e:Effective diffusivity (m2/s)
Φ:Potential (V)
𝐾𝑖𝑗𝐾𝑠:Thermal conductivity for gas mixtures and solids (W/mK).