Numerical analysis of mass and heat transport in proton-conducting SOFCs with direct internal reforming

Applied Energy 149 (2015) 161–175

Contents lists available at ScienceDirect

Applied Energy journal homepage: www.elsevier.com/locate/apenergy

Numerical analysis of mass and heat transport in proton-conducting SOFCs with direct internal reforming Vikram Menon a, Aayan Banerjee a, Julian Dailly b, Olaf Deutschmann a,⇑ a b

Institute for Chemical Technology and Polymer Chemistry, Karlsruhe Institute of Technology (KIT), 76131 Karlsruhe, Germany European Institute for Energy Research (EIFER), Emmy-Noether-Str. 11, 76131 Karlsruhe, Germany

h i g h l i g h t s  Computationally efficient numerical model for H-SOFCs.  Model is validated with two sets of experimental data for H2O–H2 inlet mixtures.  Novel electrochemical model to represent charge transfer.  Temperature profiles, species transport, and electrochemistry are studied.  The effect of anode zoning is investigated.

a r t i c l e

i n f o

Article history: Received 26 November 2014 Received in revised form 26 February 2015 Accepted 7 March 2015

Keywords: Solid Oxide Fuel Cell (SOFC) Proton conducting Direct internal reforming Numerical modeling Reaction kinetics

a b s t r a c t A computational model to investigate proton-conducting Solid-Oxide Fuel Cells (SOFCs) with direct internal reforming is developed. The numerical framework employs a 42-step elementary heterogeneous mechanism for Ni catalysts, using mean-field approximation. Mass transport through the porous media is described by the dusty gas model (DGM). Electrochemical parameters are deduced by reproducing two sets of experimental data, via the non-linear Butler–Volmer equation. A simple 1-D energy balance model is used to predict temperature profiles. The performance of the cell is analyzed by assuming the co-flow planar cell to be adiabatic. Simulations are carried out to understand the influence of various operating conditions on temperature distribution, species transport, and electrochemistry in the cell. The effect of dividing the anode into four zones, with different specific catalytic areas, on macroscopic performance parameters is investigated. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Fuel cells deliver an efficient approach for energy conversion via usage of chemical energy. Unlike low temperature fuel cells, high temperature SOFCs have enhanced kinetics and reduced resistivity of the solid electrolyte. Additionally, high temperatures also facilitate the usage of cheaper catalysts that allow direct internal reforming within the cell [1]. Direct internal reforming reduces system complexity and cost (lack of need for an external reforming unit), and promotes intrinsic thermal-coupling between endothermic and exothermic reactions within the cell. This coupling results in enhanced system efficiencies due to a shift in the reforming reaction towards hydrogen production for oxide-ion-conducting SOFCs, as steam is also produced through electrochemical reactions in the anode side. Furthermore, the produced carbon monoxide ⇑ Corresponding author. Tel.: +49 721 608 43064; fax: +49 721 608 44805. E-mail address: [email protected] (O. Deutschmann). http://dx.doi.org/10.1016/j.apenergy.2015.03.037 0306-2619/Ó 2015 Elsevier Ltd. All rights reserved.

enables added hydrogen production via the water–gas shift reaction (COðgÞ þ H2 OðgÞ CO2 ðgÞ þ H2 ðgÞ). However, some of the disadvantages of internal reforming over Ni/YSZ anodes include coking that leads to catalyst deactivation, and large thermal gradients that cause thermal mismatch/stress due to variable expansion rates of different materials [2]. In the former case, carbon deposition/removal mainly occurs through the disproportionation of carbon monoxide (Boudouard reaction – 2COðgÞ CðsÞ þ CO2 ðgÞ), cracking of methane (CH4 ðgÞ CðsÞ þ 2H2 ðgÞ) and oxidation of carbon (CðsÞ þ O2 ðgÞ CO2 ðgÞ) [3]. The interplay of physico-chemical phenomena in SOFCs is complex. Modeling the system requires immense understanding of the coupled interactions at work. Traditional SOFCs employ an oxideion conducting electrolyte, such as yttria-stabilized zirconia (YSZ), which facilitates the electrochemical production of H2O at the anode three-phase boundary (TPB). There is a great deal of literature available on the modeling of these systems, ranging from simple 0-D electrochemical cells to quasi-2-D unit cells to complex

162

V. Menon et al. / Applied Energy 149 (2015) 161–175

Nomenclature Aact Ac As Asp Bg Cp dp dpore D Dekl Dh Dek;Kn Ecell Erev F Gz H Hc DH H_ I Jk Kg ne n_ p Pe Pr Q_ R Rct Re s_ DS t T V VOC Vtn

active anode surface area for electrochemistry (m2) area of cross section of flow channels (m2) area of the solid MEA (m2) specific catalytic area (m1) permeability (m2) specific heat (J kg1 K1) particle diameter (m) pore diameter (m) diffusivity (m2 s1) effective binary diffusion (m2 s1) hydraulic diameter (m) effective Knudsen diffusion (m2 s1) cell voltage (V) reversible cell potential (V) Faraday constant (C mol1) Graetz number heat transfer coefficient (J m2 K1 s1); specific enthalpy (J kg1) channel height (m) enthalpy of formation (J mol1) mixture enthalpy (J mol1) current density (A cm2) species flux (mol m2 s1) number of gas-phase species number of electrons involved molar flow rate (mol s1) pressure (Pa) MEA active perimeter (m) Prandtl number heat source term (J m3 s1) gas constant (J mol1 K1) charge transfer resistance (X cm2) Reynolds number molar production rate (mol m2 s1, mol m3 s1) entropy change (J K1) time (s) temperature (K) volume (m3) open-circuit voltage (V) thermo-neutral voltage (V)

3-D stacks [4–7]. Oxide-ion-conducting SOFCs can also be operated in ‘‘reverse’’ mode at high temperatures as solid oxide electrolysis cells (SOECs) for the production of H2 or H2/CO by means of H2O electrolysis or H2O/CO2 co-electrolysis, respectively [8–10]. Operation of SOFCs with proton-conducting electrolytes, such as BaCeO3-based ceramics, have emerged as systems that provide higher theoretical energy efficiencies than its oxide-ion-conducting counterpart due to higher average EMF under the same conditions [11,12]. Here, the electrochemical production of H2O occurs at the cathode three-phase boundary (TPB), which leads to exclusive utilization of H2 in the fuel channel. Thus, unless the system is run with fuels other than pure hydrogen, the problem of gas separation at the exit of the fuel channel can be avoided. Also, direct recycling of the anode tail-gas to the inlet is possible. The sphere of proton-conducting solid-oxide fuel cell technology offers promising new ways of dealing with the energy crisis. Extant literature gives insight into the progress made in this field. An electrochemical model was developed to study the influence of micro-structural parameters on cell performance [13]. The computed overpotentials were compared to its oxide-ion-conducting counterpart to understand the effect of using different electrolytes. The same model was extended to investigate methane fed SOFCs

Wk [X] x, y, z X Y

molecular weight of kth species (kg mol1) concentration (mol m3) co-ordinate direction (m) mole fraction mass fraction

Greek letters d Kronecker delta symbol g overpotential k thermal conductivity (J m1 s1) h surface coverage fraction l viscosity (kg m1 s1) q density (kg m3) r conductivity (S m1) s tortuosity t velocity (m s1) u porosity Subscripts ac air channel el electrolyte/electrode f fluid fc fuel channel k species index s solid (porous media – MEA) Abbreviation ASR area-specific resistance B–V Butler–Volmer DGM dusty-gas model H-SOFC proton-conducting solid oxide fuel cell LHV lower-heating value MEA membrane electrode assembly OCV open-circuit voltage RWGS reverse water–gas shift SOC solid-oxide cell TPB three-phase boundary WGS water–gas shift

with a proton-conducting electrolyte, and dispelled previous notions about oxide-ion-conducting SOFCs having lower performance [14]. This was mainly due to higher ohmic resistance of the proton-conducting electrolyte. A 2-D CFD model was used to probe the effect of inlet mass flow rates on electrochemical performance parameters of SOFCs using a proton-conducting electrolyte [15]. All of the aforementioned studies revolved around protonconducting BaCeO3-based ceramics. A detailed electrochemical model was used for the performance analysis of a SOFC composed of a proton-conducting SrCeO3-based ceramic electrolyte and Pt electrodes [16]. The performance of the system was found to improve with a reduction in the thickness of the cathode and electrolyte, and by increasing temperature and pressure. A similar study was carried out in the same system with direct internal reforming, for methane-based inlet fuel compositions [17]. The effect of operating conditions on V–I behavior and species distribution in the channels was presented. Operation of proton-conducting SOFCs with ammonia has also been reported [18]. However, none of the aforementioned studies concerning protonconducting SOFCs take into account a detailed elementary mechanism to describe the thermo-catalytic chemical interactions in the anode, and neither do they describe the interactions in the cell

163

V. Menon et al. / Applied Energy 149 (2015) 161–175

during non-isothermal operation with direct internal reforming. Also, global B–V equations were used to represent charge transfer. Moreover, most reported models restricted validation to a single set of experiments, at one specified temperature. In this paper, a numerical model which accounts for coupled interactions between fluid flow, mass transport, heterogeneous chemistry, porous media transport, energy transport and electrochemistry is reported. The detailed electrochemical model, comprising of non-linear B–V equations, is validated with two sets of experimental data [19,20]. A modified B–V equation is derived to represent charge transfer at the anode/electrolyte interface, while a global B–V equation is used at the electrolyte/cathode interface. In addition, the estimated electrochemical parameters are employed in a quasi-2-D planar cell model for parametric analysis. The dependence of temperature distribution on cell voltage, cell length, specific catalytic area of the anode, and anode zoning is illustrated. This analysis is also extended to investigate species transport within the whole cell. The code is available as a part of the DETCHEM™ software package [21].

2.1. Channel flow Plug flow is assumed in the fuel and air channels as shown in Fig. 1 [22]. The equation for species continuity for the one-dimensional laminar flow is given by

@ðqf Y k Þ @ðqf tY k Þ Pe ¼ þ Jk W k ; @t @z Ac

Kg X @ðqf tÞ @ðqf ttÞ Pe ¼ þt Jk W k @t @z A k¼1 c

The basic principle of a proton-conducting SOFC is the same as that of an oxide-ion-conducting SOFC, except for the type of ion conducted by the electrolyte. The representation of the planar SOFC under consideration, in this study, is shown in Fig. 1. The hydrogen oxidation produces protons and electrons (Eq. (1)), while oxygen reduction requires these protons and electrons to produce steam (Eq. (2)). They take place through the following reactions: At the anode–electrolyte interface,

with W given as

At the cathode–electrolyte interface,

1 2Hþ ðelÞ þ O2 ðgÞ þ 2e H2 OðgÞ 2

ð2Þ

The net reaction reduces to,

1 H2 ðgÞ þ O2 ðgÞ H2 OðgÞ 2

ð3Þ

In the following sub-sections, the various models capturing the involved physics in H-SOFCs are described.

ð5Þ

where Pe is the perimeter associated with the electrochemically active membrane electrode assembly (MEA), qf is the fluid density, Yk is the species mass fraction of species k, t is the velocity, z is the axial position, Wk is the species molecular weight, Kg is the number of gaseous species, and Ac is the cross-sectional area of the channel. Assuming constant pressure in the channels, the density is calculated from the ideal gas equation

pW ¼ qf RT

ð1Þ

ð4Þ

The velocity in the channel can be calculated from the momentum equation,

2. Modeling approach

H2 ðgÞ 2Hþ ðelÞ þ 2e

k ¼ 1; . . . ; K g

W¼

Kg X Xk W k

ð6Þ

ð7Þ

k¼1

In Eqs. (4) and (5), Jk is the flux at the electrode channel interface, which is calculated using the dusty-gas model (DGM). Species molar fluxes depend on the heterogeneous chemistry within the porous-electrode structure and local current density i(z). The energy balance equation for the flow channels is given by [7]

@ðqf C pf T f Þ @ðtqf C pf T f Þ 4 ¼ þ hðT s  T f Þ Dh @t @z

ð8Þ

Here, Cpf is the specific heat capacity of the fluid. Tf and Ts are the temperatures of the fluid stream and solid phase (MEA), respectively. The geometry dependent Nusselt number Nu, which is used to determine the heat transfer coefficient h(z), is represented as an empirical formula [23]

Nu ¼

hDh k

Fig. 1. Schematic representation of a planar co-flow proton-conducting Solid Oxide Fuel Cell (SOFC).

ð9Þ

164

V. Menon et al. / Applied Energy 149 (2015) 161–175

 0:5386   1000 6:7275 Nu ¼ 3:095 þ 8:933 exp  Gz Gz

ð10Þ

where Dh is the hydraulic diameter, k is the thermal conductivity of the fluid, Gz is the Graetz number which is given in terms of the Reynolds number, Re, and Prandtl number, Pr, as

Dh RePr z

Gz ¼

ð11Þ

2.2. Porous media transport

used to estimate the binary diffusivities Dkl. Solution of Eq. (12) requires the reaction source terms s_ k and boundary conditions at the electrode–gas chamber and electrode–electrolyte interfaces. At the electrode–gas chamber interface the inlet mass fractions serve as the boundary condition, while at the electrode–electrolyte interface the chemical species fluxes are zero. The electrochemical reaction source terms, for electrochemically active species, are calculated from the current density and are accounted along with the chemical source terms. Thus, the fluxes of all species but hydrogen are zero at the anode–electrolyte interface, and is given by

i 2F

The reaction–diffusion equation for species transport in the electrodes is solved one-dimensionally along the electrode thickness, transverse to the direction of axial flow through the channel. The transient form of the equation is given by

At the cathode–electrolyte interface, the electrochemical fluxes amount to

@ð/qf Y k Þ @ðJ W k Þ ¼ k þ s_ k W k Asp ; @y @t

J H2 O ¼

k ¼ 1; . . . ; K g

ð12Þ

J H2 ¼ 

ð20Þ

i i ; J ¼ 2F O2 4F

ð21Þ

The total density of the fluid within the porous structure can be computed from Kg Kg X @ð/qf Þ @ðJ k W k Þ X s_ k W k Asp ¼ þ @t @y k¼1 k¼1

2.3. Solid-phase energy balance

ð13Þ

Here, s_ k is the heterogeneous molar production rate of the chemical species k, y is the independent spatial variable along the thickness, / is the porosity, and Asp is the specific catalyst area available for surface reactions. The ideal gas equation can then be used to calculate the pressure. The species molar flux Jk in the porous bed is evaluated using the DGM equation as

"K ! # Kg g X X DDGM ½X l  Bg DGM kl Jk ¼  Dkl r½X l  þ rp l Del;Kn l¼1 l¼1

ð14Þ

The DGM is written as an implicit relationship between the pressure gradient, concentration gradients, molar fluxes and molar concentrations. It neglects the effect of external forces and thermodiffusion [24,25]. The first term on the right-hand side of Eq. (14) represents the diffusive flux and the second term represents the viscous flux. DDGM is a matrix of DGM diffusion coefficients and is kl formulated as [22]

DDGM ¼ H1 kl

ð15Þ

where the elements of the H matrix are [22]

" hkl ¼

1

Dek;Kn

X Xj Xk e dkl þ ðdkl  1Þ e D D kj kl j–k

ð16Þ

Assuming the catalyst bed is made up of spherical particles, the permeability Bg in Eq. (14) is given by the Kozeny–Carman relationship [26] 2

Bg ¼

/3 dp

72sð1  /Þ2

ð17Þ

Here, dp is the particle diameter and s is the tortuosity. The effective Knudsen diffusion coefficient Dek;Kn in Eq. (14) is given by

Dek;Kn

/ dpore ¼ s 3

sffiffiffiffiffiffiffiffiffiffiffi 8RT pW k

ð18Þ

The effective binary diffusivity Dekl is given by

Dekl ¼

/

s

Dkl

qs C ps

  @T s @ @T s Ac Aact þ Q_  ¼ keff;s Pout @z @t @z As Vs

ð19Þ

The porous medium is considered to be a stationary component of the mixture in which the Chapman–Enskog kinetic theory is

ð22Þ

where t is the time, Ts is the temperature, qs is the density, Cps is the heat capacity, keff,s is the effective heat conductivity, Aact is the electrochemically active anode surface area, Vs is the volume of the solid, Pout is the power density output and Q_ is the heat source term arising from the interaction with the channels (this changes along the channel axis with species composition and by means of heat exchange with the solid-phase). Subscript ‘s’ refers to the solidphase or combined porous media. Here, Ac and As are the areas of cross section of the channel and solid-phase, respectively. As the channels experience negligible pressure change, Qs can be expressed in terms of an enthalpy flux.

Q_ ¼ 

#

þ

In this paper, the solid-phase refers to the combined thickness of the anode, electrolyte and cathode. The variation of temperature along the thickness is neglected because the length of the porous media is much greater than its thickness. Hence, the equation in 1-D is given as

Kg X @ðqf thk Y k Þ k¼1

@z

ð23Þ

In 1-D, effective solid thermal conductivity keff,s is obtained from the harmonic mean (decreases the influence of larger values and strongly tends towards smaller values in the element set) of the effective thermal conductivities of the three different porousmedia materials. The effective thermal conductivity of each material is obtained by averaging the sum of thermal conductivities over all radial nodes by the number of radial nodes, before the harmonic mean is taken. The effective conductivity at each radial nodal point of a specified material is given by [27,28]

keff ¼ /kfluid þ ð1  /Þksolid

ð24Þ

where / is the porosity, kfluid is the thermal conductivity of the gas and ksolid is the solid thermal conductivity. It is important to note that heat transfer to/from the interconnects is not modeled in this study and hence, surface to surface phenomena like radiation is not taken into account. Although radiation becomes important at high temperatures to accurately determine surface temperatures, we neglect it due to its significant influence being restricted only to the interconnects [1].

165

V. Menon et al. / Applied Energy 149 (2015) 161–175

Ecell ¼ Erev  ga ðiÞ  jgc ðiÞj  gohm ðiÞ  gconc ðiÞ

Table 1 Cell parameters and properties used for model validation/analysis.

where ga and gc are the activation overpotentials at the anode and cathode respectively, gohm is the ohmic overpotential, and gconc is the concentration overpotential. The concentration overpotential is not treated explicitly as porous media transport is modeled in detail, i.e., the reversible potential is calculated using gas-phase concentrations at the electrode–electrolyte interface. Erev is the ‘reversible’ cell voltage, which is the maximum possible potential that can be derived from a cell operating reversibly, and is given by the Nernst equation as

Model validation [19]

Model validation [20]

Parametric analysis (base case)

Gas channels (planar) Length (cm) Height (mm) Width (mm) Air inlet velocity (m/s) Fuel inlet velocity (m/s)

– – – – –

– – – – –

5, 7.5 1 1 2.5 0.3

Anode Thickness (lm) Porosity (%) Tortuosity Particle diameter (lm) Pore diameter (lm) Specific area (m1)

650 35 3.8 2.5 1.0 1.025  105

850 35 3.8 2.5 1.0 1.025  105

500 35 3.8 2.5 1.0 1.025  105

Electrolyte Thickness (lm)

75

85

25

Cathode Thickness (lm) Porosity (%) Tortuosity Particle diameter (lm) Pore diameter (lm) Specific area (m1)

35 35 3.8 2.5 1.0 1.025  105

50 35 3.8 2.5 1.0 1.025  105

30 35 3.8 2.5 1.0 1.025  105

Operating conditions Pressure (bar) Temperature (°C)

1.0 600, 700

1.0 550, 600, 650

1.0 650 (air channel), 800 (fuel channel)

3% H2O + 97% H2 3% H2O + 97% O2

3% H2O + 97% H2 3% H2O + 21% O2 + 76% N2

20.37% CH4 + 32.2% H2 + 13.67% CO + 0.87% CO2 + 32.89% H2O 3% H2O + 21% O2 + 76% N2

– – – – –

– – – – –

1.86 2.16 5.84 3515.75 452.63

Parameter

Inlet gas composition Fuel channel

Air channel

Thermal properties [1] kanode kelectrolyte kcathode Cps

qs

1=2

Erev

¼ 0;

¼0

ð29Þ

Rtot ¼ Rel þ Rcontact þ Rc þ Ra

ð30Þ

The magnitudes of these resistances depend on the type of material used and the micro-structure of the porous electrode. In modern cells, the electronic resistances of both electrodes Rc, Ra and the contact resistances between solid–solid interfaces Rcontact are negligible compared to the ionic resistance of the electrolyte Rel, which is given by

Rel ¼

lel

ð31Þ

rel

where lel is the thickness of the electrolyte, and rel is the electrolyte conductivity, with the SI unit – S/m, which varies as a strong function of temperature,

rel ¼ r0 T 1 exp

  Eel RT

ð32Þ

Here, Eel is the activation energy for ion transport (29.5 kJ/mol) and

r0 is the pre-exponential factor (339.34  102 S/m) [19]. The activa-

Similarly, for the solid-phase temperature (Eq. (22)) adiabatic boundary conditions are used at both ends, i.e., at z = 0 and z = L,

z¼0

ð28Þ

where Rtot is given by

ð25Þ

 dT s  dz 

!

gohm ¼ Rtot i

The inlet boundary conditions at time t = 0 are given in Table 1. At the exit of the reactor channels, adiabatic Neumann boundary conditions are used,

 dT s  dz 

pH2 ;a pO2 ;c RT ¼E þ ln 2F pH2 O;c 0

where E0 is the electromotive force (EMF) at standard pressure, and pi represents the partial pressures of O2, H2O (at the cathode TPB), and H2 (at the anode TPB). The temperature dependent E0 is calculated from thermodynamic data (DG(T)/2F). The source of the thermodynamic species data is the JANAF tables. The ohmic overpotential in Eq. (27) is given by

2.4. Boundary conditions

 dT f  ¼0 dz z¼L

ð27Þ

ð26Þ

z¼L

2.5. Electrochemistry The charge transfer chemistry occurs at the three-phases boundaries (TPB), which are basically interfaces formed by the electrocatalyst, electrolyte and gas-phase boundaries. In this study, we only consider charge transfer occurring at the electrode–electrolyte interface (interfacial charge transfer) and not across the utilization region of the electrodes (distributed charge transfer). The potential balance equation relates the cell voltage to all overpotential losses that occur during operation, and is given as

tion overpotential at the electrode–electrolyte interfaces is related to the current density, in implicit form, by the non-linear B–V equation. For the electrochemical oxidation of H2, the modified B–V equation is used (detailed derivation is given in Appendix A). The multi-step reactions that are considered for the derivation of the modified B–V equation are [29,30]

H2 ðgÞ þ 2ðNiÞ 2HðNiÞ

ð33Þ

HðNiÞ þ ðelÞ Hþ ðelÞ þ e ðNiÞ þ ðNiÞ

ð34Þ

Eq. (34) is assumed to be rate limiting. Consequently, the modified B–V equation is given by

     Fne Fn i ¼ i0;a exp ba ga  exp ð1  ba Þ e ga RT RT

ð35Þ

where i is the current, i0 is the exchange current density, ne is the number of electrons involved per reaction (equal to 1 since a single electron transfer reaction is considered), ga is the anode activation overpotential, F is the Faraday constant, T is the local MEA temperature and b is the symmetry factor [22]. At OCV, the rate of charge transfer in the anodic direction is equal to the rate of charge transfer in the cathodic direction, and is equal to the exchange current density. For the electrochemical reduction of O2, a global B–V equation

166

V. Menon et al. / Applied Energy 149 (2015) 161–175

is used because the kinetic parameters needed to calculate the equilibrium partial pressure of O2 were not available in the literature for the specific material being considered. This takes the form

     Fne Fn i ¼ i0;c exp ba gc  exp ð1  ba Þ e gc RT RT

ð36Þ

where ne is the number of electrons involved per reaction (equal to 2) and gc is the cathode activation overpotential [31]. The exchange current density is expressed as a function of temperature, and partial pressures of products and reactants (obtained by coupling the electrochemical model with the micro-kinetic model) participating in the charge transfer chemistry, although it makes more physical sense to express it as a function of open surface coverage and surface coverage of electrochemically active species. The exchange current density is given by,

  ðpH2 =pH2 Þ1=4 EH i0;a ¼ kH2 exp  2 RT 1 þ ðpH =pH Þ1=2 2

ð37Þ

2

at the anode, and

   pO2 ;c 1=4 EO i0;c ¼ kO2 exp  2 RT pamb

ð38Þ

at the cathode. Here kH2 , EH2 , kO2 and EO2 are the electrochemical fit parameters, and pH2 is the equilibrium partial pressure. The order dependency of the oxygen partial pressure is obtained by assuming surface diffusion of oxygen intermediates to be the rate limiting step, rather than the charge transfer reaction [32–34]. However, the effect of foreign impurities in the TPB and degradation effects on the Ni catalyst, which are neglected here, should also be considered for a more robust model. 2.6. Thermo-catalytic chemistry In this analysis, we use an elementary surface reaction mechanism that consists of 42 reactions with 6 gas-phase species. The detailed multi-step heterogeneous reaction mechanism, used in this study, takes into account the adsorption/desorption of H2, O2, CH4, CO, CO2 and H2O from the surface of Ni [35]. It also encompasses the water–gas shift reaction, formation of carbon monolayer, methanation reactions, steam reforming, dry reforming, and partial and total oxidation of C1 species. It was made thermodynamically consistent and extended to temperatures between 493.15 and 1973.15 K. Furthermore, it was validated and tested with experimental data obtained from methane reforming over nickel/alumina monoliths in the temperature range of 900–1350 K, and with additional data from literature [36]. Comparison with equilibrium calculations demands further work to predict surface carbon deposition rates. For the fuel channel composition considered in this paper, gas-phase reactions in the fuel channel can be neglected [37]. The computational framework used in modeling these reactions via the mean-field approximation can be found elsewhere [1,9,10,22,38]. 2.7. Computational scheme Eqs. (4)–(6), (8), (12), (13), (22) and (27) form a system of coupled non-linear equations, which can be treated as a differential– algebraic system mathematically. Their residual form can be written as

FðUÞ ¼ 0

ð39Þ

where the vector U is given by

_ TÞfc ; ðY; q; T; hÞa;1 . . . ðY; q; T; hÞa;n ; ðTÞel ; ðY; q; TÞc;1 ; . . . U ¼ ½ðY; m; _ TÞac T ðY; q; TÞc;m ; ðY; m;

ð40Þ

Here, the indices (fc), (a, 1), (a, n), (el), (c, 1), (c, m), (ac) stand for fuel channel, first discretized cell in the anode, nth discretized cell in the anode, cell in the electrolyte, first discretized cell in the cathode, mth discretized cell in the cathode, and air channel, respectively. They are first cast in finite-volume form, with the specification of required number of axial and transverse nodes/cells. The entire solution procedure follows a space marching algorithm: at each axial position the transient system of equations is solved until a steady-state solution is obtained. The initial condition at each axial position assumes the converged solution from the previous finite volume cell. The solid-phase temperature is decoupled from that of the fluid-phase, due to the fact that heat transfer has a higher time constant compared to mass transfer, diffusion and kinetics. At every time step of the solid-phase temperature calculation, space marching of the entire channel is carried out. At every axial position, the solid phase temperature is obtained, and subsequently used to calculate the reaction rates and gas phase temperature. The equation system is solved using the differential algebraic equation (DAE) solver LIMEX [39]. A damped Newton iteration algorithm is employed to solve the system of algebraic model equations – Eqs. (27), (35) and (36), to obtain the current density [1,9,10,22,38]. The entire program is written in FORTRAN and is a part of the software package DETCHEM™ [21]. The flow chart illustrating the implementation of the solution algorithm is shown in Fig. 2. 3. Results and discussion A comparison between the simulated polarization curves and the first set of experimental data obtained from literature [19] is depicted in Fig. 3a. In the experiments, the electrochemical characteristics of a single SOFC with a dense Sm-doped BaCeO3 electrolyte (BaCe0.8Sm0.2O2.90 – BCSO) was studied. The anode and cathode were made of porous NiO–BCSO and Ba0.5Sr0.5Co0.8Fe0.2O3–BCSO (BSCF–BCSO), respectively. Conductivities of BCSO at 773.15 K, 873.15 K and 973.15 K were 0.416 S/m, 0.662 S/m and 0.936 S/m, respectively. These conductivity values correspond to ohmic resistances of 1.132 X cm2 and 0.799 X cm2 for a 75 lm BCSO film at 873.15 K and 973.15 K, respectively. The cell was operated with humidified hydrogen (3% H2O) as fuel and pure oxygen as oxidant. Button cell simulations are used to reproduce experimental data and derive the electrochemical fit parameters, which are the preexponential factors and activation energies. At low current densities, the non-linearity of the V–I curves increases with decreasing temperature due to the dominance of activation losses. The OCV calculated by the model is higher than what is predicted by the experiments. One has to bear in mind that the OCV calculated in the code, using the Nernst potential equation, only predicts the maximum theoretical electromotive force when no external current is flowing. This differs from the measured OCV obtained in the experiments due to i. Partial contribution from electron–hole charge carriers/defects via the progressive reduction of Ce4+ ions with increasing temperatures. These BaCeO3-based ceramics are known to exhibit hole conduction under oxidizing atmospheres, as pure oxygen is fed into the cathode chamber in this case [40]. ii. Apart from electron–hole conduction, Samarium-doped Barium Cerate (BCS) is a hybrid oxide-ion and proton conductor [40]. In order to obtain the right OCV, one should take into account the oxide-ion-conductivity and proton-conductivity via the measured ionic transport numbers of the oxide-ion and proton, respectively. From Ref. [40], it is seen that the proton transport number decreases with increasing temperature and dopant concentration, while the oxide-ion transport

167

V. Menon et al. / Applied Energy 149 (2015) 161–175

Start t=0:me

Inial condions (Table 1) For N=0:Naxialcells (ll residuals < 10-6)

Jk

.

Yelectrolyte-electrode Ts

Ts

At TPB, solve electrochemical model (Eqs. 27, 35 & 36) i, Ecell, ηa

1-D channel conservaon (Eqs. 4, 5 & 8)

Porous media transport (Eqs. 12, 13 & 14) For N=0:Ntransversecells t=t+∆t

Q, Pout

Solid temperature ‘Ts‘ calculaon (Eq. 22)

Recalculate Ts

Ychannel-electrode

(ll residuals < 10-6)

Ts No

Check if t > me Yes Stop

Fig. 2. Schematic representation of the overall solution algorithm using a flowchart.

0.3

δanode = 650 μm δelectrolyte = 75 μm δcathode = 35 μm

Cell voltage (V)

1 0.8

0.2

700 °C

0.6 0.4

0.1

0.2

600 °C

0

Power denisty (W/cm2)

1.2

0 0

0.2

0.4

0.6

0.8

Current density (A/cm2) Fig. 3a. Model validation with experimental data obtained from Ranran et al. [19].

number increases with rising temperature and dopant content for BaCe1xSmxO3a materials. Interested readers are directed to Ref. [41]. iii. The densification of the electrolyte by Ni-addition can affect the OCVs. An increase in Ni content decreases the OCV due to an increase in electronic conductivity. This causes a reduction in the ionic transport numbers, which in turn reduces the ionic conductivities [42]. iv. The usage of the Nernst equation, which considers the partial pressure of steam at the cathode TPB in the denominator of the logarithmic term (Eq. (28)). For practical simulation of the numerical model, a small amount of steam needs to be fed at the inlet of the air channel along with the gas composition that is actually used in the experiments. For the model simulations, the inlet composition at the air channel was 3% H2O + 97% O2, while the experiments used pure O2. Nevertheless, this initial guess was found to have a very minor affect on the OCV.

For our model validation, we assume proton conduction to be predominant below 1027 K, i.e., proton transport number tH = 1 [40]. To obtain OCV values corresponding to those that are experimentally measured, one has to characterize the electrolyte to attain the ionic transport numbers along with its electronic, protonic and oxide-ionic conductivities and incorporate these values in OCV calculation. Nevertheless, our model predicts the theoretical maximum OCV [43]. While it is right that the OCV and concentration losses are coupled, it is important to remember that BaCeO3-based ceramics are mixed-ionic electronic conductors. At conditions where the contributions of all charge carrying defects become important, the net current density should be calculated as the sum of the individual oxide-ionic, protonic and electronic current densities. In such cases, steam is produced at both the anode–electrolyte and cathode–electrolyte interfaces, thereby affecting the mass transfer within the electrodes. This influences the concentration overpotentials, even if the proton transport number tH is close to 1 for our model conditions. Moreover, this effect becomes more pronounced at high current densities. Hence, future work could involve the consideration of an equivalent circuit modeling these individual current densities with their respective potential balance equations involving Nernst potentials and irreversibilities. Nevertheless, this phenomenon was also observed by Ni et al. [13]. Hence, the model is in reasonable agreement with experiments. The values of thicknesses of the electrodes and electrolyte, along with cell properties used in model validation are listed in Table 1. The electrochemical parameters used for reproducing the experimental data are shown in Table 2. The model is also validated with a second set of experiments from literature [20]. The anode, electrolyte and cathode were composed of Ni–BCY10, yttrium-doped barium cerate (BCY10 – BaCe0.9Y0.1O2.95) and Pr2NiO4+d, respectively. The materials were chosen because of their good electrochemical and morphological properties, although electrode–electrolyte interfaces

V. Menon et al. / Applied Energy 149 (2015) 161–175

0.12

0.95

0.09

0.8

0.06

650 °C 0.65

0.03

550 °C

0.5 0

0.05

600 °C

0.1

0.15

Power density (W/cm2)

Cell voltage (V)

0.15

δanode = 850 μm δelectrolyte = 85 μm δcathode = 50 μm

1.1

0 0.25

0.2

Current density (A/cm2) Fig. 3b. Model validation with experimental data obtained from Taillades et al. [20].

Table 3 Electrochemical model/input parameters (Fig. 3b).

a

Anode symmetry factor (ba) Cathode symmetry factor (ba)

0.5 0.5

Exchange current density parameters Pre-exponential for H2 oxidation (kH2 ) (A/cm2) Pre-exponential for O2 reduction (kO2 ) (A/cm2) Activation energy for H2 oxidation (EH2 ) (J/mol) Activation energy for O2 reduction (EO2 ) (J/mol)

288176.65a 528582.89a 100.0  103 120.0  103

Fitted with experimental data.

1080

Inlet velocity: Fuel channel - 0.3 m/s Air channel - 2.5 m/s

Flow

Temperature (K)

suffered from delamination post-operation. Ohmic resistances of the cell at 823.15 K, 873.15 K and 923.15 K were 2.38 X cm2, 1.84 X cm2 and 1.36 X cm2, respectively. The model is able to reproduce the experiments well, as shown in Fig. 3b. Although 823.15 K seems to be low from a SOFC standpoint, it was considered because the conductivity of the proton-conducting electrolyte at this temperature was 0.357 S/m. This is on par with its oxide-ion-conducting counterparts, and is important due to its contribution to the ohmic overpotential. Also, the conductivity of the proton-conductor, used in the experiments of Fig. 3a [19], was found to be 0.416 S/m at 773.15 K and is higher than the conductivity of the oxide-ion-conductor used in Refs. [1,10,22,38]. The validated electrochemical model parameters are described in Table 3. The pre-exponential factors in Tables 2 and 3 are used as fit parameters, while the activation energies are taken from literature [31]. It is important to note that the micro-structural properties are assumed to be within a realistic range, due to its unavailability in the aforementioned literature. In order to understand the multi-physics in a methane-fed proton conducting SOFC, simulations are carried out using the electrochemical model parameters listed in Table 2. All other relevant information for the parametric analysis is given in Table 1. The fuel is assumed to enter the fuel channel at 1073.15 K and 0.3 m/s, while air enters the air channel at 973.15 K and 2.5 m/s. These conditions represent the base case. Fig. 4 shows the temperature profiles of the fuel channel, solid (MEA) and the air channel, at t ? 500 s (steady-state distribution). The temperature profiles result from the combination of heterogeneous chemical reactions in the anode, heat transfer between the channels and the solid, exothermic electrochemical reactions at the TPB, and ohmic heat produced by the resistances of the MEA materials. Close to the inlet section, the temperature in the fuel channel drops due to heat lost to the air channel, as air enters at a lower temperature. This leads to a rise in temperature in the air channel. As one moves towards the exit of the gas channels, the temperature starts to increase due to heat released within the MEA. At steady-state (t > 500 s), the temperatures in the gas channels approach the temperature of the solid. The multi-physics phenomena affecting temperature profiles in the cell will be discussed hereafter. The variations in fuel channel temperature with length and cell voltage are shown in Fig. 5. All temperature profiles are shown at steady-state. When operated at OCV, the fuel channel continuously loses heat to the air channel by convection due to higher air inlet velocities. The fuel channel attains an exit temperature of 941.4 K, due to the occurrence of the exothermic water–gas shift reaction. When a voltage of Ecell = 0.6 V is applied, there is a rise in fuel channel temperature after its initial drop. This can be attributed to the exothermic electrochemical reactions occurring at the TPB, which also produces ohmic heat. The current density output at Ecell = 0.9 V is lower than that at Ecell = 0.6 V. This results in a slower exothermic electrochemical reaction rate, thereby causing a temperature rise that is smaller in magnitude as one moves towards the exit of the fuel channel. The net current density output of the cell at Ecell = 0.9 V is 0.056 A/cm2, while that at Ecell = 0.6 V is

1040

Solid (porous MEA) 1000

Fuel channel Ecell = 0.6 V Tfc, in - 1073.15 K Tac, in - 923.15 K

960

Air channel 920 0

0.01

0.02

0.03

0.04

0.05

Axial posion (m) Fig. 4. Temperature profiles within the fuel channel, air channel, and solid (MEA) along the length of the cell.

Inlet velocity: Fuel channel - 0.3 m/s Air channel - 2.5 m/s

1120

Temperature (K)

168

1080

z = 5 cm

Flow

1040

Ecell = 0.6 V 1000

OCV

960

Ecell = 0.9 V

Tfc, in - 1073.15 K Tac, in - 923.15 K

920 0

0.015

0.03

0.045

0.06

0.075

Axial posion (m) Table 2 Electrochemical model/input parameters (Fig. 3a).

a

Anode symmetry factor (ba) Cathode symmetry factor (ba)

0.5 0.5

Exchange current density parameters Pre-exponential for H2 oxidation (kH2 ) (A/cm2) Pre-exponential for O2 reduction (kO2 ) (A/cm2) Activation energy for H2 oxidation (EH2 ) (J/mol) Activation energy for O2 reduction (EO2 ) (J/mol)

174919.71a 58016.6a 100.0  103 120.0  103

Fitted with experimental data.

Fig. 5. The effect of cell length and applied voltage on the temperature profile of the fuel channel. The dashed line represents the temperature profile for a channel length of 7.5 cm.

0.428 A/cm2. An increase in the length of the cell causes further increase in temperature due to higher electrochemical utilization of available hydrogen in the fuel channel. The effect of specific catalytic area on fuel channel temperature is illustrated in Fig. 6. The characteristic initial temperature drop in

169

V. Menon et al. / Applied Energy 149 (2015) 161–175

0.7

1.06

isothermal 0.6

1.04

isothermal 0.5

1.02

0.4

1

non-isothermal

0.98

0.2 0

0.01

2.5625x105 m-1

1.025x105 m-1

1000

960

Ecell = 0.6 V Tfc, in - 1073.15 K Tac, in - 923.15 K

4.1x105 m-1

0.03

0.04

0.05

Fig. 7. Reversible potential and current density for the cell operating at 0.6 V. The temperature for the isothermal case is 1073.15 K.

decline in the concentration of H2 occurs, when a cell voltage of 0.6 V is applied, due to its electrochemical oxidation. We speculate that the initial rise in CH4 axial concentration is due to the hydrogenation of surface carbon C(Ni). The decrease in CH4 axial concentration after 3.25 cm is due to weak steam reforming that leads to a slightly higher concentration of CO2 than CO (due to the WGS reaction) towards the exit. Predicting accurate surface carbon coverage is vital for the correct determination of OCVs [38]. The equilibrium composition (computed using DETCHEMEQUIL) at the exit of the fuel channel, at 1044.96 K, was found to be 1% CH4 + 66.96% H2 + 21.32% CO + 2.99% CO2 + 7.73% H2O. It serves as a means to understand if the catalytic activity is optimized to give maximum possible selectivity and/or yield. In specific, this shows slow reaction kinetics over the considered catalyst surface, which can be improved by increasing the specific catalytic area and residence time of the gases in the fuel channel. Also, simulations with a simple plug flow model (computed using DETCHEMPLUG) employing the same specific catalytic area and operating conditions, at 941.4 K (OCV exit temperature – Fig. 5), yielded a similar outlet gas composition as Fig. 9 when the cell is operated at OCV. The coverage of C(Ni) adsorbate is slightly higher at 0.6 V than at OCV, as shown in Fig. 10. The key global reaction that governs this phenomena involves the hydrogenation of CO (CO + H2 M C + H2O). Also, it is important to bear in mind that the flux of H2 towards the anode TPB rises with increasing current densities, while the amount of H2O decreases due to its electrochemical formation on the cathode side. This promotes the hydrogenation of surface carbon by shifting its equilibrium towards product (CH4) formation. Moreover, the above findings are specific to the inlet fuel composition used in this study. The gas concentrations of various species and surface coverages of major species within the anode at three different positions along

Ecell = 0.6 V

ηc

0.27

0.14

ηohm

0.25

ηohm

ηc

0.23 0.21

0.12 0.1

isothermal non-isothermal

ηa

0.19

0.08

ηa

0.17

0.06

0.15

920 0

0.01

0.02

0.03

0.04

0.05

Axial posion (m) Fig. 6. The effect of specific catalytic area on the temperature profile of the fuel channel.

0.16

Overpotenal (V)

1040

0.02

Axial posion (m)

Overpotenal (V)

Inlet velocity: Fuel channel - 0.3 m/s Air channel - 2.5 m/s

0.3

Current density (A/cm2)

0.8

Ecell = 0.6 V

non-isothermal

0.29

1080

Temperature (K)

1.08

Reversible potenal (V)

the inlet region appears to be similar in magnitude for all the three specific areas. However, the temperature is lower for a catalyst with higher specific catalytic area, as one moves towards the exit of the fuel channel. In the first 1.5 cm of the channel, H2 and CO2 production occurs through the exothermic water–gas shift reaction that leads to the consumption of H2O and CO. But, the heat loss to the air channel dominates and hence, governs the early drop in temperature. After 1.5 cm, CH4 and H2O are consumed by the endothermic steam reforming process, as a result of which concentrations of H2 and CO continue to increase. Nonetheless, exothermic electrochemical reactions dominate and govern the overall temperature rise, though higher specific areas can lead to more reforming, thereby decreasing the exit temperature. The current density and reversible potential profiles, for both isothermal and non-isothermal cases, are depicted in Fig. 7. The isothermal temperature is taken to be equal to the temperature at the fuel channel inlet, i.e., 1073.15 K. In the non-isothermal case, the reversible potential increases close to the inlet due to the temperature drop in the solid. The combined effect of increasing solid temperature as well as fuel dilution causes the reversible cell potential to decrease towards the exit. However, in the isothermal case, the linear drop in reversible cell potential is only due to fuel dilution. The current density initially increases near the inlet, then decreases due to the temperature drop in the solid and continues to increase due to the rise in solid temperature, when operating non-isothermally. Current density seems to be a weak function of fuel dilution but a strong function of solid temperature. This phenomenon is similar to what happens in an oxide-ion-conducting SOFC when operated with humidified methane [1]. In the isothermal case, the current density simply varies as a function of hydrogen concentration in the fuel channel. Sufficient hydrogen is produced in the case considered here. An increase in the length of the fuel channel would further bring down the current density, due to enhanced hydrogen utilization. The axial variation of overpotential losses along the cell length is shown in Fig. 8, for both isothermal and non-isothermal operation. The cathode activation overpotential is higher than the anode activation and ohmic overpotentials. As expected, the ohmic overpotential varies as a function of current density and temperature. The axial species distribution resulting from the combination of thermo-catalytic chemical reactions and electrochemical fluxes in the anode is shown in Fig. 9, for the cell operating at OCV and 0.6 V. At OCV, species profiles proceed pre-dominantly through the equilibrium between the WGS/RWGS reactions. Due to the WGS reaction, the consumption of H2O and CO, and the production of H2 and CO2 follow similar trends. The decrease in axial CH4 concentration can be attributed to weak steam reforming that leads to the added production of H2 and CO. Subsequently, the CO generates more CO2 via the WGS reaction, as depicted by the axial profiles. A

0.04 0

0.01

0.02

0.03

0.04

0.05

Axial posion (m) Fig. 8. Overpotential loses as a function of axial position along the cell. The temperature for the isothermal case is 1073.15 K.

170

V. Menon et al. / Applied Energy 149 (2015) 161–175

OCV 0.6 V 0.22

0.5

Flow 0.215

H2

0.21

0.3

CH4

0.2

0.205

H2O

xCH4

Mole fracon

0.4

0.2

CO 0.1

0.195

CO2

0 0

0.01

0.19 0.02

0.03

0.04

0.05

Axial posion (m) Fig. 9. Gas-phase species distribution in the fuel channel along the length of the cell.

Anode depth (µm)

(a)

At Ecell = 0.6 V,

450

Anode / Channel interface

400

Coverage: C(Ni) 4.4E-05 4.2E-05 4E-05 3.8E-05 3.6E-05 3.4E-05 3.2E-05 3E-05 2.8E-05

350 300 250 200 150 100 50

Electrolyte / Anode interface 0

0.01

0.02

0.03

0.04

0.05

Axial length (m)

At OCV,

Anode depth (µm)

(b) 450

Anode / Channel interface

400

Coverage: C(Ni) 4.25E-05 4.2E-05 4.15E-05 4.1E-05 4.05E-05 4E-05 3.95E-05 3.9E-05 3.85E-05

350 300 250 200 150 100 50

Electrolyte / Anode interface 0

0.01

0.02

0.03

0.04

0.05

Axial length (m) Fig. 10. The influence of cell voltage on surface carbon coverage within the anode when operated at: (a) Ecell = 0.6 V, and (b) OCV.

the cell length are shown in Figs. 11–13. Fig. 11 describes this close to the inlet of the cell. The major adsorbates that determine the open surface coverage of the Ni catalyst are CO(Ni) and H(Ni). The trend in these coverages is similar to the concentration of its corresponding gas species in the anode. At Ecell = 0.6 V, the decrease in concentration of H2 towards the TPB due to electrochemical consumption is offset by its production from the WGS reaction. This is reflected by decreasing H2O and CO, and increasing CO2, concentrations towards the TPB away from the anode–channel interface. An interesting phenomenon to note in Figs. 12 and

13 is the increase in the concentration of adsorbate H(Ni) to a point where it exceeds the concentration of adsorbate CO(Ni), as one moves down the channel length. This can be attributed to current density which also drives the direction of the H2 gas within the anode. One can note a more prominent decrease in H2 gas concentration towards the TPB, with increasing axial length, due to an increasing current density profile along the axial length (Fig. 7). However, this is not the case at OCV due to the absence of electrochemical flux. The amount of catalyst covered increases as we move towards the exit, although only marginally at OCV, due to higher surface coverage of adsorbate CO(Ni) and increasing surface coverage of adsorbate H(Ni). The gas composition of the air channel as a function of axial distance is described in Fig. 14. The mole fractions of O2 and N2 decrease, while that of H2O increases due to electrochemical reactions at the cathode TPB (Eq. (2)). The corresponding pattern is reflected in the cathode gas concentrations at three different axial locations of the air channel. The profiles at OCV act as a yardstick to measure deviation in species distribution when a voltage is applied. In order to understand the effect of catalyst zoning on temperature profiles in the cell, the anode is divided into four equal zones with increasing and decreasing specific catalytic areas as presented in Fig. 15. Case I represents a step decrease in specific area, while case II represents a step increase in specific area. The difference in exit temperatures between cases I and II is 18.5 K. In the initial 5.1 cm of the cell, case I has a lower temperature due to higher amount of methane consumed in the endothermic steam reforming process as compared to case II. After 5.1 cm, case II has a lower temperature/higher temperature drop due to higher endothermic steam reforming rates as opposed to case I. Thus, the amount of H2 and CO at the exit of the fuel channel is greater in case II, i.e.,

xH2;case I xH2;case II

¼ 0:782 and

xCOcase I xCOcase II

¼ 0:841. However, the net

output current density for both cases do not differ significantly from each other (0.02 A/cm2) due to the lack of significant fuel depletion, and temperatures in both sections (before and after 5.1 cm) balancing each other out. The effect of different operating conditions on exit temperature of the solid (at steady-state) and electrochemical performance parameters are listed in Table 4. The base case (Table 1) is compared with three other cases: (1) inlet temperature at the air channel Tac,in = 873.15 K, (2) inlet temperature at the fuel channel Tfc,in = 1223.15 K, and (3) inlet velocity at the fuel channel tfc,in = 2.5 m/s. A predictable trend was observed for all the cases. A decrease in the inlet temperature of the gases at the air channel leads to a higher temperature drop in the porous MEA, which causes the attainment of a lower current density due to prominent overpotential losses. This is accompanied by low hydrogen and oxygen utilization in the gas channels. The power density of the cell can be improved by raising the inlet gas temperature at the fuel channel, although this would include the supply of more heat energy to the system and steepen the temperature gradients. As expected, an increase in current density raises the electrochemical flux that leads to higher consumption of active species participating in charge transfer. Highest power density was obtained for case (3), which only involved raising the inlet velocity of gases at the fuel channel. This results in the highest solid temperature amongst the four cases. The utilization of electrochemically active H2 is the least due to the dominance of the WGS reaction and higher supply of fuel at the inlet. Nevertheless, the oxygen consumption is the highest. Although the ohmic resistance of the electrolyte decreases with an increase in temperature, its overpotential is driven strongly by rising current densities. In comparison, the impact of MEA materials on H-SOFC performance, obtained by other research groups, is listed in Table 5. A good compilation of this can be found in Lin et al. [45].

171

V. Menon et al. / Applied Energy 149 (2015) 161–175

Electrolyte / Anode interface 0

Anode / Channel interface

Anode thickness (μm) 100

200

300

400

500

0.45

0.45

z = 5.36 mm 0.4

CO(Ni) H(Ni)

0.35

0.35

0.3

0.3

Coverage

Coverage

0.4

(Ni) 0.25

0.25

0.2 0.42

0.2 0.22

H2

0.17

H2O CH4

0.22

0.12

CO 0.12

Mole fracon

Mole fracon

0.32

OCV 0.6 V

0.07

CO2

0.02

0.02 0

Electrolyte / Anode interface

100

200

300

400

500

Anode / Channel interface

Anode thickness (μm)

Fig. 11. Species profiles and surface coverages within the anode near the reactor inlet (z = 5.36 mm). Top panel shows the surface coverage and bottom panel shows species profiles.

Electrolyte / Anode interface 0

Anode / Channel interface

Anode thickness (μm) 100

200

300

400

500

0.45

0.45

z = 2.54 cm 0.4

0.4

0.35

0.35

H(Ni)

0.3

0.3

(Ni)

0.25

Coverage

Coverage

CO(Ni)

0.25

0.2 0.42

OCV 0.6 V

0.2 0.16

Mole fracon

0.32

H2O CH4

0.22

0.11

CO2

Mole fracon

H2

0.12

CO

0.02 0

Electrolyte / Anode interface

100

200

300

Anode thickness (μm)

400

0.06 500

Anode / Channel interface

Fig. 12. Species profiles and surface coverages within the anode halfway from the inlet (z = 2.54 cm). Top panel shows the surface coverage and bottom panel shows gaseous species profiles.

172

V. Menon et al. / Applied Energy 149 (2015) 161–175

Electrolyte / Anode interface

Anode / Channel interface

Anode thickness (μm)

0

100

200

300

400

500 0.46

0.46

Coverage

0.41

H(Ni)

0.41

CO(Ni)

0.36

0.36 0.31

0.31

0.26

0.26

Coverage

z = 4.54 cm

(Ni) 0.21

0.21

0.16 0.42

0.16 0.21

H2O

0.16

CH4

0.22

0.11

CO

0.12

Mole fracon

H2

0.32

Mole fracon

OCV 0.6 V

CO2

0.02

0.06 0

Electrolyte / Anode interface

100

200

300

400

500

Anode / Channel interface

Anode thickness (μm)

Fig. 13. Species profiles and surface coverages within the anode near the reactor exit (z = 4.54 cm). Top panel shows the surface coverage and bottom panel shows gaseous species profiles.

xH2O

xH2O

25

0

0.03

0

0.04

0.08

O2 H2O

H2O

20

O2

15 10

O2

5 0 0.205 0.209 0.213

H2O 0.196 0.205 0.214

0.185

xO2

xO2

z = 5.36 mm

0.2

0.215

OCV 0.6 V

xO2

z = 2.54 cm

0.75

z = 4.54 cm

0.07

N2

Flow Mole fracon

xH2O 0.06

0.06

0.6 0.05 0.45

0.04

xH2O

Cathode depth (μm)

0.027 0.031 0.035 30

H2O

0.3

0.03

O2 0.15

0.02 0

0.01

0.02

0.03

0.04

0.05

Axial posion (m) Fig. 14. Gas-phase species distribution within the air channel and cathode, as a function of cell length. Top panels depict the gas concentrations within the cathode, while the bottom panel shows the gaseous concentrations in the air channel.

173

V. Menon et al. / Applied Energy 149 (2015) 161–175

Case I

Case II

1140

Flow

Temperature (K)

4.0E+05

Asp (m-1)

3.0E+05 2.0E+05 1.0E+05 0.0E+00

z = 7.5cm

1100

z = 5.62 cm

z = 1.87 cm 1060

z = 3.75 cm

1020

Case II

980

Case I

940

0.0-0.018

0.018-0.037 0.037-0.056 0.056-0.075

0

0.015

Axial posion (m)

0.03

0.045

0.06

0.075

Axial posion (m)

Fig. 15. The effect of catalyst zoning on the temperature distribution in the fuel channel. The left panel shows two cases: (i) case I represents a step-wise decrease in Asp and, (ii) case II represents a step-wise increase in Asp. The right panel shows its impact on the temperature profiles in the fuel channel.

Table 4 Effect of operating conditions on system performance parameters. Case

i (A cm2)

Erev (V)

Pout (W cm2)

ga (V)

gc (V)

gohm (V)

Ts (K)

xH2% (at exit)

xO2% (at exit)

Base (1) (2) (3)

0.4285 0.2858 0.4681 0.5593

1.043 1.065 1.037 1.029

0.2571 0.1715 0.2809 0.3356

0.073 0.0927 0.0685 0.0608

0.2611 0.2868 0.2534 0.2402

0.1105 0.0869 0.116 0.1286

1043.75 970.34 1059.32 1079.4

28.85 29.15 27.46 34.8

18.98 19.71 18.79 18.38

Table 5 Impact of MEA on peak power density for H-SOFCs. Ref.

MEA T (°C) (anode/electrolyte/cathode)

Hibino et al. [44]

3 wt% Pd-loaded FeO/ BCY25/Ba0.5Pr0.5CoO3a Ni–BCY10/BCY10/BSCFb

Lin et al. [45] Dailly et al. [46]

Ni–BCY10/BCY10/ Pr2NiO4+dc

Fuel Ppeak (mW/cm2)

600

H2

134

600 700 600

H2

380 550 180

H2

a

BCY25–25 mol% Y3+-doped BCY (0.5 mm thick). Oxidant-air. Electrolyte (BCY10) – BaCe0.9Y0.1O2.95 (50 lm thick), (BSCF) – Ba0.5Sr0.5Co0.8Fe0.2O3d. Oxidant-air. c Electrolyte – 50 lm thick. Oxidant-compressed air. b

cathode

parameters constant. It was determined that a higher inlet velocity at the fuel channel results in higher current and power densities but lower fuel utilization rates and efficiencies. Better heat management configurations and innovative materials are important in enhancing performance and sustainability of the system. Future work could involve the qualitative comparison between the co-flow, counter-flow and cross-flow configurations to study heat generation/utilization in the cell. The electrochemical model used in this paper does not consider defect transport within the proton-conducting electrolyte. A detailed algorithm that describes the transport of multiple charge-carrying defects can be found in Kee et al. [47]. Nevertheless, the numerical model presented in this paper can be considered as one of the first steps towards understanding the interplay between various physico-chemical phenomena in a proton-conducting SOFC.

4. Conclusion A numerical investigation of proton-conducting SOFCs with direct internal reforming has been carried out. The model is fitted against two sets of experimental data to derive electrochemical parameters. Further parametric analysis is carried out to study the influence of operating conditions on mass and heat transport in the cell. It was found that the WGS reaction played an important role in driving temperature distributions in the cell. The air channel was mainly used to regulate temperature gradients, and for heat management purposes. The equilibrium between the WGS/RWGS reactions was found to dominate the occurrence of weak endothermic steam reforming. Since the electrochemical formation of H2O occurs at the cathode side, further steam reforming of CH4 is not possible unlike its oxide-ion-conducting counterpart. Moreover, slow chemical kinetics were exhibited under the used operating conditions. However, it should be noted that the findings were limited to the specific fuel composition used for this analyses. The effects of specific catalytic area and cell voltage on temperature and species distribution in the system were also investigated. A qualitative study of the influence of cell voltage on surface carbon coverage shows its concentration to be the highest in the first half of the anode. Two different configurations of catalyst zoning in the anode does not result in a significant difference in temperature gradients, although it could be used to enhance cell efficiencies and fuel utilization rates. A brief analysis of the impact of fuel and air inlet conditions on the cell was conducted by keeping all other

Acknowledgements We thank the Steinbeis GmbH für Technologietransfer (STZ 240 Reaktive Strömungen) for a cost free academic license of DETCHEM™. We deeply value all the rewarding exchanges on this topic with Prof. Vinod M. Janardhanan (Indian Institute of Technology Hyderabad). Financial support by the Helmholtz Research School Energy-Related Catalysis is gratefully acknowledged. Appendix A The appendix serves to provide details about the derivation of the modified Butler–Volmer equation from the following elementary reaction steps considered for the electrochemical oxidation of H2 at the anode–electrolyte three phase boundary (TPB),

H2 ðgÞ þ 2ðNiÞ 2HðNiÞ

ðA:1Þ

HðNiÞ þ ðelÞ Hþ ðelÞ þ e ðNiÞ þ ðNiÞ

ðA:2Þ

Reaction (A.1) describes the dissociative adsorption of gas-phase hydrogen, H2(g), into two empty active sites on the Ni surface, (Ni). For the charge-transfer step, as illustrated in reaction (A.2), the H2 spillover mechanism is assumed in accordance with the general consensus in the literature. The adsorbed atomic hydrogen, H(Ni), vacates its Ni site, (Ni), and spills over to the electrolyte to

174

V. Menon et al. / Applied Energy 149 (2015) 161–175

release a proton in the electrolyte, H+(el), and an electron in the Ni anode, e(Ni). Assuming the charge-transfer step to be rate limiting, the fast adsorption–desorption reaction proceeds to equilibrium. Therefore, from the law of mass action,

K1 ¼

h2H

ðA:3Þ

pH2 h2Ni

where K1 is the equilibrium constant, hH and hNi are the site fractions or surface coverages of adsorbed hydrogen and empty active sites on the Ni surface and pH2 is the partial pressure of hydrogen at the TPB. For the charge-transfer reaction, the rate is dependent on both the surface coverages of the reacting species and the potential difference between the electrode and the electrolyte, Ea. Therefore, the current density, i, for the rate-limiting step is given by

     b2;a FEa b2;c FEa  k2;c hNi exp i ¼ lTPB F k2;a hH exp RT RT

ðA:4Þ

where lTPB is the TPB length, F is the Faraday constant, k2,a and k2,c are the reaction rate constants in the anodic and cathodic direction respectively, b2,a and b2,c are the anode and cathode charge transfer coefficients, R is the universal gas constant and T is the cell temperature. It is assumed here that there is no change in the site vacancies in the electrolyte, (el), and H+ ion concentrations at the TPB. Consequently, the surface coverages in the electrolyte are neglected in Eq. (A.4). Noting that b2,a + b2,c = 1, at equilibrium, when the forward reaction rate equals the backward reaction rate (i.e. i = 0), the equilibrium potential difference, Eeq a , is given by

exp



 FEeq a RT

¼

Here, the parameter pH2 is obtained from a balance between adsorption and desorption rates of H on Ni at equilibrium,

c0

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   RT pH2 2 Edes 2 2 hNi ¼ Ades exp hH C 2pW H2 RT RT

which yields

pH2 ¼ pH2

ðA:5Þ

 2   hNi Edes C2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ Ades exp 2pW H2 RT hH RT c0

ðA:14Þ

where C = 2.6  109 mol/cm2, the pre-exponential factor Ades = 5.59  1019 cm2 mol1 s1, sticking coefficient for H2 adsorption c0 = 0.01, the activation energy Edes = 88.12 kJ/mol, and W H2 is the molecular weight of H2. Finally substituting Eqs. (A.11) and (A.12) in Eq. (A.9), and grouping the leading constants in Eq. (A.9) into a single parameter, the exchange current density can be expressed as 

i0 ¼ iH2

ðpH2 =pH2 Þð1b2;a Þ=2

ðA:15Þ

1 þ ðpH2 =pH2 Þ1=2

  l Fk  The first term, iH2 ¼ TPBK b 2;a , is a constant for a given electrode 2

microstructure and is a function of temperature only. It is adjusted to fit the measured fuel cell performance. The second term which shows how the exchange current density depends on the partial pressures of the reactant, divulges the apparent reaction order of the charge transfer reaction. For b2,a = 0.5, Eq. (A.15) reduces to 

i0 ¼ iH2

1 hNi K 2 hH

ðA:13Þ

ðpH2 =pH2 Þ1=4 1 þ ðpH2 =pH2 Þ1=2

ðA:16Þ

k

where the equilibrium constant K 2 ¼ k2;a . Substituting Eq. (A.5) in 2;c

References

Eq. (A.4) leads to the Butler–Volmer form for current density,

     b2;a F ga b2;c F ga  exp i ¼ i0 exp RT RT

ðA:6Þ

where the exchange current density,

i0 ¼ lTPB Fk2;a hH exp

  b2;a FEeq a RT

ðA:7Þ

and the anodic activation overpotential is defined as

ga ¼ Ea  Eeq a

ðA:8Þ

The expression for exchange current density can be further simplified by raising Eq. (A.5) on both sides by b2,a and substituting the result in Eq. (A.7):

i0 ¼

lTPB Fk2;a K b2

1b2;a b2;a hNi

hH

ðA:9Þ

Now, by definition, the site fractions on the Ni anode surface must sum up to unity:

hH þ hNi ¼ 1

ðA:10Þ

Substituting Eq. (A.10) in Eq. (A.3) and defining K 1 ¼ 1=pH2 , hH and hNi can be expressed in terms of the gaseous partial pressures of hydrogen,

hH ¼

ðpH2 =pH2 Þ1=2 1 þ ðpH2 =pH2 Þ1=2

ðA:11Þ

and,

hNi ¼

1 1 þ ðpH2 =pH2 Þ1=2

ðA:12Þ

[1] Janardhanan VM, Deutschmann O. Numerical study of mass and heat transport in solid-oxide fuel cells running on humidified methane. Chem Eng Sci 2007;62:5473–86. [2] Steele BCH. Fuel-cell technology: running on natural gas. Nature 1999;400: 619–21. [3] Finnerty CM, Coe NJ, Cunningham RH, Ormerod RM. Carbon formation on and deactivation of nickel-based/zirconia anodes in solid oxide fuel cells running on methane. Catal Today 1998;46:137–45. [4] Hajimolana SA, Hussain MA, Daud WMAW, Soroush M, Shamiri A. Mathematical modeling of solid oxide fuel cells: a review. Renew Sustain Energy Rev 2011;15:1893–917. [5] Bhattacharyya D, Rengaswamy R. A review of solid oxide fuel cell (SOFC) dynamic models. Indust Eng Chem Res 2009;48:6068–86. [6] Janardhanan VM, Deutschmann O. Modeling of solid-oxide fuel cells. Zeitschrift fur Physikalische Chemie 2007;221:443–79. [7] Menon V, Janardhanan VM, Tischer S, Deutschmann O. A novel approach to model the transient behavior of solid-oxide fuel cell stacks. J Power Sources 2012;214:227–38. [8] Laguna-Bercero MA. Recent advances in high temperature electrolysis using solid oxide fuel cells: a review. J Power Sources 2012;203:4–16. [9] Menon V, Janardhanan VM, Deutschmann O. A mathematical model to analyze solid oxide electrolyzer cells (SOECs) for hydrogen production. Chem Eng Sci 2014;110:83–93. [10] Menon V, Fu Q, Janardhanan VM, Deutschmann O. A model-based understanding of solid-oxide electrolysis cells (SOECs) for syngas production by H2O/CO2 co-electrolysis. J Power Sources 2015;274:768–81. [11] Demin AK, Tsiakaras PE, Sobyanin VA, Hramova SY. Thermodynamic analysis of a methane fed SOFC system based on a protonic conductor. Solid State Ionics 2002;152–153:555–60. [12] Demin A, Tsiakaras P. Thermodynamic analysis of a hydrogen fed solid oxide fuel cell based on a proton conductor. Int J Hydrogen Energy 2001;26:1103–8. [13] Ni M, Leung MKH, Leung DYC. Mathematical modelling of proton-conducting solid oxide fuel cells and comparison with oxygen-ion-conducting counterpart. Fuel Cells 2007;7:269–78. [14] Ni M, Leung DYC, Leung MKH. Modeling of methane fed solid oxide fuel cells: comparison between proton conducting electrolyte and oxygen ion conducting electrolyte. J Power Sources 2008;183:133–42. [15] Ni M. Modeling of a planar solid oxide fuel cell based on proton-conducting electrolyte. Int J Energy Res 2010;34:1027–41.

V. Menon et al. / Applied Energy 149 (2015) 161–175 [16] Patcharavorachot Y, Brandon NP, Paengjuntuek W, Assabumrungrat S, Arpornwichanop A. Analysis of planar solid oxide fuel cells based on protonconducting electrolyte. Solid State Ionics 2010;181:1568–76. [17] Arpornwichanop A, Patcharavorachot Y, Assabumrungrat S. Analysis of a proton-conducting SOFC with direct internal reforming. Chem Eng Sci 2010;65:581–9. [18] Ni M, Leung DYC, Leung MKH. Electrochemical modeling of ammonia-fed solid oxide fuel cells based on proton conducting electrolyte. J Power Sources 2008;183:687–92. [19] Ranran P, Yan W, Lizhai Y, Zongqiang M. Electrochemical properties of intermediate-temperature SOFCs based on proton conducting Sm-doped BaCeO3 electrolyte thin film. Solid State Ionics 2006;177:389–93. [20] Taillades G, Dailly J, Taillades-Jacquin M, Mauvy F, Essouhmi A, Marrony M, et al. Intermediate temperature anode-supported fuel cell based on BaCe0.9Y0.1O3 electrolyte with novel Pr2NiO4 cathode. Fuel Cells 2010;10:166–73. [21] Deutschmann O, Tischer S, Kleditzsch S, Janardhanan VM, Correa C, Chatterjee D, et al. DETCHEMTM software package, 2.4 ed. Karlsruhe; 2012. [22] Zhu H, Kee RJ, Janardhanan VM, Deutschmann O, Goodwin DG. Modeling elementary heterogeneous chemistry and electrochemistry in solid-oxide fuel cells. J Electrochem Soc 2005;152:A2427–40. [23] Hayes RE, Kolaczkowski ST. Introduction to catalytic combustion. 1st ed. Amsterdam: Gordon and Breach Science Publishers; 1997. [24] Jackson R. Transport in porous catalysts. New York: Elsevier Scientific Publishing Co.; 1977. [25] Mason EA, Malinauskas AP. Gas transport in porous media: the dusty-gas model. New York: Elsevier; 1983. [26] Bear J. Dynamics of fluids in porous media. New York: American Elsevier; 1972. [27] Beale SB. Transport phenomena in fuel cells: numerical models for planar solid oxide fuel cells. Massachusetts: WIT Press; 2005. [28] VanderSteen JDJ, Kenney B, Pharoah JG, Karan K. Mathematical modeling of the transport phenomena and the chemical/electrochemical reactions in solid oxide fuel cells: a review. In: Proceedings of hydrogen and fuel cells. Toronto, Canada; 2004. [29] Vogler M. Elementary kinetic modelling applied to solid oxide fuel cell pattern anodes and a direct flame fuel cell system. Ph.D. Dissertation. Heidelberg: Universität Heidelberg; 2009. [30] Bessler WG, Warnatz J, Goodwin DG. The influence of equilibrium potential on the hydrogen oxidation kinetics of SOFC anodes. Solid State Ionics 2007;177:3371–83. [31] Costamagna P, Selimovic A, Del Borghi M, Agnew G. Electrochemical model of the integrated planar solid oxide fuel cell (IP-SOFC). Chem Eng J 2004;102:61–9. [32] He F, Wu T, Peng R, Xia C. Cathode reaction models and performance analysis of Sm0.5Sr0.5CoO3d–BaCe0.8Sm0.2O3d composite cathode for solid oxide fuel cells with proton conducting electrolyte. J Power Sources 2009;194:263–8.

175

[33] Zhao L, He B, Gu J, Liu F, Chu X, Xia C. Reaction model for cathodes cooperated with oxygen-ion conductors for solid oxide fuel cells using proton-conducting electrolytes. Int J Hydrogen Energy 2012;37:548–54. [34] Peng R, Wu T, Liu W, Liu X, Meng G. Cathode processes and materials for solid oxide fuel cells with proton conductors as electrolytes. J Mater Chem 2010;20:6218–25. [35] Maier L, Schädel B, Herrera Delgado K, Tischer S, Deutschmann O. Steam reforming of methane over nickel: development of a multi-step surface reaction mechanism. Top Catal 2011;54:845–58. [36] Hecht ES, Gupta GK, Zhu H, Dean AM, Kee RJ, Maier L, et al. Methane reforming kinetics within a Ni–YSZ SOFC anode support. Appl Catal A: Gen 2005;295: 40–51. [37] Sheng CY, Dean AM. Importance of gas-phase kinetics within the anode channel of a solid-oxide fuel cell. J Phys Chem A 2004;108:3772–83. [38] Janardhanan VM, Deutschmann O. CFD analysis of a solid oxide fuel cell with internal reforming: coupled interactions of transport, heterogeneous catalysis and electrochemical processes. J Power Sources 2006;162:1192–202. [39] Deuflhard P, Hairer E, Zugck J. One-step and extrapolation methods for differential–algebraic systems. Numer Math 1987;51:501–16. [40] Iwahara H, Yajima T, Hibino T, Ushida H. Performance of solid oxide fuel cell using proton and oxide ion mixed conductors based on BaCe1xSmxO3a. J Electrochem Soc 1993;140:1687–91. [41] Huang J, Yuan J, Mao Z, Sundén B. Analysis and modeling of novel lowtemperature SOFC with a co-ionic conducting ceria-based composite electrolyte. J Fuel Cell Sci Technol 2009;7:011012. [42] Shimada H, Li X, Hagiwara A, Ihara M. Proton-conducting solid oxide fuel cells with yttrium-doped barium zirconate for direct methane operation. J Electrochem Soc 2013;160:F597–607. [43] Chiodelli G, Malavasi L. Electrochemical open circuit voltage (OCV) characterization of SOFC materials. Ionics 2013;19:1135–44. [44] Hibino T, Hashimoto A, Suzuki M, Sano M. A solid oxide fuel cell using Y-doped BaCeO3 with Pd-loaded FeO anode and Ba0.5Pr0.5CoO3 cathode at low temperatures. J Electrochem Soc 2002;149:A1503–8. [45] Lin Y, Ran R, Zheng Y, Shao Z, Jin W, Xu N, et al. Evaluation of Ba0.5Sr0.5Co0.8Fe0.2O3d as a potential cathode for an anode-supported proton-conducting solid-oxide fuel cell. J Power Sources 2008;180:15–22. [46] Dailly J, Marrony M, Taillades G, Taillades-Jacquin M, Grimaud A, Mauvy F, et al. Evaluation of proton conducting BCY10-based anode supported cells by co-pressing method: up-scaling, performances and durability. J Power Sources 2014;255:302–7. [47] Kee RJ, Zhu H, Hildenbrand BW, Vøllestad E, Sanders MD, O’Hayre RP. Modeling the steady-state and transient response of polarized and nonpolarized proton-conducting doped-perovskite membranes. J Electrochem Soc 2013;160:F290–300.