Modeling and Bifurcation Analysis of a Coionic Conducting Solid

Article pubs.acs.org/IECR

Modeling and Bifurcation Analysis of a Coionic Conducting Solid Oxide Fuel Cell Mona Bavarian,† Ioannis G. Kevrekidis,‡ Jay B. Benziger,‡ and Masoud Soroush†,§ †

Department of Chemical and Biological Engineering, Drexel University, Philadelphia, Pennsylvania 19104 Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544

Downloaded via UNIV OF SOUTH DAKOTA on August 21, 2018 at 08:45:24 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

ABSTRACT: A macroscopic first-principles mathematical model of a solid oxide fuel cell (SOFC) with a BaCe1−xSmxO3‑α type electrolyte is developed. This class of electrolytes exhibits both proton and oxygen-anion conductivity. The existence of steadystate multiplicity with respect to bifurcation parameters such as cell outlet voltage, cell current, ohmic external load, and cell power density is investigated. The analyses with respect to the first two bifurcation parameters represent potentiostatic and galvanostatic modes of operation, respectively. The cell can have up to three steady states with respect to the external load resistance and cell outlet voltage, and a unique steady state with respect to the cell current. The cell can have four steady states with respect to cell power density (when power density is the bifurcation parameter). The same qualitative steady-state behavior had been observed in oxygen-ion conducting and proton conducting SOFCs. Furthermore, thermal and concentration multiplicities coexist in this class of SOFCs; ignition in the solid temperature is accompanied by (a) extinction in the fuel and oxygen concentrations, and (b) ignition and extinction in concentrations of water in the anode and cathode sides, respectively. The cell exhibits steady-state multiplicity with respect to the fuel and air flow rates under fixed ohmic load, potentiostatic, and fixed cell power density modes of operation, but does not show this behavior under a galvanostatic mode of operation. The dynamic response of the cell also shows that the cell is highly nonlinear and some of the cell variables can show inverse response.

1. INTRODUCTION In typical solid oxide fuel cells the electrolyte is an oxygen ion conductor, commonly yttria-stabilized zirconia (YSZ).1,2 During the last two decades, several oxides have been found with potential application as the electrolytes in SOFCs.3 Among these electrolytes are the ones that have high proton conductivity at lower temperatures or have mixed proton and oxygen ion conductivity under certain conditions. It has been shown that SOFCs with proton conducting electrolytes can have higher efficiencies compared to oxygen ion conducting ones.3,4 Among various perovskites, barium cerates are potential candidates for use in SOFCs. BaCeO3 by itself has a low conductivity. When it is doped with a trivalent rare-earth element such as Gd, Nd, and Sm, a higher conductivity can be achieved. Substitution of Ce4+ with a trivalent rare-earth element causes the oxygen vacancy (VÖ ), which is responsible for oxygen ion conductivity, to appear in the lattice to retain the electrical neutrality.5 The oxygen vacancy can combine with O2 and H2O, leading to the formation of the lattice oxygen (Oxo), an electron hole (h·) and a proton (H+) according to the following postulated reaction mechanism:4−6 VÖ +

BaCeO3-base electrolytes exhibit proton conduction under hydrogen-containing atmospheres and are primarily considered as proton conducting ceramics. However, an increase in temperature can change these proton conductors to oxygenanion conductors.7 BaCe1−xSmxO3−α electrolytes have been observed to behave as a mixed proton and oxygen ion conductor under some fuel cell conditions, as water formation was detected in both cathode and anode gas channels of the cell. As can be seen in Figure 1,8 at temperatures above 1000 K, the proton conductivity of BaCe1−xSmxO3−α decreases but its oxygen-anion conductivity increases as temperature increases. Above 1073 K, its oxygen-anion conductivity is higher than its proton conductivity. Its oxygen-anion conductivity changes more with temperature than its proton conductivity does. To predict the behavior of a fuel cell system with a BaCeO3-base electrolyte and determine the rate of electrochemical reactions, both oxygen ion and proton conduction should be taken into account. To achieve this, protonic and oxide-ionic currents are calculated from the rates of water production at the cathode and anode sides of the cell, respectively. There have been several studies9−12 on electrochemical transport and performance behavior of “mixed conducting” oxides for utilization in fuel cells, batteries, and sensors. Herein, oxides that have ionic and electronic conductivity (allow one mobile ionic species [e.g., oxygen ion] and electronic carriers [electrons and holes] to go through) are called mixed

1 O2 → Oox + 2h· 2

VÖ + H 2O → Oox + 2H+ Special Issue: Process Engineering of Energy Systems

VÖ + H 2O → Oox → 2OH·o

h·+H 2O → 2H+ +

Received: Revised: Accepted: Published:

1 O2 2 © 2012 American Chemical Society

3165

April 16, 2012 August 24, 2012 September 18, 2012 September 19, 2012 dx.doi.org/10.1021/ie3009814 | Ind. Eng. Chem. Res. 2013, 52, 3165−3177

Industrial & Engineering Chemistry Research

Article

presents and examines the dynamic response of the system. Finally, concluding remarks are given in section 7.

2. MODEL DEVELOPMENT The SOFC system under study is a single planar SOFC. Figure 2 shows the schematic of the SOFC. As can be seen, the cell is

Figure 1. Component conductivities in BaCe1−xSmxO3−α8. Reprinted with permission from ref 8. Copyright 1993 The Electrochemical Society.

conducting oxides, and oxides that conduct two different ionic species and electronic charge carriers are referred to as coionic conducting oxides. Compared to the studies on fuel cells with mixed conducting oxides,9−12 there have been only a very few theoretical studies on electrochemical transport and performance behavior of coionic conducting oxides in solid-state electrochemical devices. For example, Huang et al.13 presented a mathematical model to describe the performance of a coionic conducting ceria-based composite electrolyte fuel cell. Their focus was more on transport of ionic species and electronic carriers, and calculation of the ionic and electronic flux by utilizing the intrinsic ionic and electronic transport properties of the oxide layer. They did not account for mass transfer in the gas channels or heat transfer between the solid part of the cell and the gas streams in the cathode and anode gas channels. Comprehensive reviews of recent publications on mathematical modeling of SOFCs can be found in refs 2,14 and 15. The number of studies on mathematical modeling of SOFCs with alternative types of oxides (proton conducting or coionic conducting) is very limited. As many researches explore the utilization of these alternative oxides for application in SOFCs, there is also a need to theoretically model and capture the behavior of these fuel cells. Prior to this work, no study on steady-state multiplicity in coionic-conducting SOFCs was reported. This paper presents a study on mathematical modeling and steady-state analysis of a coionic-conducting solid oxide fuel cell with a Sm-doped barium cerate electrolyte, Ni as the anode and Ni0.95Li0.05O1‑α as the cathode. The model is validated using the experimental data reported in ref 8, wherein a hydrogen−air SOFC with BaCe0.9Sm0.1O3‑α as electrolyte, Ni as the anode, and various oxides as the cathode was studied experimentally. The existence of multiple steady states with respect to ohmic external load, cell outlet voltage, cell current, and cell power density is studied. Detailed steady-state analysis of this type of cells is indispensible to determine whether the cells can exhibit steady-state multiplicity. The dynamic response of the cell in two different operating regions is also examined, which confirms the multitime scale and highly nonlinear nature of the cell. Organization of the rest of this paper is as follows. Section 2 describes the model development. Section 3 presents the model parameter estimation and validation results. Section 4 explains the simulation method used to generate the results. Section 5 discusses the steady state behavior of the system. Section 6

Figure 2. Planar SOFC divided into five subsystems.

divided into five subsystems: (1) anode gas channel, (2) anode side diffusion layer, (3) solid part of the cell consisting of the anode, electrolyte and cathode, (4) cathode side diffusion layer, and (5) cathode gas channel. The following assumptions are made: (i) Each subsystem is locally homogeneous; each subsystem is considered as a stirred tank. (ii) The electron hole conductivity is ignored compared to the ionic conductivities (proton and oxygen ion), as the electronic conductivity is one order of magnitude less than those of oxygen ion and proton. (iii) Heat transfer by radiation between the electrodes and the gases in the anode and cathode side channels is ignored. For nonpolar gases, the radiation heat transfer can be ignored since the gases do not emit radiation.16 Three of the species, O2, N2 and H2 are nonradiating, while water vapor is radiating. Due to the low partial pressure of water vapor, the emissivity factor of water vapor is low, and accordingly heat transfer by radiation is negligible in comparison to heat transfer by convection. (iv) All gases are considered to be ideal gases. (v) The fuel cell is adiabatic. (vi) The inlet and outlet velocities of the gas stream in each channel are equal. 2.1. Electrochemical Submodel. As illustrated in Figure 3, the electrochemical reactions occur in the triple-phase

Figure 3. Triple phase boundary in a coionic conducting SOFC. 3166

dx.doi.org/10.1021/ie3009814 | Ind. Eng. Chem. Res. 2013, 52, 3165−3177


Article

The actual cell voltage, V, which is lower than both reversible cell voltages (electric potentials) due to the corresponding activation, ohmic and concentration overpotentials, is calculated from

boundary zone where the oxygen-anions/protons, electrons, and gas phase are in contact. The hydrogen oxidation and oxygen reduction result in the formation of protons and oxygen ions, respectively. Unlike in SOFCs with an oxygen ion-conducting or a proton-conducting electrolyte in which water is produced either in the anode or in the cathode gas channel, in coionic conducting SOFCs water is formed in both anode and cathode gas channels. The electrochemical reactions in the anode and cathode sides of a coionic conducting SOFC are as follows: (1)

H 2 + O2 − → H 2O + 2e−

(2)

Cathode side: 1 O2 + 2H+ + 2e− → H 2O 2

(7)

V = VOC,H − ηPol,H

(8)

where ηPol,O and ηPol,H are, respectively, total overpotentials for the oxygen-ion and proton conducting subcells, given by

Anode side: H 2 → 2H+ + 2e−

V = VOC,O − ηPol,o

ηPol,O = ηan,O + ηcath,O + ρE,O dE

IO + ηconc,O LB

(9)

ηPol,H = ηan, H + ηcath,H + ρE,H dE

IH + ηconc,H LB

(10)

In Figure 4, RO,PolIO = ηPol,O and RH,PolIH = ηPol,H.The terms ηan,H and ηan,O represent the anode activation polarizations related to the reactions in eqs 1 and 2, respectively, and ηcath,H and ηcath,O are the cathode activation polarizations related to the reactions in eqs 3 and 4, respectively. The activation overpotentials are calculated from the following Butler−Volmer equations considered for the anode and cathode sides:18,19

(3)

1 O2 + 2e− → O2 − (4) 2 These reactions reflect the postulation that (a) hydrogen molecules release electrons at the anode and the resulting hydrogen cations (protons) leave the anode, go through the electrolyte and arrive at the cathode side to receive electrons and react with oxygen molecules to form water molecules, and (b) oxygen molecules receive electrons at the cathode, and the resulting oxygen ions move in the reverse direction and go to the anode side to release electrons and react with hydrogen molecules to form water molecules. On the basis of our second assumption that the electron hole conductivity is negligible compared to the ionic conductivities (proton and oxygen ion) and the postulation made above, we propose the equivalent circuit shown in Figure 4 for the cell.

⎛ − E A ⎞⎧ ⎛ A F ⎞ IO = γ ApHTPB pHTPB,an ηan,O⎟ exp ⎜ ⎟⎨exp⎜θa O 2 2 LB ⎝ RTs ⎠⎩ ⎝ RTs ⎠ ⎪

⎪

⎛ ⎞⎫ F − exp⎜ −θcA ηan,O⎟⎬ RTs ⎝ ⎠⎭

⎪

⎪

(11)

⎛ − E C ⎞⎧ ⎛ C F ⎞ IO 0.25 = γ CpOTPB exp⎜ ηcath,O⎟ ⎟⎨exp⎜θa 2 LB ⎝ RTs ⎠⎩ ⎝ RTs ⎠ ⎪

⎪

⎛ ⎞⎫ F − exp⎜ −θcC ηcath,O⎟⎬ RTs ⎝ ⎠⎭ ⎪

⎪

(12)

⎧ ⎛ F ⎞ ⎛ ⎞⎫ IH F = i0,an ⎨exp⎜α ηan,H⎟ − exp⎜ − (1 − α) ηan,H⎟⎬ LB RTs ⎠ ⎝ ⎠⎭ ⎩ ⎝ RTs ⎪

⎪

⎪

⎪

(13)

⎧ ⎛ F ⎞ IH = i0,cath ⎨exp⎜α ηcath,H⎟ LB ⎠ ⎩ ⎝ RTs ⎪

Figure 4. Equivalent circuit of the fuel cell.

⎪

⎛ ⎞⎫ F − exp⎜ −(1 − α) ηcath,H⎟⎬ RTs ⎝ ⎠⎭ ⎪

The circuit includes an electromotive force due to the oxygen ion current serially connected to the corresponding polarization “resistance” (charge transfer, ohmic and diffusion resistances), and these two elements are connected in parallel to two similar elements related to the proton current. The reversible cell voltages, VOC,O and VOC,H, are calculated using the Nernst equation: ⎛ ⎞ bulk,an RTS ⎜ pH2O ⎟ VOC,O = E − ln nF ⎜⎜ p bulk p bulk 0.5 ⎟⎟ ⎝ H2 O2 ⎠

(5)

⎛ bulk,cath ⎞ RTS ⎜ pH2O ⎟ =E − ln nF ⎜⎜ p bulk p bulk 0.5 ⎟⎟ ⎠ ⎝ H2 O2

(6)

⎪

where θ is the charge transfer coefficient and α is the symmetry factor, whose value is an indicator of the symmetry of the activation barrier.20 The terms ρE,OdEIO/(LB) and ρE,HdEIH/(LB) in eqs 9 and 10 represent the ohmic overpotentials, which are potential losses due to the resistance of the electrolyte with respect to oxygen ion and proton flows, where ρE,O (ρE,H) is the oxygen ion (proton) electrolyte resistivity, dE is the electrolyte thickness, and IO (IH) is the oxygen ion (proton) current. As shown in Figure 1, the electrolyte resistivity with respect to oxygen ions decreases as temperature increases. To describe the oxygen-ion resistivity, we curve-fitted an Arrhenius type equation to the experimental data in ref 8, leading to the correlation:

0

VOC,H

(14)

0

3167


Industrial & Engineering Chemistry Research ρE,O =

Article

⎛ 8266 ⎞ 1 exp⎜ ⎟ 3834 ⎝ TS ⎠

R H2O,an = (15)

The proton resistivity decreases slightly as temperature decreases in the range of 700−1000 °C. Using the experimental data in ref 8, and curve-fitting, we arrived at the following correlation for the proton resistivity: 1 ρE,H = 3701 4437 0.002894 exp T − 0.001201 exp T

( ) S

S

The terms ηconc,O and ηconc,H in eq 9 and 10 represent the concentration overpotentials (losses), which account for the potential losses due to the difference between the concentrations of each species in the bulk and the triple phase boundary. This loss is the difference between the Nernst potentials in the bulk and the triple phase boundary:17

⎛ 1 ⎞ ⎜ ⎟I ⎝ 2F ⎠ H

⎜c cathy cath k Ocath O2 2 ⎜ air ⎝

an an k Han2⎜c fuel yH ⎜ 2

⎝

NH2O,an =

pOTPB ⎞ − 2 ⎟ RTs ⎟⎠

(26)

pHTPB ⎞ 2 ⎟ − RTs ⎟⎠

(27)

an an k Han2O⎜ −c fuel yH O ⎜ 2

⎝

⎞ pHTPB,an O + 2 ⎟ RTs ⎟⎠

(28)

(29)

where kj is the mass transfer coefficient of species j, given by k Han2 = k Hj 2O =

k Ocath = 2

De H2

=

Δan

De H2O Δj

DeO2 Δcath

εan DH τanΔan 2 =

=

εj τjΔj

D H 2O ,

εcath DO τcathΔcath 2

(30)

j = cath, an (31)

(32)

Dej is the effective diffusion coefficient of species j, εcath and τcath are the porosity and tortuosity of the cathode, respectively, and εan and τan are the porosity and tortuosity of the anode electrode. Dj is the total diffusion coefficient of species j that accounts for both molecular and Knudsen diffusion mechanisms:21,22 1 1 1 = + Dj Dj , i /mix D Kj (33) The Knudsen diffusion, DKj, depends on the pore mean radius, cell temperature, and species molecular weight. It is calculated from23

(21)

DK j = 97δj

where RH2 and RO2 are the rates of consumption of hydrogen and oxygen by the electrochemical reactions, respectively, and RH2O is the rate of generation of water, obtained from Faraday’s law:

R H2O,cath =

(25)

⎛

where j = cath, an. ⎛ TPB ⎞ Δan d ⎜ pH2 ⎟ 1 = NH2 − RH R dt ⎜⎝ Ts ⎟⎠ LB 2

⎛ 1 ⎞ ⎜ ⎟I ⎝ 4F ⎠

⎛ p TPB,cath ⎞ ⎜ −c cathy cath + H2O ⎟ NH2O,cath = k Hcath air H 2O 2O⎜ RTs ⎟⎠ ⎝

At low currents, the activation overpotentials are dominant. As the cell current increases, ohmic losses (mainly because of the electrolyte resistance rather than the resistances of the electrodes) lower the cell potential. Lastly at high current densities, mass-transport is the dominant cause in lowering the cell potential; at high currents, the concentrations of the reactants and products in the electrode−electrolyte interface are significantly different from those in the bulk. 2.2. Mass-Transfer Submodel. Mole balances on hydrogen, water, and oxygen in each subsystem of the cell yield this submodel. Mole balances on (a) hydrogen and water in the anode-side diffusion layer and (b) oxygen and water in the cathode-side diffusion layer lead to

(20)

R O2 =

NH2 =

(18)

TPB, j Δj d ⎛ pH2O ⎞ 1 ⎜ ⎟ = −N R H O, j H 2O, j + R dt ⎜⎝ Ts ⎟⎠ LB 2

(24)

⎛

⎛ bulk bulk 0.5 ⎞ ⎛ TPB TPB 0.5 ⎞ RTs ⎜ pH2 pO2 RTs ⎜ pH2 pO2 ⎟ ⎟ ln⎜ bulk,cath ⎟ − ln⎜ TPB,cath ⎟ = ⎟ ⎜ ⎟ nF ⎜ pH O nF p ⎠ ⎝ H 2O ⎝ ⎠ 2

(19)

⎛ 1 ⎞ ⎜ ⎟I ⎝ 2F ⎠

NO2 =

⎛ bulk bulk 0.5 ⎞ ⎛ TPB TPB 0.5 ⎞ RTs ⎜ pH2 pO2 RTs ⎜ pH2 pO2 ⎟ ⎟ ln = ⎟⎟ − nF ln⎜⎜ p TPB,an ⎟⎟ nF ⎜⎜ pHbulk,an H 2O ⎝ ⎠ ⎝ ⎠ 2O

⎛ TPB ⎞ Δcath d ⎜ pO2 ⎟ 1 = NO2 − RO R dt ⎜⎝ Ts ⎟⎠ LB 2

R H2 =

⎛

(17)

ηconc,H

(23)

where the total current, I = IO + IH. Nj is the molar flux of species j transferred from the gas channel into the diffusion layer or vice versa. The mass flux of the diffusing species can be defined in terms of the mass transfer coefficient using

( )

(16)

ηconc,O

⎛ 1 ⎞ ⎜ ⎟I ⎝ 2F ⎠ O

Ts Mj

(34)

where j = cath, an. Dj,i/mix is the binary diffusion coefficient of species j in i or in the case of multiple components in the system is the diffusion coefficient of species j in the mixture. In this study, we have three species, O2, N2, and H2O in the cathode gas stream, so the molecular diffusion is the ternary diffusion and is given by the Wilke formula:24,25

(22) 3168


Industrial & Engineering Chemistry Research Dj ,mix =

1 − xj x n ∑i = 1, i ≠ j Di j,i

Article

Table 1. Fuel Cell Model Parameter Values (35)

The binary diffusion coefficient Dj,i is calculated using the Chapman and Enskog equation:26 Dj , i =

0.00266TS1.5 pMji 0.5σji 2 Ω Dji

(36) 2

where Dj,i is the diffusion coefficient (cm /s), p is pressure (bar), σji is the characteristic length (Å), σji = (σj+σj)/2, and ⎛ 1 ⎞−1 1 ⎟ Mji = 2⎜⎜ + Mi ⎟⎠ ⎝ Mj

(37)

Mj and Mi are molecular weights of j and i species. ΩDji is the diffusion collision integral of j and i species, given by26 Ω Dji =

AΩ CΩ EΩ + + B Ω * exp(DΩTij ) exp(FΩTij*) (Tij*) +

GΩ exp(HΩTij*)

(38)

where Tij*= kT/εji, εji = (εjεi) , k is Boltzmann’s constant, and ε is the characteristic Lennard-Jones energy. The values of the parameters in eq 38 are given in Table 1. Mole balances on (a) hydrogen and water in the anode gas channel and (b) oxygen and water in the cathode gas channel yield 1/2

LBW

cath dcair cath cath = nair,in cair WB − No2(LB) ̇ cath − vair dt

+ NH2O,cath(LB) LBW

cath cath d(cair yo ) 2

dt

= nair,in ̇ cath yo

LBW

(39)

2,in

cath cath cath − vair cair yo WB − No2(LB) 2

cath cath d(cair yH O ) 2

dt

(40)

cath cath cath = −vair cair yH O WB + NH2O,cath(LB) 2

(41)

LBW

an dc fuel an an = n fuel,in c fuelWB − NH2(LB) ̇ an − vfuel dt

+ NH2O,an(LB) LBW

an an d(c fuel yH ) 2

dt

=

n fuel,in ̇ an yH 2,in

(42)

−

− NH2(LB)

an an an vfuel c fuelyH WB 2

(43)

2.3. Heat-Transfer Submodel. Energy balances for the anode and cathode gas channels as well as the solid part of the cell are written to describe the temperatures of the subsystems. For the cathode and anode gas channels, the enthalpies of the inlet fuel and air streams, heat transfer by convection between the gas channels and the solid part, and the enthalpies of outlet fuel and air streams are accounted for. The governing equations for heat transfer in the subsystems SS1 and SS5 (shown in Figure 2) are as follows:

parameter

value

AΩ B BΩ CΩ DΩ EΩ EA EC FΩ GΩ HΩ L Tcath air,in Tan fuel,in W dA dC dE i0,an i0,cath ncath air,in nan fuel,in α ρsĈ cp γA γC θAa θAc θCa θCc Δan Δcath εcath εan τan τcath δan δcath υH2

1.06036 0.005 m 0.15610 0.19300 0.47635 1.03587 253648 J mol−1 198954 J mol−1 1.52996 1.76474 3.89411 0.01 m 890 K 890 K 5 × 10−3 m 5 × 10−5 m 5 × 10−5 m 5 × 10−4 m 1330 A m−2 1330 A m−2 5.6 × 10−6 mol s−1 5.6 × 10−6 mol s−1 0.5 106 J m−3K−1 3.7 × 107 A m−2 1.5 × 109 A m−2 2 1 1.4 0.6 5 × 10−6 m 5 × 10−6 m 4 × 10−1 4 × 10−1 5 5 5 × 10−7 5 × 10−7 7.07

26 this 26 26 26 26 this this 26 26 26 this this this this 19 19 8 this this 33 33 19 29 this this 18 18 18 18 this this 19 19 19 19 23 23 21

υO2

16.6

21

υN2

17.9

21

υH2O

12.7

21

σH2

2.968

34

σO2

3.433

34

σN2O

2.65

34

σN2

3.681

34

yO2,in

0.2

this work

yN2,in

0.8

this work

yH2,in

0.97

this work

yH2O,in

0.03

this work

(LBW ) =

source

work work

work work work work

work work

work work

work work

an an Cp,fuelTan) d(c fuel

dt

an n fuel,in ̇ an Hfuel,in

an an an c fuelHfuel (BW ) + han(Ts − Tan) − vfuel

(LB) − NH2HHan2(LB) + NH2O,anHHan2O(LB)

3169

work

(44)


Industrial & Engineering Chemistry Research (LBW ) =

Article

cath cath d(cair Cpair Tcath)

3.2. Steady-State Multiplicity. To investigate the steadystate behavior of the SOFC, the model is solved under nonisothermal conditions to consider the thermal behavior of the system as well. To this end, the heat generation and heat removal rates (QP and QR) are calculated under steady-state conditions and plotted versus the solid temperature, as is common in studying steady-state multiplicity in chemical reactors. On the basis of eq 47, the rates of heat generation and heat removal in the solid part of the cell are:

dt

cath nair,in ̇ cath Hair,in

cath cath cath cair Hair (BW ) + hcath(Ts − Tcath) − vair

(LB) − No2Hocath (LB) + NH2O,cathHHcath (LB) 2 2O

(45)

where han and hcath are the convective heat transfer coefficients for the fuel and the air streams, respectively. With the assumption of laminar boundary layer, the heat convection coefficients are calculated using the mean Nusselt number:21 hcathL = 0.664Prair 0.3ReairL 0.5 , kair h L Nu fuel = an = 0.664Prfuel 0.3RefuelL 0.5 k fuel

Nuair =

Q P = ΔHR LB − VI

(48)

Q R = LBhan(Tan − Ts) + LBhcath (Tcath − Ts)

(49)

The intersection points of the two plots represent the steadystates of the system. Bifurcation analysis is conducted under four operating modes: fixed ohmic load, potentiostatic, galvanostatic, and fixed cell power density (with bifurcation parameters being the external load resistance, cell voltage, cell current and cell power density, respectively). The inlet fuel and air flow rates are also considered as bifurcation parameters. To plot bifurcation diagrams, the nonlinear algebraic equations are solved at a given sequence of the bifurcation parameter values. For a given value of a bifurcation parameter, a numerical root finding method can find one steady state value depending on the chosen initial guess. One way to capture all the steady state values systematically is to find the roots by varying the initial guess. This method is not efficient, as it needs many trials and errors. On the other hand, for a given solid temperature, one load resistance always exists. Therefore, if we solve the algebraic equations for a given value of solid temperature rather than the external load resistance, one suitable initial guess is adequate to systematically find all steady states. The MATLAB function fsolve is used to solve the set of nonlinear algebraic equations. After the root corresponding to a given value of the solid temperature was obtained, for the next value of the solid temperature, the initial guess vector is chosen to be the root corresponding to the previous value of the solid temperature (zero order continuation). This approach reduces the convergence time and the number of numerical error messages.

(46)

An energy balance for the solid partaccounting for the flow of thermal energy in and out through the transport of reactants to the reaction sites and transport of products (water from the anode side and cathode side triple phase boundary) to the gas channel, respectively, as well as the convective heat transfer between the solid part and the fuel and air streams, and the electric power supplied to the external loadleads to ̂ dTs (dA + dE + dC)ρs Cps dt IV = ΔHR − + han(Tan − Ts) + hcath(Tcath − Ts) LB (47)

where ΔHR = NO2HOcath + NH2HHan2 − NH2O,cathHHcath 2 2O − NH2O,anHHan2O

Depending on the mode of operation of the cell, one should add another algebraic equation to the system of equations given above to arrive at a complete cell mathematical model. In the case of a constant load resistance operation, V = IRload, where the load resistance is known; for a potentiostatic operation, V = a known value; for a galvanostatic operation, I = a known value; and for a fixed-power operation, VI = a known value. The model parameter values are listed in Table 1.

4. PARAMETER ESTIMATION AND MODEL VALIDATION The polarization (I−V) curve predicted by the model is validated using the experimental data of Iwahara et al.8 obtained under isothermal operation at 1273 K. Figure 5 shows the experimental data and the model prediction. The unknown parameters, i0,an, i0,cath, γA, γC, EA, and EC are estimated from the experimental data. The exchange current densities, i0,an and i0,cath are assumed to be equal. To estimate the parameters, the MATLAB function fminsearch, which is based on the Nelder−Mead simplex algorithm, is used to calculate the parameter values that minimize the objective function:

3. SIMULATION METHODS 3.1. I−V Curve at Isothermal Conditions. The I−V curve of an actual fuel cell is usually obtained by linear sweep voltammetry using a potentiostat. This experiment is generally conducted in an isothermal furnace, thus isothermal conditions can be assumed.27 The prediction of the current−voltage (I−V) curve of the cell is an aim of this study. Here, to simulate the I− V curve, energy balance equations are omitted from the model and all temperatures in the remainder of the model are set to the constant temperature at which the experiment was conducted. Under steady-state conditions, the remaining model reduces to a set of algebraic equations. The MATLAB routine fsolve is used to find the root(s) of the algebraic equations. To calculate the cell current corresponding to a given cell outlet voltage (applied potential of the potentiostat), say Vg, we set V = Vg, we then solve the set of coupled algebraic equations. By repeating this step for different cell outlet voltages, we are able to generate the polarization curve.

i=1

J(Δ⃗) =

∑ (Imod ,i(Δ⃗) − Imea,i)2 + (Vmod ,i(Δ⃗) − Vmea,i)2 i=n

(50)

where Δ⃗ = (γ A , γ C , E A , E C , i0) 3170



Article

removal rates intersect at one point only, which exhibits a unique steady state. When the load resistance is 2.6 Ω, the lower and upper steady states are stable and the middle one is unstable. The bifurcation analysis is conducted under fixed ohmic load operation. In this case, the bifurcation parameter is the external load resistance. Figure 7 shows the bifurcation diagram of the

Figure 5. The I−V curve of the fuel cell at 1273 K (circles = experimental data,8 solid line = model prediction).

Imod,i and Imea,i are model predicted and experimental values of the cell current, respectively, and Vmod,i and Vmea,i are model predicted and experimental values of cell outlet voltage, respectively. Estimated values of the six model parameters are given in Table 1. Figure 7. Bifurcation analysis (fixed ohmic load mode of operation).

5. STEADY STATE BEHAVIOR Ignition, extinction, and hysteresis phenomena can occur in polyelectrolyte membrane fuel cells as well as SOFCs with an oxygen-conducting or proton-conducting electrolyte.28−31 A review of theoretical and experimental studies of steady-state multiplicity in polyelectrolyte membrane fuel cells and oxygenconducting SOFCs can be found in refs 2 and 32. Prior to this work, no study on steady-state multiplicity in co-ionic conducting SOFCs was reported. Detailed steady-state analysis of this type of cells is indispensible to determine whether steady-state multiplicity in the cell temperature can occur. Figure 6 depicts the heat generation and heat removal rates versus the solid temperature. The inlet cathode and anode channel gas temperatures are 890 K. As can be seen in Figure 6a, QP curve intersects QR line at three points when the applied load resistance is 2.6 Ω. As the load resistance is increased from 2.6 to 5 Ω, the cell current and consequently the heat generation rate decrease, and the heat generation and heat

solid temperature versus the load resistance. It can be seen that there is a region of three steady states between 2.4 and 2.8 Ω. Points A and C are the ignition and extinction points, respectively. At point A (ignition point), a small decrease in the load resistance leads to a jump in the cell temperature to the upper branch temperature, point B. The ignition phenomenon can have detrimental effects on the cell system; sudden large changes in the solid temperature can ruin the balance of plant components. As can be seen, the upper branch temperature is outside the range within which the model is valid. However, the results show the potential instability and ignition−extinction phenomenon in the system even if the ignited state is not approachable in reality.29 Mole fractions of different species in the anode and cathode gas channels versus the bifurcation parameter, the load resistance, are shown in Figure 8. It can be observed that concentration and temperature multiplicities coexist. The bifurcation curves for reactant species, O2 and H2, and water formed at the cathode side are in the form of an S-shaped curve, while that of the water formed at the anode side, is in the form of an inverse S-

Figure 8. Bifurcation analysis under fixed ohmic load: species mole fractions (yO2, yH2, yH2O,cath and yH2O,an) in the gas channels.

Figure 6. Steady-state heat removal and heat generation rates under two different external load resistances, (a) 2.6 Ω and (b) 5 Ω. 3171



Article

shaped curve. Figure 9 shows the cusp catastrophe steady-state manifold. One can observe that at higher inlet gas temperatures

Figure 11. Solid temperature versus load power density.

Figure 9. The bifurcation manifold of the cusp catastrophe: TS versus Rload and the inlet gas temperatures.

the multiple steady-state region disappears. The multiplicity region can be determined by projecting the turning points (ignition and extinction points) locus onto the parameter space. Figure 10 shows the cusp shape (shaded) region where three steady states are found.

Figure 12. Solid temperature versus the load resistance and equal inlet fuel and air flow rates.

increase in the flow rates decreases the anode and cathode channel gas temperatures and hence increases the heat dissipation by convection. Thus, the load resistance should decrease to compensate for the increase in the amount of heat dissipation. Under four operation modes of fixed ohmic load resistance, galvanostatic, potentiostatic, and fixed cell power density, the bifurcation diagrams of solid temperature versus the inlet fuel and air flow rates are shown in Figure 13. As can be seen, in the galvanostatic mode, there is a monotonic relationship between the solid temperature and the cell current, while under the three other modes, the presence of multiple steady states can be observed. The steady state multiplicity is also studied under potentiostatic mode (fixed cell voltage). This mode of operation requires voltage control and is used to characterize a cell. In this case, the cell voltage is the bifurcation parameter. Figure 14 shows the bifurcation diagram drawn by plotting steady state solid temperature as a function of cell voltage at equal inlet gas temperatures of 890 K. One can observe that more than one steady state exists in certain potential regimes. The middle one of the three steady states, shown by the dashed line, is unstable.

Figure 10. Multiplicity region in the fixed ohmic load mode.

The bifurcation diagram can be drawn with respect to the cell power density instead of the load resistance as shown in Figure 7. Figure 11 shows the bifurcation diagram of the solid temperature versus the cell power density. For given cell operating conditions, the cell power is maximum at a certain ohmic load. The range of the load power is 0−160 mW/cm2 for the single cell under study. The effect of fuel and air flow rates is shown in Figure 12. As can be seen, the range of the load for the single cell under study. The effect of fuel and air flow rates is shown in Figure 12. As can be seen, the range of the load resistance within which the multiplicity exists widens as the inlet fuel and air flow rate increases, and the inverse S-shaped curve moves toward lower values of the load resistance. The 3172



Article

Figure 13. Bifurcation diagrams: solid temperature versus equal inlet fuel flow rates under four different operation modes of fixed ohmic load, galvanostatic, potentiostatic, and fixed load power density.

behavior has also been reported for an oxygen ion-conducting SOFC.29,30

6. DYNAMIC BEHAVIOR The dynamic response of the cell to two step changes in the external load resistance is investigated. One step change is from 3 to 5 Ω (step size of 66%), and the other is from 2.6 to 3 Ω (step size of 15%). As shown in Figure 7, in the former case the cell does not undergo ignition or extinction, while in the latter case the cell undergoes extinction. The model presented in section 2 is a system of differential algebraic equations (DAEs), which is solved using the ode15s function of MATLAB to simulate the cell responses. The solid temperature, outlet voltage, total current, and reactants (oxygen and hydrogen) and product (water at the cathode and anode sides) outlet concentration responses of the cell to the step change from 3 to 5 Ω at time t = 500 s (prior to which the cell has been at steady state) are shown in Figure 16. As can be seen, the cell outlet voltage increases while the total current decreases as in a typical fuel cell.35 As the total current decreases, the rate of heat release by the electrochemical reactions decreases and thus the solid temperature decreases. A decrease in the solid temperature increases the electrolyte oxygen-ion resistivity, leading to a decrease in the oxygen-ion current. On the other hand, a decrease in the solid temperature decreases the electrolyte proton resistivity, resulting in an increase in the proton current. The more decrease in the oxygen-ion current and the less increase in the proton current lead to a decrease in the total current and a lower steady-state solid temperature. Two time scales can be seen clearly in the cell responses; an initial very quick change in each variable is followed by a slow change in the variable until the variable reaches its new steady state value. The quick change is governed by Ohm’s law, which has no dynamics (has a time constant of zero). The slow change is primarily governed by the cell thermal process, which has an apparent time constant of about 375 s. The solid temperature has a strong effect on the other variables; the other

Figure 14. Bifurcation diagram: solid temperature versus cell outlet voltage.

In the galvanostatic operation mode, the cell current is under control (is the bifurcation parameter). As shown in Figure 15, there is a unique steady state for all species concentrations, outlet gas temperatures, and cell temperature. This same

Figure 15. Solid temperature versus cell current. 3173



Article

Figure 16. Response of the SOFC to a step change from 3 to 5 Ω in the external load resistance.

variables do not settle (reach steady state) until the solid temperature does so; responses of variables other than the solid temperature and the solid temperature have the same long apparent time constant. As can be seen, it takes about 1500 s for the solid temperature to go from its initial steady-state value of 1164 K to its new steady-state value of 1108 K; the solid temperature response has an apparent time constant of about 375 s. The initial quick rise in the outlet voltage is due to the initial quick drop in the cell current, which quickly decreases the ohmic loss and increases the outlet voltage. The subsequent small slow drop in the outlet voltage is due the slow drop in the solid temperature, which slowly increases in the electrolyte oxygen-ion resistivity and increases the ohmic overpotential. The outlet reactant concentrations increase, as the reactants are consumed less at a lower total current. The concentration of water at the anode side (a reaction product) decreases, as it is expected (because as Figure 17 shows, the oxygen current decreases as the external load resistance increases). However, the concentration of water at the cathode side first slightly decreases but it then increases to a higher steady state value (it shows an inverse response). This inverse response can be explained by referring to Figure 17. The quick initial drop in the proton current is governed by Ohm’s law, which has no dynamics (has a time constant of zero). The proton current then increases slowly, because the proton conductivity of the electrolyte (shown in Figure 1) increases as the solid temperature drops. As can be seen, the slow increase in the proton current has an apparent time constant (about 375 s) equal to that of the solid temperature response, indicating that the slow increase is primarily governed by the cell thermal process. Thus, the slow increase in the proton current and the

Figure 17. Responses of the oxygen-ion and proton currents to the step change from 3 to 5 Ω in the external load resistance.

lower temperature of the cathode gas stream together are responsible for the inverse response in the concentration of water at the cathode side. The solid temperature, outlet voltage, total current, and reactants (oxygen and hydrogen) and product (water at the cathode and anode sides) outlet concentration responses of the cell to the step change from 2.6 to 3.0 Ω in the external load resistance at time t = 500 s (prior to which the cell has been at steady state) are shown in Figure 18. As can be seen in Figure 7, when the external load is 2.6 Ω, the cell operates on the 3174



Article

Figure 18. Response of the SOFC to a step change from 2.6 to 3 Ω in the external load resistance.

upper branch. The step change from 2.6 to 3 Ω in the external load resistance forces the cell to undergo extinction in the solid temperature (which drops from the upper branch [1795 K] to the lower branch [1166 K]). Unlike the cell response to the step change from 3 to 5 Ω in the external load resistance, here the total current and outlet voltage both decrease. This behavior is in agreement with the dynamic behavior of an oxygen-ion conducting SOFC subjected to a step change of the same kind,30 but it is not a typical behavior of a SOFC. The outlet voltage exhibits inverse response. The initial quick small rise in the outlet voltage is due to the initial quick small drop in the cell current, which quickly decreases the ohmic loss and increases the outlet voltage. The subsequent large slow drop in the outlet voltage is due the large slow drop in the solid temperature, which slowly increases in the electrolyte oxygenion resistivity and increases the ohmic overpotential. Indeed, the inverse response of the outlet voltage is a result of a completion between Ohm’s law (electronics) and heat transfer. The decrease in the steady-state solid temperature is caused by the lower rate of heat release by the electrochemical reactions at the lower total current. As the total current decreases, the rate of heat release by the electrochemical reactions decreases and thus the solid temperature decreases. A decrease in the solid temperature increases the electrolyte oxygen-ion resistivity, leading to a decrease in the oxygen-ion current. On the other hand, a decrease in the solid temperature decreases the electrolyte proton resistivity, resulting in an increase in the proton current. The more decrease in the oxygen-ion current and the less increase in the proton current lead to a decrease in the total current and a lower steady-state solid temperature. The solid temperature drops from its initial steady-state value of 1795 K to its new steady-state value of 1166 K in about 7500 s. Again, two time scales can be seen clearly in these responses; an initial very quick change in each variable is followed by a

slow change in the variable until the variable reaches its new steady state value. As before, the quick change is governed by Ohm’s law, which has no dynamics. The slow change is primarily governed by the cell thermal process, which has an apparent time constant of about 1875 s. This apparent time constant is significantly higher than the 375 s time constant that was observed for the solid temperature response to the 3 to 5 Ω step change in the external load resistance. As the solid temperature has a strong effect on the other variables, we see again that the other variables do not settle (reach steady state) until the solid temperature does so. The significantly different apparent time constants (375 vs 1875 s) and steady-state gains (−28 vs −1572.5 K/Ω) observed in the solid temperature responses to the two step changes in the external load resistance are indicative of the very high nonlinearity of the SOFC.

7. CONCLUDING REMARKS A mathematical model was developed to investigate the presence of multiplicity in a coionic conducting SOFC. Simulation results showed that a multiple steady state region exists in this type of SOFCs. Ignition and extinction phenomena in the solid temperature and outlet gas concentrations were observed. This result is in agreement with behavior reported for oxygen ion-conducting SOFCs and proton-conducting SOFCs, attributing the existence of multiplicity to the dependence of the electrolyte oxygen ion conductivity to the electrolyte temperature. Under fixed ohmic load and potentiostatic operations, it was found that concentration and temperature multiplicities coexist. Ignition and extinction behavior in the cell temperature is reversed for the reactants concentrations. While the reactants concentrations and water product at the cathode side ignite, the cell temperature and the water product concentration at the anode side become extinguished and vice versa. Qualitatively, coionic 3175



Article

ρs = density of the solid (kg m−3) σji = characteristic length (Å), j, i: H2, H2O, O2, N2 τ = tortuosity ΩD = diffusion collision integral

conducting SOFCs show steady-state behavior similar to that of oxygen ion-conducting SOFCs and proton-conducting SOFCs When the cell power density is the bifurcation parameter, the presence of steady state multiplicity was found. This multiplicity can pose a challenge when power control is needed. Also, it was shown that as the inlet fuel and air flow rates increase, the region of steady-state multiplicity expands and shifts toward lower values of load resistance. The study of the dynamic behavior of the cell also showed that the cell can exhibit inverse response. Significantly different responses of the cell in terms of steady-state gain and apparent time constant in the two different operating regions are indicative of the high nonlinearity of the fuel cell.

■

Subscripts

an = anode cath = cathode E = electrolyte S = solid Superscripts

■

AUTHOR INFORMATION

TPB = triple phase boundary

REFERENCES

(1) Singhal, S. C.; Kendall, K. High Temperature Solid Oxide Fuel Cells: Fundamentals, Design, And Applications; Elsevier Science: Oxford, UK, 2003. (2) Bavarian, M.; Soroush, M.; Kevrekidis, I. G.; Benziger, J. B. Mathematical modeling, steady-state and dynamic behavior, and control of fuel cells: A review. Ind. Eng. Chem. Res. 2010, 49 (17), 7922−7950. (3) Demin, A.; Tsiakaras, P.; Gorbova, E.; Hramova, S. A SOFC based on a co-ionic electrolyte. J. Power Sources 2004, 131 (1−2), 231−236. (4) Gorbova, E.; Maragou, V.; Medvedev, D.; Demin, A.; Tsiakaras, P. Investigation of the protonic conduction in Sm doped BaCeO3. J. Power Sources 2008, 181 (2), 207−213. (5) Jiang, K.; He, Z.; Meng, J.; Ren, Y.; Su, Q. Low temperature preparation and fuel cell properties of rare earth doped barium cerate solid electrolytes. Sci. China Ser. B 1999, 42 (3), 298−304. (6) Glöckner, R.; Islam, M.; Norby, T. Protons and other defects in BaCeO3: A computational study. Solid State Ion 1999, 122 (1−4), 145−156. (7) Iwahara, H.; Hibino, T.; Okada, T. In proton and oxide-ion mixed conduction in perovskite-type oxide ceramics based on cerates. Proc. 2nd Int. Symp. Ionic Mixed Conduct. Ceram. 1994, 1−8. (8) Iwahara, H.; Yajima, T.; Hibino, T.; Ushida, H. Performance of solid oxide fuel cell using proton and oxide ion mixed conductors based on BaCe1‑xSmxO3‑α. J. Electrochem. Soc. 1993, 140, 1687. (9) Riess, I. Current-voltage relation and charge distribution in mixed ionic electronic solid conductors. J. Phys. Chem. Solids 1986, 47 (2), 129−138. (10) Ross, P., Jr; Benjamin, T. Thermal efficiency of solid electrolyte fuel cells with mixed conduction. J. Power Sources 1977, 1 (3), 311− 321. (11) Yuan, S.; Pal, U. Analytic solution for charge transport and chemical potential variation in single layer and multilayer devices of different mixed conducting oxides. J. Electrochem. Soc. 1996, 143, 3214. (12) Riess, I. Mixed ionic-electronic conductorsmaterial properties and applications. Solid State Ion. 2003, 157 (1−4), 1−17. (13) Huang, J.; Yuan, J.; Mao, Z.; Sunden, B. Analysis and modeling of novel low-temperature SOFC with a co-ionic conducting ceriabased composite electrolyte. J. Fuel Cell Sci. Technol. 2010, 7 (1), 011012−7. (14) Bhattacharyya, D.; Rengaswamy, R. A review of solid oxide fuel cell (SOFC) dynamic models. Ind. Eng. Chem. Res. 2009, 48 (13), 6068−6086. (15) Hajimolana, S. A.; Hussain, M. A.; Daud, W. M. A. W.; Soroush, M.; Shamiri, A. Mathematical modeling of solid oxide fuel cells: A review. Renew. Sustain. Energy Rev. 2011, 15 (4), 1893−1917. (16) Incropera, F. P. Fundamentals of Heat and Mass Transfer, 6th ed.; John Wiley: Hoboken, NJ, 2007. (17) Bove, R.; Ubertini, S. Modeling Solid Oxide Fuel Cells: Methods, Procedures and Techniques; Springer Verlag: Great Britain, 2008. (18) Costamagna, P.; Honegger, K. Modeling of solid oxide heat exchanger integrated stacks and simulation at high fuel utilization. J. Electrochem. Soc. 1998, 145, 3995.

Corresponding Author

§E-mail: [email protected]. Tel.: (215) 895-1710. Fax: (215) 895-5837. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS M.B. and M.S. acknowledge the National Science Foundation (NSF) through Grant CBET-0932882. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.

■

NOMENCLATURE B = height of cell (m) Cp = specific heat capacity of the solid (J mol−1 K−1) Ĉ ps = specific heat capacity of the solid (J kg−1 K−1) dA = thickness of anode (m) dC = thickness of cathode (m) dE = thickness of electrolyte (m) F = Faraday constant (A s mol−1) h = heat transfer coefficient (W m−2 K−1) H = enthalpy (J mol−1) I = cell current (A) i0,an = exchange current density at anode (A m−2) i0,cath = exchange current density at cathode (A m−2) L = length of cell (m) ṅ = molar flow rate (mol s−1) Rload = external load resistance (Ω) T = temperature (K) V = outlet cell voltage (V) v = velocity (m s−1) W = width of gas channels (m) yj = mole fraction of species j

Greek Letters

α = symmetry factor δ = pore mean radius (m) Δ = thickness of the diffusion layer (m) εan = porosity of anode εcath = porosity of cathode εji = the characteristic Lennard-Jones energy (J), j, i: H2, H2O, O2, N2 ηconc = concentration overpotential (V) ηan = anode activation overpotential (V) ηcath = cathode activation overpotential (V) θ = charge transfer coefficient ρ = resistivity (Ω m) 3176



Article

(19) Patcharavorachot, Y.; Brandon, N. P.; Paengjuntuek, W.; Assabumrungrat, S.; Arpornwichanop, A. Analysis of planar solid oxide fuel cells based on proton-conducting electrolyte. Solid State Ion 2010, 181, 1568−1576. (20) Bard, A. J.; Faulkner, L. R. Electrochemical Methods: Fundamentals and Applications, 2nd ed.; John Wiley: New York, 2001. (21) Welty, J.; Wicks, C.; Rorrer, G.; Wilson, R. Fundamentals of Momentum, Heat, And Mass Transfer; Wiley: New York, 2009. (22) Ni, M. Modeling of a planar solid oxide fuel cell based on proton conducting electrolyte. Int. J. Energy Res. 2010, 34, 1027−1041. (23) Chan, S.; Khor, K.; Xia, Z. A complete polarization model of a solid oxide fuel cell and its sensitivity to the change of cell component thickness. J. Power Sources 2001, 93 (1−2), 130−140. (24) Wilke, C. R. Diffusional properties of multicomponent gases. Chem. Eng. Prog 1950, 46 (2), 96−104. (25) Wang, C. Modeling and Diagnostics of Polymer Electrolyte Fuel Cells; Springer Verlag: New York, 2010; Vol. 49. (26) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases and LiquidsMcGraw-Hill: New York, 1987. (27) Suzuki, M.; Shikazono, N.; Fukagata, K.; Kasagi, N. Numerical analysis of coupled transport and reaction phenomena in an anodesupported flat-tube solid oxide fuel cell. J. Power Sources 2008, 180 (1), 29−40. (28) Debenedetti, P.; Vayenas, C. Steady-state analysis of high temperature fuel cells. Chem. Eng. Sci. 1983, 38 (11), 1817−1829. (29) Mangold, M.; Krasnyk, M.; Sundmacher, K. Theoretical investigation of steady state multiplicities in solid oxide fuel cells. J. Appl. Electrochem. 2006, 36 (3), 265−275. (30) Bavarian, M.; Soroush, M. Steady-state multiplicity in a solid oxide fuel cell: Practical considerations. Chem. Eng. Sci. 2012, 67 (1), 2−14. (31) Bavarian, M.; Soroush, M. Mathematical modeling and steadystate analysis of a proton-conducting solid oxide fuel cell. J. Process Control 2012, 22 (8), 1521−1530. (32) Hanke-Rauschenbach, R.; Mangold, M.; Sundmacher, K. Nonlinear dynamics of fuel cells: a review. Rev. Chem. Eng. 2011, 27 (1−2), 23−52. (33) Iwahara, H.; Uchida, H.; Morimoto, K. High temperature solid electrolyte fuel cells using perovskite type oxide based on BaCeO. J. Electrochem. Soc. 1990, 137, 462. (34) Hirschfelder, J. O.; Curtiss, C. F.; Bird, R. B., Molecular Theory of Gases and Liquids; Wiley: New York, 1954. (35) Hajimolana, S. A.; Soroush, M. Dynamics and control of a tubular solid-oxide fuel cell. Ind. Eng. Chem. Res. 2009, 48 (13), 6112− 6125.

3177


Modeling and Bifurcation Analysis of a Coionic Conducting Solid

Recommend Documents