A model for SOFC anode performance

Electrochimica Acta 54 (2009) 6686–6695

Contents lists available at ScienceDirect

Electrochimica Acta journal homepage: www.elsevier.com/locate/electacta

A model for SOFC anode performance A.A. Kulikovsky a,b,∗,1 a b

Institute for Energy Research – Fuel Cells (IEF-3), Research Centre “Jülich”, D-52425 Jülich, Germany Moscow State University, Research Computing Center, 119991 Moscow, Russia

a r t i c l e

i n f o

Article history: Received 12 March 2009 Received in revised form 22 June 2009 Accepted 22 June 2009 Available online 28 June 2009 Keywords: SOFC anode Activation resistivity Modeling

a b s t r a c t A model for anode performance of a planar anode-supported SOFC is developed. The model includes Butler–Volmer relation for the hydrogen oxidation, Ohm’s law for ionic current and equation of hydrogen mass balance in the anode channel. We show that the regime of anode operation depends on the relation between the cell current density j and the critical current density jcrit . Analytical solutions to the system of governing equations for the case of “low” (j < jcrit ) and “high” (j  jcrit ) currents are derived. In the “low-current” regime the anode polarization voltage is proportional to cell current, which justifies the notion of anodic activation resistivity Ra . Full hydrogen utilization increases the value of Ra by a factor of 2. In the “high-current” regime polarization voltage depends on cell current logarithmically, with the effective Tafel slope being twice the kinetic value (doubling of Tafel slope). In this regime 100% hydrogen utilization leads to a constant  230-mV shift of polarization curve as a whole. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction Polarization curve of a solid oxide fuel cell (SOFC) is close to linear [1]. This justifies introduction of cell area-specific resistivity Rcell , which is the slope of the line Vcell (J) Vcell = Voc − Rcell J

(1)

Here Vcell is the cell voltage, Voc is the open circuit voltage and J is the mean current density. In spite of its simplicity (or rather because of it), Eq. (1) raises several questions. To a good approximation kinetics of electrochemical reactions on both sides of the cell follow Butler–Volmer law, which establishes exponential dependence of cell current on overpotential. Polarization curve of a low-temperature hydrogen cell exhibits logarithmic shape in the region of small currents (activation region) and a rapid drop to zero close to the limiting current [2]. Why the polarization curve of SOFC is linear? Does it mean that the resistive losses in SOFC dominate and the contribution of activation polarizations to the overall voltage loss is small? A model of SOFC developed by Imperial College group [3] includes Butler–Volmer kinetics and transport processes on both sides of the cell. The results [3] show that for cell temperature 800 + 273 K hydrogen diffusion to the catalyst sites has minor effect

∗ Corresponding author at: Institute for Energy Research – Fuel Cells (IEF-3), Research Centre “Jülich”, D-52425 Jülich, Germany. Tel.: +49 02461 61 5396; fax: +49 02461 61 6695. E-mail address: [email protected]. 1 ISE member. 0013-4686/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.electacta.2009.06.054

on cell performance. However, the anode activation polarization calculated in [3] is quite significant, it is comparable to the resistive voltage loss in a cell. Thus, understanding the nature of anode activation loss is of large interest for SOFC technology. Over the past decade a number of CFD models for SOFC performance have been developed. Some of these models rely on empirical notion of ASR, not detailing the kinetics of electrochemical reactions [4,5]. The other utilize Butler–Volmer equation for calculation of activation losses [6,7,3,8–11]. In recent years a novel class of micro-structural electrode models has emerged. These models are based on conservation equations in either random [12] or stochastically reconstructed electrode structures [13,14]. The conservation equations in such a structure are solved using either direct numerical simulation [13] or Lattice–Boltzmann method [14,15]. These models reveal the structure–function relations for real electrodes (e.g., Bruggemann exponent), which in conventional hydrodynamic models are taken as empirical values. However, all these models are numerical and they do not give the irrefragable answer to the questions above. Much less has been done in analytical description of SOFC polarization curve. A through-plane model, which includes equations for diffusive transport of reactants and current conservation equations coupled by Butler–Volmer terms is developed in [16]. An approximate analytical (asymptotic) solution to model equations is obtained. However, the solution [16] is rather cumbersome and difficult for interpretation. Analytical model [17] describes performance of SOFC electrode in the limit of small overpotential . Smallness of  allows to linearize Butler–Volmer equation, which greatly simplifies the

A.A. Kulikovsky / Electrochimica Acta 54 (2009) 6686–6695

problem. This approximation, however, does not give a picture of electrode function at large , neither it describes the polarization curve in the transition region from small to large anodic overpotentials. Besides, the model [17] is one-dimensional and it ignores the effect of hydrogen depletion along the channel on anode performance. 1D analytical model of SOFC sandwich [18] leads to implicit expressions for the electrode voltage–current relation in the limits of overpotential. This model is designed to assist the multidimensional CFD simulations of SOFC; in the CFD environment solution of implicit equations poses no problems. However, the implicit relations [18] are not suitable for analysis of the effect of hydrogen utilization on cell performance. In this work we report an analytical model, which clarifies the situation with the anode activation polarization voltage in the whole range of current densities. The model is close in spirit to [17], though here we report a solution to non-linear governing equations, which follow from the non-linearized Butler–Volmer law. We show that there exists a critical current density jcrit , which separates “low-” and “high-current” regimes of SOFC anode operation. If cell current density j is lower than jcrit , the dependence of  on j is linear, which justifies the notion of polarization resistivity Ra ; a simple relation for Ra is derived. When j is much larger than jcrit ,  depends on cell current logarithmically, with the effective Tafel slope being twice the kinetic value. In this regime there is no simple proportionality between current and activation polarization. Then we rationalize the effect of hydrogen utilization u on the anodic polarization. In the “low-current” limit the anode resistivity Ra increases with u; this growth is dramatic for u  80%, leading to doubling of Ra at 100% utilization. In the “high-current” mode the effect of hydrogen utilization is also quite significant: 100% utilization shifts the cell polarization curve down by the value 2b ln 2  235 mV, where b  170 mV is the characteristic Tafel slope. 2. Through-plane model

Fig. 1. The shape of ionic current across the anode.

the transfer coefficient,  is the polarization voltage of the anode side and i is the anode ionic conductivity. Eq. (2) expresses the decay of ionic current as one moves from the anode/electrolyte interface to the anode depth (Fig. 1). The right side of this equation is the Butler–Volmer rate Q of the electrochemical conversion. Eq. (3) is Ohm’s law relating ionic current with the gradient of overpotential. Based on reaction stoichiometry, the rate of hydrogen oxidation is taken to be proportional to the hydrogen concentration cH2 . Impedance measurements show that Q may also depend on water concentration cw (see [21] and the literature cited therein). Model−0.5 ing study [21] indicates that Q ∼cw ; however, this dependence is weak and to a first approximation it can be ignored. Following [22], the transfer coefficient ˛ for the reaction with the rate-determining step involving single-electron transfer is 0.5. Measurements [19] gave ˛  0.7; however, to simplify calculations we will adopt the value ˛ = 1 − ˛ = 0.5. Note, however, that this choice of ˛ is important only in the region of small currents, when both exponents in Butler–Volmer equation contribute to the activation polarization. At larger currents the second exponent in Eq. (2) is small and the results below are valid for arbitrary values of ˛. We introduce dimensionless variables

2.1. Basic equation for the shape of ionic current in the anode x˜ = Measurements show that at large anodic overpotentials hydrogen oxidation follows Tafel law [19]. Hydrogen diffusion coefficient in the anode material is large and we will ignore hydrogen transport losses. Detailed model [20] shows that the variation of hydrogen concentration across the anode thickness is marginal. The estimate of voltage loss due to hydrogen transport is given in Appendix B. However, ionic conductivity of cermet is small and voltage loss due to ion transport through the anode should be taken into account. Local anode performance is thus governed by two equations: ∂j(x) = i∗ ∂x



cH2 cref

∂(x) ∂x



exp

 ˛F  RT

 (1 − ˛)F 

− exp −

RT

(2)

b=

 ˜=

 , b

˜j = jla i b

(4)

RT ˛F

(5)

is the Tafel slope. With these variables Eqs. (2) and (3) take the form ε2

∂˜j = sinh() ˜ ∂˜x

(6)

˜ ˜j = ∂ ∂˜x

(7)

where ε=

Here the axis x is directed from the interconnect to the electrolyte (Fig. 1), j(x) is the local ionic current density,2 i∗ is the volumetric exchange current density (a number of charges produced in unit volume per second, A m−3 ), cH2 is the molar concentration of hydrogen in the anode, cref is the reference hydrogen concentration, ˛ is

and

2 In this Section we do not consider the dependence of j on the distance along the anode channel; this will be done in Section 3.

x , la

where la is the anode thickness and

(3)

j(x) = i

6687

l∗ =

l∗ la

(8)

 i b 2i∗



cref cH2

 (9)

is the characteristic reaction penetration depth into the anode i.e., the thickness of the reaction zone near the electrolyte. In a cell of the anode-supported type, the anode thickness la largely exceeds l∗ and thus ε is a small parameter.

6688

A.A. Kulikovsky / Electrochimica Acta 54 (2009) 6686–6695

The system of Eqs. (6) and (7) can be transformed to a single equation for ionic current. Differentiating Eq. (6) with respect to x˜ we get ε2

∂2˜j ∂ ˜ = cosh() ˜ ∂˜x ∂˜x2

Using Eq. (7) and the identity cosh2 () ˜ − sinh2 () ˜ = 1 we come to



ε2

∂2˜j = ˜j ∂˜x2

1 + sinh2 () ˜

Excluding now sinh2 () ˜ with Eq. (6) we obtain



∂2˜j ε2 2 = ˜j ∂˜x

 1 + ε4

∂˜j ∂˜x

2 (10)

Eq. (10) determines the shape of ionic current across the anode active layer. Solution to this equation subject to the boundary conditions ˜j(0) = 0,

˜j(1) = ˜j1

(11)

where ˜j1 is the ionic current density at the electrolyte interface i.e., the cell current density. We see that the problem is governed by two parameters, ε and ˜j1 . The general analytical solution to Eq. (10) hardly exists. However, for the anode-supported design ε  1. With ε  1 we have just one free parameter: ˜j1 . In the limiting cases of small and large ˜j1 Eq. (10) can be solved. 2.2. The shapes of local current and overpotential across the anode 2.2.1. “Low-current” regime Consider first the case of small cell current ˜j1  1. In that case the term ε4

 2 ∂˜j ∂˜x

in Eq. (10) can be neglected (a rigorous condition for

this approximation will be obtained below). The resulting equation ε2

∂2˜j = ˜j, ∂˜x2

˜j(0) = 0,

˜j(1) = ˜j1

(12)

can easily be integrated to yield ˜ ˜j(˜x) = j1 (exp(˜x/ε) − exp(−˜x/ε)) exp(1/ε) − exp(−1/ε)

(13)

Since ε is small, 1/ε is large and the exponents with the negative power can be dropped off.3 We finally obtain ˜j(˜x) = ˜j1 exp

 −(1 − x˜ ) 

(14)

ε

Now the limits of validity of Eq. (14) can be established. From (14) we have ∂˜j/∂˜x = ˜j/ε and the right side of (10) takes the form



˜j

 1 + ε4

∂˜j ∂˜x

2



= ˜j

 1 + ε2˜j2  ˜j

1+

ε2˜j2 2

 (15)

We see, that the term with ε can be neglected if ε2˜j2 /2

 1. Peak of local current density is achieved at the electrolyte interface. Thus, the latter inequality is equivalent to ε2˜j12 2

1

Fig. 2. Dimensionless ionic current density vs. distance across the anode. Interconnect is at x/la = 0 and electrolyte is at x/la = 1. Points and dashed curve: the exact numerical solution to Eq. (10). Solid lines: the approximate analytical solution (14). Dotted line in the bottom frame: the large-current solution (31).

(16)

Eq. (16) determines the limits of validity of Eq. (14).

The exact numerical solution to the full equation Eq. (10) and approximate expression (14) are compared in Fig. 2 for the three values of ˜j1 : 1, 10 and 100. Note that in all three cases ε = 0.1 i.e., the reaction penetration depth is the same: l∗ = la /10. For ε = 0.1 and ˜j1 = 10 we have ε˜j1 = 1 and strictly speaking Eq. (16) is not fulfilled. However, due to the factor 1/2 in Eq. (16), Eq. (14) still approximates the exact result with a good accuracy (Fig. 2). This allows us to mitigate Eq. (16) replacing it with ε˜j1 < 1

3 Omitting exp(−˜x/ε) we distort the shape of ˜j(˜x) near x˜ = 0 (interconnect). However, ionic current at x˜ = 0 is vanishingly small and this distortion is also negligible.

(17)

It is advisable to rewrite Eq. (17) in dimension form. Using (4), (8) and (9) we obtain

A.A. Kulikovsky / Electrochimica Acta 54 (2009) 6686–6695

 j1 < jcrit

cH2

(18)

cref

jcrit =

2i∗ i b

(19)

is the critical current density, below which Eq. (14) is valid (here for the estimate we put cH2 = cref ). High exchange current density and/or ionic conductivity and/or Tafel slope increase the

value of jcrit . Note that the right side of Eq. (18) is proportional to cH2 , that is low hydrogen concentration may distort the exponential shape of ionic current density in the anode. Using ∂˜j/∂˜x = ˜j/ε in Eq. (6) we get equation for overpotential exp() ˜ − exp(−) ˜ = 2ε˜j

(20)

Solving this for  ˜ we find



 ˜ = ln ε˜j +



1 + ε2˜j2



(21)

ε3˜j3 + ... 6

(22)

˜j2 ∂˜j ˇ2 = − 2 2 ∂˜x



˜j = ˇ tan

(29)

 2

−y

−1

˜j1 2 + ˜j1

and thus

This allows us to neglect 1 under the root sign in Eq. (10). This equation takes the form



Using this in Eq. (29) and solving the resulting equation for ˇ we get

 ˜ = ε˜j

2.2.2. “High-current” regime Consider now the case of large cell current, ε˜j1  1. Rapid variation of ionic current and overpotential occurs in the active sublayer near the electrolyte interface 1 − ε ≤ x˜ ≤ 1. Below we will see that in the “high-current” regime ∂˜j/∂˜x  ˜j2 /2. Thus, if ε˜j1  1, the

2 inequality ε4 ∂˜j/∂˜x  ε4˜j4 /4  1 is guaranteed near electrolyte.4

ˇ˜x 2

The local current density near the electrolyte is, therefore, described by tan-function.5 For large ˜j the argument of tan-function in (29) is close to /2 and the following approximation holds:

ˇ

Thus, the shape of overpotential across the anode is also exponential. The solid curves in Fig. 2 hence depict also (˜ ˜ x), provided that the y-axis in this figure is re-scaled by a factor ε.

(28)

Solution to this equation subject to the boundary condition ˜j(0) = 0 is

If (17) is fulfilled, we can safely neglect the second term on the right side of Eq. (22). This finally gives (23)

(27)

and using this in (25) we get

tan(y) 

(the other root is negative). Expanding (21) over ε˜j we obtain  ˜ = ε˜j −

Below we will show that the right side of (26) is positive. Introducing the notation exp  ˜ 1 ˜2 − j1 = ˇ2 ε2

where



6689



˜j =

˜j1 2 + ˜j1

(30)



 tan

˜j1 x˜ 2(2 + ˜j1 )

 (31)

This solution for ˜j1 = 100 is shown in the bottom frame in Fig. 2. The x-dependence of local polarization voltage is obtained from Eq. (6). Using (29) in (6) and solving for  ˜ we come to



ε2 ˇ2 2

(˜ ˜ x) = ln



 2

1 + tan

x˜ ˇ 2

  +

ε4 ˇ4 1+ 4



 1 + tan

2

x˜ ˇ 2

2



(32)

where ˇ is given by (30). The plots (˜ ˜ x) for the three values of ˜j1 are depicted in Fig. 3. Note the rapid growth of (1) ˜ with the increase in ˜j1 .

∂˜j ∂2˜j = ˜j 2 ∂˜x ∂˜x

2.3. Anode polarization curve

or

2.3.1. “Low-current” regime Setting in (23) x˜ = 1 we find a “low-current” polarization curve of the anode

∂2˜j 1 ∂˜j2 − =0 2 2 ∂˜x ∂˜x

(24)

 ˜ 1 = ε˜j1

Integrating once we get



˜j2 ∂˜j  ∂˜j − = 2 ∂˜x ∂˜x 

− x˜ =1

˜j2 1

(25)

2

Consider the right side of Eq. (25). In view of (6) we can transform it to



∂˜j  ∂˜x 

˜j2 1

x˜ =1

˜j2 sinh  ˜1 1 1 = = − − 2 2 2 ε2

 exp ˜

1

ε2

− ˜j12



We see that the activation voltage loss is proportional to current density. This justifies introduction of anode activation resistivity. In dimension variables Eq. (33) is 1 = Ra j1

We are interested in asymptotic behavior of  ˜ 1 vs. ˜j1 and we, therefore, can tolerate rather high absolute error in the shapes ˜j(˜x) and (˜ ˜ x) far from the electrolyte interface.

(34)

where the anode activation resistivity is given by



(26)

The last equality holds since large current means that  ˜ 1 is also large and the exponent with the negative power in the expression for sinh( ˜ 1 ) can be neglected.

4

(33)

Ra =

5

b 2i∗ i



cref



cH2

(35)

From Eq. (28) we have

˜j2 ∂˜j 1 = (ˇ2 + ˜j2 )  2 2 ∂˜x since ˇ2  2  ˜j2 (see below). This justifies the estimates in the beginning of this section.

6690

A.A. Kulikovsky / Electrochimica Acta 54 (2009) 6686–6695

Fig. 3. Dimensionless anodic overpotential  ˜ (32) as a function of the distance across the anode for indicated values of cell current density ˜j1 . The electrolyte is at x/la = 1, ε = 0.1.

2.3.2. “High-current” regime Polarization curve of the anode side in the “high-current” regime results from the following arguments. Setting x˜ = 1 in Eq. (29) we obtain

 

3. Hydrogen utilization and the anode voltage loss

ˇ 2

˜j1 = ˇ tan

3.1. Hydrogen concentration in the channel and cell voltage

Since ˜j1 is large, the argument of tan-function must be on the order of /2: ˇ   2 2 Using here ˇ (27) and solving for  ˜ 1 we find6  ˜ 1  ln(ε2 (2 + ˜j12 ))

(36)

Since ε˜j1  1, we have ˜j12  1/ε2 i.e., ˜j12  100 and thus 2 in Eq. (36) can be neglected. This yields



 ˜ 1  2 ln ε˜j1

(37)

Eq. (37) in dimension form is 1  2b ln

Fig. 4. The general polarization curve of SOFC anode. In the low-current region ˜ 1 = ε˜j1 . In the high-current region (ε˜j1 ≥ 10) the curve is (ε˜j1 < 1) it is linear,  non-linear,  ˜ 1 = 2 ln(ε˜j1 ). In the intermediate region (1 < ε˜j1 < 10) the “patching” function (39) connects the two analytical solutions. Note that the point ε˜j1 = 1 corresponds to the critical current density. Points: the exact numerical polarization curve, which follows from numerical solution to a problem (10); voltage loss has been determined from (6).

j R  1 a

All relations above include hydrogen concentration cH2 . This concentration is constant along the anode channel if the cell is run at low hydrogen utilization. However, real cells are usually run at high fuel utilization and cH2 decreases towards the anode channel outlet. To rationalize the effect of H2 utilization on anode resistivity we consider now SOFC with the linear anode channel (Fig. 5). The assumptions we will use are as follows: (1) Cell temperature is uniform. This assumption is justified when air stoichiometry is large, so that the cooling effect of air is also high. (2) Flow velocity in the anode channel is constant (plug flow).

(38)

b

where Ra is the anode resistivity (35). We see that in the high-current regime the anode polarization curve is logarithmic; moreover, factor 2 on the right side indicates doubling of apparent Tafel slope. In this regime the notion of anode polarization resistivity is no longer valid: there is no simple proportionality between 1 and j1 . The general dimensionless polarization curve of SOFC anode is depicted in Fig. 4. At ε˜j1 < 1 the curve is linear; above ε˜j1 = 10 this curve is given by the function 2 ln(ε˜j1 ). In the intermediate region the “patching” function  ˜ 1 = (t + (1 −  2.2 )2 ln t),

1 ≤ t ≤ 10

(39)

connects the two analytical solutions. Here  = (t − 10)/9 and t ≡ ε˜j1 . Points in Fig. 4 show the “exact” numerical curve, which results from solution of Eqs. (10) and (6).7

6 With (36) we have ˇ2 = (exp  ˜ 1 )/ε2 − ˜j12  2 > 0, which justifies introduction of ˇ. 7 At ε˜j > 4 solution to (10) was obtained introducing the stretched coordinate X = x˜ /ε2 and using the parametric continuation method.

Fig. 5. Anode channel and schematic of the hydrogen concentration profile.

A.A. Kulikovsky / Electrochimica Acta 54 (2009) 6686–6695

Under these assumptions H2 consumption in the channel is described by a simple mass balance equation:

v0a

∂cH2

j =− 2Fhc ∂z

(40)

where z is the coordinate along the channel (Fig. 5), v0a is the anode flow velocity, hc is the height of the anode channel. Hereinafter the subscript “1” is omitted (j ≡ j1 ); the superscript “0” indicates the values at the channel inlet. It is convenient to introduce dimensionless variables: z˜ =

z , L

ˆj =

j jcrit

,

c˜ =

cH2

(41)

0 cH

2

centration and jcrit is given by (19).8 With these variables Eq. (40) takes the form ˆJ

∂˜c = −ˆj ∂˜z

(43)

In the low-current mode the anode polarization curve is given by Eq. (34). With (41), Eq. (34) takes the form



 ˜ = ˆj

c˜ref

(45)

c˜

The system of Eqs. (43)–(45) determines the shapes of hydrogen molar fraction and local current along the channel. Using (45) in (44) and solving the resulting equation for ˆj we find E˜

ˆj =

˜+ R



Substituting this into (43) we come to equation for hydrogen concentration c˜ ˆJ

∂˜c E˜ , =−

∂˜z ˜ R + c˜ref /˜c

c˜ (0) = 1

(47)

Solution to this equation is 1 c˜ (˜z ) = ˜2 R

0 2Fhc v0a cH

(46)

c˜ref /˜c



where, by the definition, hydrogen stoichiometry =

3.2. Low current: the z-shapes

2

0 is the inlet hydrogen molar conwhere L is the channel length, cH

6691



−



 ˜ c˜ref + R

c˜ref +



˜ +2 R

c˜ref

˜z E˜ − ˆJ

2 (48)

2

LJ

Cell voltage V˜ cell can be written as ˜ ˆj V˜ cell = V˜ oc −  ˜ −R ˜ accumulates all resistivities, Here V˜ oc is open-circuit voltage and R ˜ is scaled as except the anodic one. R Rj ˜ = crit R b ˜ includes the cathodic resistivity, which is assumed to be Note that R independent of local cell current and thus independent of z˜ . Typically, due to small cathode thickness the respective concentration polarization for currents below 1 A cm−2 is negligible [23]. The situation with the activation polarization on the cathode side is more subtle. The assumption of constant activation resistivity is justified for cathodes with the high exchange current density i∗ORR . For sufficiently large i∗ORR the critical current density (19) largely exceeds working current density and the cathode operates in the linear mode, with the current-independent activation resistivity. It is worth mentioning that the literature data on i∗ORR vary by several orders of magnitude (see e.g. [18] and [24]), which reflects large sensitivity of this parameter to cathode composition and preparation technique. The conductivity of interconnect is high i.e., the cell electrodes are equipotential. Thus, the sum ˜ ˆj = E˜  ˜ +R

1

0

Since

∂˜c d˜z = − ∂˜z

1 0



1

ˆj d˜z 0

ˆj d˜z = ˆJ we obtain

(1 − c˜ 1 ) = 1

(49)

c˜ 1

where ≡ c˜ (1) is hydrogen concentration at the channel outlet. Solving (49) for c˜ 1 we find9 c˜ 1 = 1 −

1 

(51)

Eq. (51) relates outlet hydrogen concentration with the H2 stoichiometry factor . Clearly, the solution (48) must obey Eq. (51). Setting in (48) z˜ = 1 and equating the right side of the resulting equation to 1 − 1/ we get



E˜ =

8 Note the change in scaling for ˆj as compared to ˜j (4) in the previous sections. Note also that the following relation holds:

cref



ˆJ



−

 c˜ref +

 ˜ c˜ref + R



˜ +2 R

c˜ref

E˜ − ˆJ

2 =1−

1 

Solving this for E˜ we finally find the anode polarization curve

is the total voltage loss.

cH2

Cell polarization curve results from the following arguments. Integrating (43) over z˜ we find

1 ˜2 R

E˜ = V˜ oc − V˜ cell



3.3. Low current: polarization curve

(44)

is constant along the channel. Here

ˆj = ε˜j .

The shape of local current along z˜ is given by (46) with c˜ (48). Parameters E˜ and ˆJ are related by the polarization curve. In the ˜ ˆJ (see the next section) and hence the functions low-current limit E∼ c˜ (˜z ) and ˆj(˜z )/ˆJ do not depend on ˆJ (this is seen is we plug E˜ = kˆJ into Eqs. (48) and (46)). Thus, c˜ (˜z ) and ˆj(˜z )/ˆJ are the universal oneparametric functions of z˜ (Fig. 6).

(42)



˜ + 2  R





c˜ref

ˆJ

(52)

9 Note that Eq. (49) determines the relation of stoichiometry and utilization. By the definition, hydrogen utilization is u = 1 − c˜ 1 .Using this in (49) we get [25]

1 u= . 

(50)

6692

A.A. Kulikovsky / Electrochimica Acta 54 (2009) 6686–6695

Fig. 7. The dependence of anode resistivity

on hydrogen utilization u in the low√ current regime (the function 2 1 − 1 − u /u). When utilization is above 80%, the resistivity dramatically grows.

In the limiting cases of zero and 100% utilization, Eq. (58) yields Ra |u=0 = Ra0 Ra |u=1 = 2Ra0 We see, that 100% hydrogen utilization increases anode resistivity by a factor of 2. √

The function fu = 2u 1 − 1 − u is depicted in Fig. 7. Below u = 0.8 the growth of anode resistivity with u is moderate; however, above u = 0.8 this growth is dramatic. Running the cell at high utilization is, thus, detrimental for anode performance. Fig. 6. (a) The shapes of hydrogen concentration and (b) local current density along the anode channel for indicated values of hydrogen stoichiometry .

3.4. High-current regime: z-shapes and polarization curve In coordinates (41), Eq. (37) for the high-current polarization voltage takes the form

where





 = 1 −

1 

1−

(53)

We see that the total voltage loss E˜ is still proportional to the mean current density ˆJ . The coefficient of proportionality (total cell resis˜ cell ) is a sum of two terms tivity R ˜ +R ˜a ˜ cell = R R where ˜ a = 2  R

(54)



c˜ref

Ra = 2  Ra0 where

     b cref 0  R = 0 cH

2i∗ i

(57)

2

Ra =

u

1−



(59)

c˜



c˜ref

= E˜

c˜

(60)

Eqs. (43) and (60) determine the shapes of hydrogen concentration and local current density along the hydrogen channel. To obtain these shapes we solve (60) for ˆj: ˆj =

c˜ exp E˜ c˜ref

Using this in (43) and solving the resulting equation with the condition c˜ (0) = 1 we get





1 − u Ra0

ˆj



is anode resistivity at zero utilization (cf. Eq. (35)). In terms of hydrogen utilization u = 1/ Eq. (56) is

 2



2 ln

(56)

ˆj

(the subscript “1” is omitted). Hydrogen consumption is still governed by Eq. (43) and Eq. (44) completes the system. In this section, to simplify the analysis we will omit the resistive term in Eq. (44). This helps to understand the effect of hydrogen stoichiometry/utilization on the polarization curve without cumbersome calculations. ˜ = 0, Eq. (44) reduces to  ˜ or With R ˜ = E,

(55)

is the anode resistivity at finite . In dimension variables the latter relation takes the form

a

 ˜ = 2 ln

 c˜ref

(58)

c˜ =

1−

z˜



2

exp E˜

c˜ref ˆJ

2 (61)

A.A. Kulikovsky / Electrochimica Acta 54 (2009) 6686–6695

6693

Fig. 10. The critical current density as a function of temperature, Eq. (73). Fig. 8. Voltage loss due

to finite hydrogen utilization u in the high-current regime √ 2 ln 2(1 − 1 − u)/u . When utilization is above 80%, the voltage loss dramatically grows.

Cell polarization curve can now be obtained using Eq. (51). Setting in Eq. (61) z˜ = 1 and equating the result to 1 − 1/ we obtain

 1−



exp E˜



2

2 =1−

c˜ref ˆJ

1 

Solving this for E˜ we finally find



E˜ = 2 ln ˆJ



c˜ref

+ 2 ln (2  )

(62)

In terms of hydrogen utilization u = 1/ Eq. (62) takes the form



E˜ = 2 ln ˆJ



c˜ref

+ 2 ln

2  u

1−



1−u



(63)

We see that the effect of hydrogen utilization reduces to a constant shift of the polarization curve as a whole along the voltage axis. The value of this shift is given by the second term in Eq. (63). As in the low-current case, this shift dramatically increases above u  80% and reaches the value 2b ln 2  235 mV at u = 1 (Fig. 8, here we take b = 170 mV, Table C.2). The effect of hydrogen utilization on the anode polarization curve is summarized in Fig. 9. This figure shows the range of anode polarization voltage variation when hydrogen utilization varies from 0 to 100%.

4. Discussion Consider first the case of low hydrogen utilization. The key parameter which determines the regime of SOFC anode operation is the critical current density jcrit . If working current density j is less than jcrit , the anode operates in a linear regime with the anodic polarization proportional to j. If j  jcrit , the anode works in a non-linear regime, when polarization depends on current logarithmically, with the apparent Tafel slope being twice the kinetic value. In the intermediate region j  jcrit analytical solution to the problem discussed hardly exists. Semi-empirical relation (66) (Appendix C) enables to express jcrit as a function of temperature (Fig. 10). Typical working current density of SOFC is about 0.5 A cm−2 ; Fig. 10 shows that jcrit  0.5 A cm−2 is achieved at T  750 + 273 K. Therefore, running the cell above 800 + 273 K secures the inequality j < jcrit and the anode operation in the linear mode. At temperature below 700 + 273 K the anode works in the intermediate- or high-current regime. The most difficult situation arises when hydrogen utilization is large. In that case c˜H2 varies from 1 at the channel inlet to almost

zero at the outlet and hence the local value of the product jcrit cH2 /cref decreases toward the channel outlet. This may lead to a situation when domain at the inlet works in the linear mode (low-current), while domain close to the outlet operates in a nonlinear (high-current) mode. Analysis of this mixed regime requires numerical calculations. Figs. 8 and 9 show that in the low- and high-current regimes hydrogen utilization of 80% or above dramatically increases the activation voltage loss. It is worth mentioning that the model [21] also indicates a strong growth of anodic polarization voltage at hydrogen utilizations above 80%. 5. Conclusions

Fig. 9. The effect of hydrogen utilization on the anode polarization curve.  ˜ is the anode voltage loss, ˆJ is the dimensionless mean current density.

A quasi-2D model for SOFC anode performance is developed. The model couples variation of ionic current across the anode with the variation of hydrogen concentration along the anode channel. Solution to a through-plane model shows that there exists critical current density jcrit , which separates linear “low-current” regime of anode operation and non-linear “high-current” regime. If cell current j is lower than jcrit , the anode polarization voltage  linearly depends on j. If j  jcrit ,  depends on ˜j logarithmically with the apparent Tafel slope being twice the kinetic value. In the linear regime 100% hydrogen utilization doubles the anode activation resistivity. In the non-linear mode total hydrogen utilization shifts the polarization curve as a whole by  230 mV. In both cases the dependencies of voltage loss on utilization indi-

6694

A.A. Kulikovsky / Electrochimica Acta 54 (2009) 6686–6695

cate dramatic worsening of anode performance for utilization above 80%.

where jD =

0 2FDeff cH

2

∼ ˆ b cH2 cref E E∗ F hc J j j∗ jcrit i∗ i∗ORR L la l∗ R Ra Rcell T u Vcell Voc

va x z

Marks dimensionless variables Dimensionless current density in Section 3 Tafel slope (V) Hydrogen molar concentration (mol cm−3 ) Reference molar concentration (mol cm−3 ) Total voltage loss (V) Activation energy (J) Faraday constant Anode channel height (cm) Mean current density (A cm−2 ) Local proton current density in Section 2; local cell current density in the rest of the paper (A cm−2 ) Superficial exchange current density (A cm−2 ) Critical current density (A cm−2 ), Eq. (19) Volumetric exchange current density on the anode side (A cm−3 ) Volumetric exchange current density on the cathode side (A cm−3 ) Anode channel length (cm) Anode thickness (cm) Reaction penetration depth into the anode (cm), Eq. (9) Area-specific resistivity of all cell layers except anode ( cm2 ) Area-specific resistivity of the anode ( cm2 ) Area-specific resistivity of a cell ( cm2 ) Cell temperature (K) Hydrogen utilization Cell voltage (V) Open-circuit cell voltage (V) Flow velocity in the anode channel (cm s−1 ) Coordinate across the anode (cm) Coordinate along the channel (cm)

Superscripts 0 Inlet of the anode channel 1 Outlet of the anode channel

is the limiting current density due to H2 transport in the anode. Here Deff is the effective hydrogen diffusion coefficient in the anode. At 1073 K the value of Deff estimated with the relations [27] amounts to 4 cm2 s−1 (here anode porosity and tortuosity are taken to be 0.3 and 3, respectively). With the anode thickness la = 0.1 cm 0 = 10−5 mol cm−3 (100% hydrogen at and molar concentration cH 2

1072 K and atmospheric pressure), from (65) we get jD  77 A cm−2 . With b = 170 mV the concentration polarization at working current density j = 1 A cm−2 is conc = 2.2 mV. Calculations of Imperial College group [3] give for conc the value  20 mV. Unfortunately, not all parameters required for the calculation are specified in [3], so that the direct comparison with the estimate above is impossible. A detailed through-plane model of SOFC membrane-electrode assembly [11] gives for the limiting current density the value jD  7 A cm−2 , when the cell operates on 97%-H2 feed. Using this value in (64), at j = 1 A cm−2 we get conc  26 mV, which correlates with the data [3]. As seen, various sources and estimates show that the contribution of diffusion transport of H2 to the anode voltage loss can be safely ignored. Appendix C. Parameters for calculations were taken from the work of Aguiar et al. [3]. Superficial exchange current density reported in [3] has the form j∗ =



= i∗

cH2

(66)



exp

cref

x=la

j1 = j∗

 ˛F  1

RT

 (1 − ˛)F  1

− exp −

RT

(67)

cH2



exp

cref

 ˛F  1

RT

 (1 − ˛)F  1

− exp −

RT

(68)

Dividing (67) by (68) we get



∂ ln j  ∂x 

= x=la

i∗ j∗

(69)

We will assume that the anode works in low-current regime, when local ionic current varies exponentially with x (14). Calculating the derivative on the left side of (69), we obtain

Appendix B. Concentration polarization due to hydrogen diffusion through the anode thickness is given by [26] j jD





i∗ =





Exchange current density j∗ appears in equation [3]

Greek symbols ˛ Transfer coefficient ˇ Dimensionless parameter, Eq. (27)  Dimensionless parameter, Eq. (39)  Anode overpotential (V)  Hydrogen stoichiometry ε l∗ /la , small parameter (8) i Ionic conductivity of the anode ( −1 cm−1 )  Dimensionless function (53)

conc = −b ln 1 −



RT Ea ka exp − 2F RT

where the pre-exponential factor ka and the activation energy Ea are listed in Table C.1. The value of j∗ (A cm−2 ) is related to the volumetric exchange current density i∗ (A cm−3 ). The latter expresses the local electrochemical activity of the electrode, a number of charges consumed in unit volume per second. In contrast, j∗ is the integral characteristic, which includes the thickness of the current-generating domain in the anode. To show this we write Eq. (2) at the anode/electrolyte interface: ∂j  ∂x 

Subscripts a Anode H2 Hydrogen 1 Anode/electrolyte interface

(65)

la

Appendix A. Nomenclature



(64)

j∗ l∗

(70)

Table C.1 Anode exchange current density data [3]. ka ( −1 m−2 )

E∗ (kJ mol−1 )

6.54 × 10

140

11

A.A. Kulikovsky / Electrochimica Acta 54 (2009) 6686–6695

6695

Table C.2 Parameters for calculations. Here we take cH2 = cref . T (K)

i∗ (A cm−3 )

l∗ (cm) −3

700 + 273 800 + 273

8.47 × 10 4.52 × 10−3

9.86 102.1

i ( cm−1 ) −3

8.44 × 10 2.26 × 10−2

With (66) we finally get RT i∗ = 2F

k  a

l∗

 E  a

exp −

(71)

RT

Ionic conductivity of the anode was taken to be equal to the conductivity of bulk electrolyte and it was calculated according to [3]

 10,300 

i = 334 exp −

T

−1 cm−1

,

(72)

Parameters used for calculations are listed in Table C.2. Eqs. (71) and (9) allow us to exclude unknown parameter l∗ from the relations above. Using (9) in (71) and taking into account (19) we come to jcrit = i∗ =

RT ka F



cH2 cref

 E  ∗

exp −

 RT 2  c   2k2  H2 a 2F

cref

i b

(73)

RT

 2E  ∗

exp −

RT

(74)

These equations express the two basic characteristics of the cell, the characteristic current density jcrit (19) and the volumetric exchange current density i∗ as functions of temperature. Eqs. (73) and (74) allow us to calculate all the other parameters; the results of these calculations are listed in Table C.2. References [1] M. Mogensen, P.V. Hendriksen, Testing of electrodes, cells and short stacks, in: S.C. Singhal, K. Kendall (Eds.), High Temperature Solid Oxide Fuel

[2]

[3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27]

b (V)

jcrit (A cm−2 )

Ra ( cm−2 )

0.168 0.185

0.167 0.924

1.0 0.2

Cells. Fundamentals, Design and Applications, Elsevier, Amsterdam, 2003, p. 261. M. Eikerling, A.A. Kornyshev, A.A. Kulikovsky, Physical modeling of fuel cells and their components, in: A.J. Bard, M. Stratmann, D. Macdonald, P. Schmuki (Eds.), in: Encyclopedia of Electrochemistry, vol. 5, Wiley, New York, 2007, p. 447, ISBN 978-3-527-30397-7. P. Aguiar, C.S. Adjiman, N.P. Brandon, J. Power Sources 138 (2004) 120. H. Yakabe, T. Ogiwara, M. Hishinuma, I. Yasuda, J. Power Sources 101 (2001) 144. X. Xue, J. Tang, N. Sammes, Y. Du, J. Power Sources 142 (2005) 211. M. Iwata, T. Hikosaka, M. Morita, T. Iwanari, K. Ito, K. Onda, Y. Esaki, Y. Sakaki, S. Nagata, Solid State Ionics 132 (2000) 297. ´ D. Larrain, J. Van Herle, F. Marechal, D. Favrat, J. Power Sources 118 (2003) 367. P. Yuan, S.-F. Liu, Num. Heat Transf. 51 (2007) 941. Y. Wang, F. Yoshiba, T. Watanabe, S. Weng, J. Power Sources 170 (2007) 101. T.X. Ho, P. Kosinski, A.C. Hoffmann, A. Vik, Chem. Eng. Sci. 63 (2008) 5356. H. Zhu, R.J. Kee, J. Electrochem. Soc. 155 (2008) B715. G. Wang, P.P. Mukherjee, C.-Y. Wang, Electrochim. Acta 51 (2006) 3151. P.P. Mukherjee, C.-Y. Wang, J. Electrochem. Soc. 153 (2006) A840. Y. Suzue, N. Shikazono, N. Kasagi, J. Power Sources 184 (2008) 52. J. Park, X. Li, J. Power Sources 178 (2008) 248. C. Bao, N. Cai, AIChE J. 53 (2007) 2968. P. Costamagna, P. Costa, V. Antonucci, Electrochim. Acta 43 (1998) 375. L. Pisani, G. Murgia, J. Electrochem. Soc. 154 (2007) B793. P. Holtappels, L.G.J. de Haart, U. Stimming, J. Electrochem. Soc. 146 (1999) 1620. M.M. Hussain, X. Li, I. Dincer, Int. J. Energy Res. 29 (2005) 1083. P. Costamagna, K. Honegger, J. Electrochem. Soc. 145 (1998) 3995. A.J. Bard, L.R. Faulkner, Electrochemical Methods. Fundamentals and Applications, Wiley, New-York, 2001. F. Zhao, A.V. Virkar, J. Power Sources 141 (2005) 79. A. Kakac, A. Pramuanjaroenkij, X.Y. Zhou, Int. J. Hydrogen Energy 32 (2007) 761. A.A. Kulikovsky, J. Electrochem. Soc. 152 (2005) A1121. A.A. Kulikovsky, Electrochem. Commun. 4 (2002) 845. B. Todd, J.B. Young, J. Power Sources 110 (2002) 186.