- Email: [email protected]
Heat flux from the catalyst layer of a fuel cell
Electrochimica Acta 56 (2011) 9172–9179
Contents lists available at ScienceDirect
Electrochimica Acta journal homepage: www.elsevier.com/locate/electacta
Heat flux from the catalyst layer of a fuel cell A.A. Kulikovsky ∗,1,2 , J. McIntyre Institute of Energy and Climate Research – Fuel Cells (IEK-3), Research Centre “Jülich”, D-52425 Jülich, Germany
a r t i c l e
i n f o
Article history: Received 11 January 2011 Received in revised form 26 July 2011 Accepted 26 July 2011 Available online 5 August 2011 Keywords: Fuel cell Catalyst layer Model Heat transport
a b s t r a c t We report exact solutions to the problem of heat transport in the catalyst layer (CL) of a fuel cell. The solutions are obtained for the low- and high-current regimes of CL operation. The approximate equation for the heat flux from the CL valid for the whole range of current densities is suggested. This equation is suitable for CFD calculations of heat transport in cells and stacks. Heat fluxes from the catalyst layers of PEMFC, HT-PEMFC and DMFC are discussed. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction Key components of fuel cells are catalyst layers, where the half-cell electrochemical reactions produce or neutralize charged particles. Even in the best cells, these reactions generate quite a substantial amount of heat. This heat affects cell performance, since most of the transport and kinetic processes in a cell depend exponentially on temperature. Thus, calculation of the temperature field in cells and stacks is of great importance. Perhaps, the most advanced models for CL function have been developed for PEM fuel cells (see [1–3] and the references cited therein). However, the vast majority of CL models are isothermal and only few works have been devoted to the heat transport in CLs. The reason is that the temperature variation in the CL is usually not large. The temperature variation across the CL is, indeed, small (see below). However, the CL is thin and the heat flux from the CL (temperature gradient) is not small. This flux is obtained from the solution of the heat transport problem coupled to the problem for proton current j and reaction rate Rreac distributions across the CL (performance problem). Fig. 1 shows a schematic of the cathode catalyst layer environment. Proton current j0 enters the CL from the membrane. In the CL, the half-cell reaction converts the proton current into the electron current, so that at the CL/GDL interface, the proton current is
∗ Corresponding author. E-mail address: [email protected] (A.A. Kulikovsky). 1 Also with: Moscow State University, Research Computing Center, 119991 Moscow, Russia. 2 ISE member. 0013-4686/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.electacta.2011.07.113
zero and the electron current reaches the value of the cell current density j0 (Fig. 1). The conversion can occur in one of the two regimes (Fig. 2) [4]. In the low-current (LC) regime, the rate of conversion is distributed uniformly along x and the proton current linearly decays with the distance (dashed curves in Fig. 2). In the high-current (HC) regime, the conversion runs mainly in a small conversion domain close to the membrane interface, where the proton current and the rate of conversion rapidly decay (solid curves in Fig. 2). This situation arises in CLs with poor proton conductivity; the reaction rate has peak at the membrane interface due to lower expenditure for proton transport there. Importantly, in the HC regime, the polarization voltage is much larger than in the LC mode. An energy conservation equation is a necessary part of numerous CFD models of cells and stacks [3,5]. Most of CFD models consider catalyst layer as a thin interface generating heat and current. The usually used expression for the heat flux q emitted by the CL is
q=
TS nF
+ 0 j0 +
j02 lt ion
(1)
Here T is the temperature, S is the entropy change in the half-cell reaction, n is the number of electrons transferred, 0 is the halfcell polarization voltage, j0 is the cell current density, lt is the CL thickness and ion is the CL proton conductivity. The first term on the right side of Eq. (1) describes the thermodynamic (the term with S) and irreversible (the term with 0 ) heat generated in the electrochemical reaction. The last term in Eq.
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
F i* j jc j0 Kevap lt l* n q qtot qlow tot high
qtot
Rreac Rion S T T x
Tafel slope molar concentration of feed molecules (mol m−3 ) reference concentration of feed molecules (mol m−3 ) Faraday constant volumetric exchange current density (A m−3 ) local proton current density (A m−2 ) critical current density (A m−2 ), Eq. (10) cell current density (A m−2 ) evaporation rate constant (atm−1 s−1 ) catalyst layer thickness (m) thickness of the conversion domain (m) number of electrons transferred heat flux (W m−2 ) total heat flux from the CL (W m−2 ) total heat flux from the CL in the low-current regime (W m−2 ) total heat flux from the CL in the high-current regime (W m−2 ) rate of the half-cell reaction, A m−3 CL proton resistivity ( cm2 ) entropy change in the half-cell reaction (J mol−1 K−1 ) cell temperature (K) temperature drop across the CL (K) coordinate across the catalyst layer (m)
Electrolyte
Gas-diffusion layer
transfer coefficient dimensionless parameter (21) enthalpy of water evaporation (J mol−1 ) half-cell overpotential (V) overpotential due to proton resistivity (V) CL thermal conductivity (W m−1 K−1 ) dimensionless reaction penetration depth, Eq. (8) CL porosity density of liquid water (kg m−3 ) CL ionic conductivity (−1 m−1 )
Catalyst layer 0 0
500
0.5
0
0
0.2
0.4
0.6
0.8
1
0
Fig. 2. The two regimes of CL operation. In the low-current regime, the proton current decays along x linearly and the rate of conversion Rreac is constant (dashed lines). In the high-current regime, due to poor proton conductivity of the CL, the proton current and the rate of conversion rapidly decay near the membrane interface (solid lines). The conversion occurs in a narrow conversion domain of a thickness l* .
(1) represents the Joule heat due to the proton current in the CL.3 This expression implicitly assumes, that the proton current is uniform across the CL. This is rather crude approximation: in fact, in the LC regime, j linearly decays along x and the Joule term in Eq. (1) is overestimated by a factor of 3 [6]. What is the heat flux from the CL working in the HC regime, when j(x) is strongly non-linear (Fig. 2)? Below we show that in the HC regime • The Joule (second) term in Eq. (1) is j0 b, where b is the Tafel slope of the half-cell reaction. Note that this changes the dependence of heat flux on cell current from quadratic to linear. • The temperature drop in the CL due to the reaction heat is twice as large as in the LC regime.
Subscripts 0 at the membrane/CL interface 1 at the CL/GDL interface evap evaporation of water w water Greek ˛ ˇ Hevap ion ε εCL w ion
1000 1
ORR rate / A cm -3
b c1 cref
Current density / A cm-2
Nomenclature
9173
1
Fig. 1. Schematic of the generic catalyst layer and the shapes of the proton current density j, the electron current density je , the feed molecules concentration c and the overpotential .
Finally, we derive a general expression for the heat flux from the CL, which interpolates between the LC and HC regimes. This relation can be used in CFD calculations of cells and stacks. In this work, our main goal is to provide cell and stack modelers with the correct equation for the heat flux emitted by the catalyst layer. Other sources and sinks of heat could be added to the cell/stack energy balance equation, taking into account specific types of cells and regimes of operation. 2. Model 2.1. The basic equation In this work, the transport of feed molecules through the CL is assumed to be ideal. The case of strong transport limitations takes place in flooded cathodes of low-temperature cells; CL performance equations for this case are considered in [7,4]. Usually, in typical operating conditions, CLs provide good transport of feed molecules to the catalyst sites. The heat transport in the CL is governed by equation 2
−
∂ T = ∂x2
TS nF
+ Rreac +
j2 ion
(2)
Here is the CL thermal conductivity, x is the coordinate across the CL, is the local polarization voltage (overpotential), Rreac is the volumetric rate of electrochemical conversion (A m−3 ) and j is the local ionic (proton) current density (Fig. 1). Eq. (2) says that the
3 Note that the CL electronic conductivity is usually much larger than the proton one, and the respective Joule heat can be neglected.
9174
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
variation of conductive heat flux (the left side) equals the rate of heating due to reaction and Joule dissipation (the right side). Rreac determines also the rate of proton current decay along x:
∂j = −Rreac ∂x
(3)
The right side of Eq. (2) includes the dominant sources of heat in catalyst layers of low- and high-temperature fuel cells. On the cathode side of low-T cells, liquid water is formed in the oxygen reduction reaction (ORR). Evaporation of liquid water is an important mechanism of fuel cell cooling. However, the question of whether evaporation predominantly occurs in the CL or in the GDL is debated (cf. [1,8]). In the first case, the evaporation term should appear on the right side of Eq. (2). In the second case, this term appears in the GDL problem, while Eq. (2) remains unchanged. Here we assume that evaporation occurs mainly in the GDL; this assumption is supported by the estimate given in Appendix A. The exact solution to Eq. (2) is rather cumbersome [9]; however, this equation can be simplified. Temperature variation across the CL is not large; thus, on the right side we may safely replace T(x) by the temperature at the CL/GDL interface T1 ≡ T(lt ) (Fig. 1). In the LC regime, the variation of with x is small and we may set (x) 0 , where 0 ≡ (0) is the overpotential at the membrane interface (Fig. 1). In the HC regime, the reaction runs close to the membrane and thus only a small domain close to the membrane contributes to the heat production (see below). In this domain 0 and hence in all cases, in Eq. (2) we may replace (x) by 0 . Taking into account these remarks and Eq. (3), we come to 2
−
∂ T T1 S =− + 0 nF ∂x2
∂j ∂x
+
j2 ion
(4)
The heat crossover through the membrane is usually small ([10], Sections 3.4.6 and 3.4.7). Thus, to a good approximation the membrane can be considered as an ideal thermal insulator separating the anode and the cathode sides. The boundary conditions for Eq. (4), therefore, are
∂T ∂x
= 0,
T (lt ) = T1
(5)
x=0
where the second condition fixes the temperature at the CL/GDL interface. It should be emphasized that though Eq. (5) determine the temperature shape in the CL, they do not affect the equation for the total heat flux from the CL. Physically, the total heat flux is determined by the rate of heat energy production inside the CL and it does not depend on the boundary conditions. The adiabatic boundary condition at x = 0 means that the total heat flux qtot leaves the CL at the CL/GDL interface. This condition is used below to calculate qtot . The general solutions to the problem (4) depend on the shape of proton current density j(x); this shape results from the electrochemical behavior of the CL [10]. So far, the general solution to this problem has not been found. However, in the limiting cases of LC and HC regimes, the explicit expressions for the shapes j(x) have been derived [4]. These shapes allow us to solve the heat transport problem (4) and (5). Since the temperature variation across the CL is small, this variation does not affect Rreac . This allows us to decouple the thermal and electrochemical (performance) problems. 2.2. The regimes of CL operation The CL performance model used in this work is based on the standard macrohomogeneous model [11,7]. The flooded agglomerate model adds the pre-exponential effectiveness factor in the Butler–Volmer equation for the half-cell reaction rate [2]. This factor takes into account transport loss due to oxygen diffusion
within the agglomerate and for the reasonable set of parameters this factor weakly depends on the local overpotential. Therefore, the results of our work can be used for the catalyst layers with the flooded agglomerates, taking into account the effective lowering of the exchange current density due to oxygen transport in the agglomerates. Note also that the performance model used below is isothermal. This is justified since the temperature variation in the CL is small. The inequalities determining the regime of CL operation have been derived assuming that the rate of electrochemical conversion follows the Butler–Volmer relation [10]. The LC regime takes place when the cell current density j0 is below the critical value jc : j0 < jc
(6)
where
jc = and
2i∗ ion b
ε=
c1 cref
coth
1 ε
ion bcref 2i∗ lt2 c1
(7)
(8)
is the dimensionless Newman’s reaction penetration depth [12]. Here i* is the volumetric exchange current density (A cm−3 ), b is the Tafel slope, c1 is the feed molecules concentration in the CL and cref is their reference concentration. For air-fed cathode, cref is the oxygen molar concentration in air. Typically, the exchange current density is small (e.g., for oxygen reduction reaction) and ε is large. Physically it means, that the Newman’s reaction penetration depth greatly exceeds the CL thickness.4 For ε 1 we have coth (1/ε) ε and Eqs.(6) and (7) simplify to j0 <
ion b lt
(9)
jc =
ion b lt
(10)
Importantly, Eq. (9) does not contain the exchange current density i* . Physically, at large ε (small i* ), the shape of Rreac is controlled by the proton transport. Regardless of the value of i* , the reaction could be either uniformly distributed in the CL (small current, good proton transport), or “concentrated” at the membrane interface (large current, poor proton transport). The critical current density separates these two regimes. The high-current regime arises if the CL proton conductivity is low, or the CL is thick. For the case of Butler–Volmer kinetics, this regime takes place if j0 > 4 ion b/lt [10]. Our calculations show that for the heat transport problem this inequality should be strengthened to j0 >
8ion b lt
(11)
If the cell current density falls into the range ion b 8ion b < j0 < , lt lt
(12)
the CL operates in the transition region between the LC and HC regimes. In this region, the exact solution of performance problem does not exist. Fortunately, in some cases the LC and HC temperature shapes are close to each other, so that the transition-region shape is “trapped” within a narrow domain between the two analytical solutions (see below). For the catalyst layers considered in
4 The inverse relation takes place on the anode side of hydrogen cells, where the exchange current density is very large.
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
9175
Table 1 Parameters for calculations. (am) – assumed, (ct) – calculated. HT-PEMFC Transfer coefficient ˛ Exchange current density i* (A m−3 ) Catalyst layer thickness lt (m) Overpotential 0 (V) Critical current density jc (A m−2 ) Thickness of the conversion domain l* (m) Parameter ˇ Proton conductivity of the electrolyte phase, p (S m−1 ) ORR entropy change, S (J mol−1 K−1 ) Thermal conductivity of CL, (W m−1 K−1 ) Current density j0 (A m−2 ) Liquid water density w (kg m−3 ) Molecular weight of water Mw (kg mol−1 ) Heat capacity of liquid water cpw (J kg−1 K−1 ) Enthalpy of water evaporation Hevap (J mol−1 ) Evaporation rate constant Kevap (atm−1 s−1 ) CL porosity εCL Liquid saturation in the CL s
PEMFC
0.777 [16] 103 [16] 1.5 × 10−4 [16] 0.665 (ct) 320 5 × 10−6 (ct) 2.858 (ct)
1.0 [13] 103 (am) 10−5 (am) 0.427 (ct) 2955 – – 1 [17] 326.36 [18] 0.27 [19] 104 (am) 103 0.018 4190 41.7 × 103 1 [15] 0.3 (am) 0.3 (am)
this work, the values of critical current density (10) are depicted in Table 1.
Setting here x = 0 we obtain the temperature drop T ≡ T0 − T1 across the CL in the LC regime:
2.3. Low-current regime: the temperature shapes
T =
In the LC regime, the proton current linearly decays with x (Fig. 2):
2.4. High-current regime: the temperature shapes
j(x) = j0 1 −
x lt
(13)
With this relation, Eq. (4) transforms to 2
∂ T − = ∂x2
T S 1 nF
+ 0
j
0
+
lt
j02 ion
x 1− lt
−
∂T = ∂x
T S 1 nF
+ 0
j x 0
+
lt
j02 lt
2
3ion
1− 1−
T S 1 nF
+ 0 j0 +
qlow tot =
nF
x lt
3
j02 lt
(15)
(16)
3ion
(17)
is the proton resistive overpotential and lt ion
is the CL proton resistivity. Integrating Eq. (15) from x to lt we find the temperature profile across the CL:
T (x) = T1 +
+
ˇ 2
ˇ
1−
x lt
(20)
j0 2jc + j0
(21)
and jc is given by (7). Note
that Eq. (21) is an approximate solution to equation j0 = jc ˇ tan ˇ/2 , which results from Eq. (20) with x = 0 [4]. Fig. 2 shows that the electrochemical conversion occurs in a narrow conversion domain at the membrane interface. The characteristic width of this domain is [4] ion b j0
(22)
In this domain, the square of proton current density can be approximated by exponent [4] j2 j02 exp −
j0 Rion = 3
Rion =
(19)
2
x
where ion
j l 0 t
In the HC regime, the proton current decays along x according to [4]
l∗ =
+ 0 + ion j0
3ion 2
where
It is advisable to rewrite this equation in the form
T S 1
+ 0 +
(14)
Setting here x = lt we get the total heat flux leaving the CL [6]: qlow tot =
nF
j = jc ˇ tan
Integrating this equation from 0 to x and taking into account the first of Eq. (5), we obtain
T S 1
T S 1 nF
j02 lt2 3ion
+ 0
1−
j l 0 t
2
1−
x
1
lt
4
1−
x2 lt2
1−
x 3 lt
Note that Eq. (23) approximates the square of the exact expression (20) only close to the membrane (Fig. 2). However, this is exactly what is needed for the heat transport problem: outside the conversion domain, the proton current and the reaction rate are small and this region does not contribute to heat production (Fig. 2). Using Eq. (23) in (4) we get an equation for CL temperature in the HC regime: 2
(23)
l∗
−
∂ T T1 S =− + 0 nF ∂x2
∂j ∂x
+
j02 ion
x
exp −
(24)
l∗
Integrating this equation from 0 to x and taking into account (5), we obtain (18)
−
∂T = ∂x
T S 1 nF
+ 0 (j0 − j) +
j02 l∗ ion
x
1 − exp −
l∗
(25)
9176
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
which describes the linear temperature shape across the CL. Though this linear shape does not satisfy the adiabatic boundary condition at x = 0, it is close to the exact result, Eq. (28) (Fig. 3). Setting in (29) x = 0, we get the temperature drop across the CL:
166 Thigh
164 Tlow, no Joule
T
GDL
Temperature / oC
Thigh, no Joule
Tlow
162
160
0
0.2
0.4
0.6
0.8
1
Fig. 3. Temperature distribution across the 150-m HT-PEMFC cathode catalyst layer calculated with the low- and high-current equations. The low-current curves are denoted as Tlow , the high-current as Thigh . Solid Thigh curve: the exact solution (28), long-dashed curve: the approximate linear solution (29). Parameters for all curves are listed in Table 1.
Setting x = lt here we find high
qtot =
T S 1 nF
+ 0 j0 +
j02 l∗
ion
l t
1 − exp −
(26’)
l∗
In the HC regime, lt /l* 1 and the exponent exp (− lt /l* ) is vanishingly small. Neglecting this exponent, we find the heat flux leaving the CL: high
qtot
T S 1 nF
+ 0 j0 +
j02 l∗
(26)
ion
Rearranging terms and taking into account (22), we come to high
qtot =
T S 1 nF
+ 0 + b j0
(27)
We see that in the HC regime, the Joule heating is equivalent to increase in the half-cell overpotential 0 by the Tafel slope b. In low-temperature fuel cells, b is in the order of 50–100 mV [13], whereas 0 is in the order of 300–500 mV. Thus, in the HC regime, the contribution of the Joule heating to the overall heat flux does not exceed 30% (see also below). In high-temperature solid oxide fuel cells (SOFCs), b 150 − 200 mV [14] and the Joule contribution to the heat flux may exceed 50%. Integrating (25) from x to lt , we get the temperature shape across the CL:
T (x) = T1 +
T S 1 nF
− exp −
+ 0
x
2l∗
x
− exp −
j l 0 t +
1−
x 2l∗ + lt lt
j0 lt b x l∗ 1− + lt lt
exp −
lt 2l∗
l t
exp −
l∗
The exponential terms in Eq. (28) contain a small pre-factor l* /lt 1. These exponents “flatten” the shape T(x) near the membrane to provide the adiabatic boundary condition at the membrane interface (Fig. 3). However, since l* lt , this flattening only occurs in a small conversion domain (Fig. 3). Therefore, to a good approximation exponents in Eq. (28) can be ignored and we come to a simple result T (x) T1 +
T S 1 nF
j l 0 t
+ 0 + b
x 1− lt
nF
j l 0 t
+ 0 + b
(30)
Comparing this to Eq. (19) we see that in the HC regime, the temperature drop due to the reaction is twice as large as in the LC regime. This result could have been foreseen. Ignoring for a moment the Joule heating we note, that the total heat flux from the CL does not depend on the distribution of reaction rate in the CL. Indeed, the total amount of heat energy released in the reaction events per unit time is independent of their spatial distribution. Further, in the HC regime, the temperature shape is close to linear, since in most of the CL thickness, Rreac and j(x) are close to zero. We, therefore, come to a simple geometric problem. The LC temperature shape is parabolic (the Joule term is ignored): Tlow (x) = T1 + a2 (1 − (x2 /lt )2 ), and the HC shape is linear: Thigh (x) =
T1 + c 2 1 − lx , where a and c are constants. Both these shapes t must have the same slope at x = lt ; elementary calculation shows that c2 = 2a2 and hence T|high = 2T|low . 2.5. The general expression for the heat flux Substituting (22) into (26’) and rearranging terms, we come to high
qtot =
T S 1 nF
+ 0 + b 1 − exp −
j0 lt ion b
j0
(31)
Eq. (31) is the exact expression for the HC heat flux from the CL. As discussed above, this expression is simplified to Eq. (27), since in the HC regime, the parameter j0 lt /( ion b) is small. Eq. (31) suggests the general expression for the heat flux from the CL, valid in the whole range of cell current density: qtot =
T S 1 nF
+ 0 + b 1 − exp −
j0 lt 3ion b
j0
(32)
Eq. (32) differs from Eq. (31) by a factor of 3 in the power of exponent. If the LC condition Eq. (9) is fulfilled, the exponent in (32) can be expanded in series; retaining the linear term, we come to the LC Eq. (15). In the HC limit of Eq. (11), the exponent in Eq. (32) can be neglected, and we come to Eq. (27). Thus, in the LC and HC limits, Eq. (32) reduces to the respective exact expression and hence this equation interpolates between the limiting solutions. Note that in Nafion-based cells, ion in Eq. (32) depends on the water content of the membrane phase in the CL. Note also that Eq. (32) does not take into account the situations with the gradient of water content through the CL. 2.6. DMFC cathode: heat flux due to fuel crossover
(28)
l∗
T S 1
(29)
A feature of a direct methanol fuel cell is large fuel crossover through the membrane. Permeated methanol is catalytically burned on the cathode side and gives a significant (sometimes dominating) contribution to the heat flux from the cathode catalyst layer (CCL). Note that methanol comes to the CCL at a temperature close to T1 and hence the methanol flux itself does not change the heat balance in the CCL. The main effect of methanol is heat released in the methanol–oxygen catalytic combustion. The respective heat flux is given by equation [10] qcross
T1 Scross = 6F
ˇ 1+ˇ
(jlim − j0 )
(33)
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
Table 2 The heat flux (kW m−2 ) from the catalyst layer calculated with the low-, high- and general relations (the rows LC, HC and General, respectively). The row “Regime” indicates the regime of CL operation. LC and HC stand for the low- and high-current regimes, respectively, T indicates the transition region between the LC and HC. The row “Reaction” shows the heat flux from the half-cell reaction without the Joule term. The row “Joule” shows the Joule heat flux. The last row indicates the Joule heat flux calculated using the conventional low-current approximation, Eq. (17). The parameters for calculation are listed in Table 1.
Low-current
10
High-current Fitting
j0 (A cm−2 ) Real regime LC equation, Eq. (17) HC equation, Eq.(27) General Eq. (32) Reaction Joule LC Joule
0 0
0.2
0.4
0.6
Current density / A
0.8
1
cm-2
Fig. 4. Heat flux from the HT-PEFC cathode calculated with the low-current Eq. (16) (top curve), high-current Eq. (27) and approximate fitting Eq. (32). The latter two curves are indistinguishable. 0 /298 is the entropy change in the Here Scross = HM 0 is the respective methanol–oxygen combustion reaction, HM enthalpy change, and
ˇ=
Dm lb Db lm
is the crossover parameter. Here Dm , Db are the diffusion coefficients of methanol in the membrane and anode backing layer (ABL), respectively, lm and lb are the membrane and ABL thickness, respectively. The limiting current density jlim is given by jlim =
6FDb cM lb
(34)
where cM is methanol molar concentration in the feed channel. In DMFC cathode, the flux (33) should be added to the right side of Eq. (32). Note that at small cell currents, the heat flux due to methanol crossover, Eq. (33) greatly exceeds all the other heat sources in DMFC. In the opposite limit of j0 jlim this flux is small and the heat sources due to ORR and Joule heating should be taken into account. 3. Numerical examples 3.1. HT-PEMFC cathode Parameters for the calculations in this section are listed in Table 1. For the CCL of HT-PEMFC, the working current density greatly exceeds 8jc (Table 1), i.e., the CCL operates in the HC regime. Nonetheless, for comparison, we also display the LC curves calculated with the same set of parameters. These curves show the temperature shapes which follow from the misconception about the regime of CL operation. Fig. 3 shows the LC and HC shapes T(x). The temperature shapes without the Joule term (the last term in Eqs. (18) and (28) are also displayed. The LC curves exhibit quite significant contribution of Joule heating to the temperature shape (Fig. 3). However, the HC curves show that in fact this contribution is small (Fig. 3). Physically, in the HC regime, the proton current rapidly disappears close to the membrane and the contribution of Joule heat is small. If the CL were operated in the LC regime, the proton current would linearly decay with the distance and fraction of the Joule heating in the overall heat flux would be much larger. This effect also is illustrated in Fig. 4, which shows the heat fluxes from the HT-PEFC cathode calculated with the low-, high- and fitting Eqs. (16), (27) and (32), respectively. As can be seen, as the CCL
HT-PEMFC
PEMFC
1.0 HC 15.31 10.79 10.79 10.31 0.48 5.0
1.0 T 7.50 7.47 7.37 7.17 0.20 0.33
operates in the HC regime, the LC Eq. (15) strongly overestimates the Joule heat flux, whereas Eqs. (27) and (32) give correct and practically indistinguishable results. The numerical data are listed in the first column of Table 2. Note that at small current densities, the contribution of Joule heating is small and all the curves in Fig. 4 coincide. Note also that the curves in Fig. 4 are plotted taking into account the dependence 0 (j0 ) (Eq. (74) in [4]). 3.2. PEMFC cathode The curves for PEMFC cathode are shown in Fig. 5. The pair j0 , jc shows that this cathode works in the transition region (jc < j0 < 8jc , Table 1). However, the gap between the LC and HC curves is small and hence both these solutions are close to the real curve. The PEMFC working current density is large (1 A cm−2 , Table 1); however, the temperature drop in the CCL does not exceed 0.2 K. The main reason is that the CCL of PEMFC is only 10 m thick, 10 times smaller, than in the HT-PEMFC (Table 1). In spite of the small temperature drop, the heat flux from the CCL is about 8 kW m−2 , which is comparable to that flux from the CCL of HTPEMFC (11 kW m−2 , Table 2). 4. Discussion Figs. 3 and 5 demonstrate the difference between the temperature profiles when the CL works in the low and high-current regimes. These shapes do not directly affect the electrochemical performance of the CL; however, they determine the heat flux leaving the CL. If the heat contact between the GDL and current collector
70.2
Temperature / oC
Heat flux / kW m-2
HT-PEFC
70.15
Thigh
Thigh, no Joule
Tlow
70.1
GDL
15
5
9177
Tlow, no Joule
70.05
70
0
0.2
0.4
0.6
0.8
1
Fig. 5. (Color online) Same as in Fig. 3 for the PEMFC cathode. Parameters for all curves are listed in Table 1.
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
is poor, part of this heat flux would be reflected and returned to the CL, thereby heating up the CL and accelerating the reaction. It should be emphasized that the total heat flux leaving the CL does not depend on boundary conditions: this flux is determined solely by the rate of heat production in the CL. Eq. (5) only affect the temperature shape across the CL. Note that Eq. (32) does not specify the direction of the heat flux. On the cathode side of low-T cells, there is a heat flux transported with the liquid water produced in the ORR. This heat flux is implicitly taken into account on the right side of Eq. (32). Physically, ORR generates water molecules in liquid state and part of the reaction heat is spent to warm up these molecules to the working temperature T. The heat flux transported with the flow of liquid water produced in the ORR is qw = Mw cpw Tj0 /(2F). With the data from Table 1 we get qw 1.34 kW m−2 . Comparing this to the total heat flux of 7.5 kW m−2 (Table 2), we see that the liquid water transports about 20% of the total heat produced in the CL.5 Note also that the contribution of water flux through the membrane to the CL heat balance is negligible, since the MEA sandwich is almost isothermal. Fig. 4 and Table 2 indicate that the general equation for the heat flux (32) has accuracy better than 2%. Note that in the case of HTPEMFC, the standard LC Eq. (17) overestimates the total heat flux by 50% (Table 2) due to incorrect (10 times higher) value of the Joule flux. Eq. (32) is recommended for use in CFD calculations. If local current over the cell surface is non-uniform, Eq. (32) allows the user not to care about the regime of CL operation. This equation returns the correct value of the heat flux from domains with low, intermediate and high local current densities. The model discussed ignores the voltage loss due to poor feed molecules transport in the CL. If this loss is high, the reaction rate Rreac has maximum at the CL/GDL interface [4]. In that case, the conversion occurs in a narrow conversion domain at the CL/GDL interface; the thickness of this domain depends on the oxygen diffusion coefficient in the CL. Qualitatively, in this regime, the drop of temperature over the CL thickness would be much smaller, since heat would be effectively removed to the GDL.
5. Conclusions Solutions to the CL performance problem are used to derive expressions for the temperature shape and for the heat flux from the catalyst layer operating in the low- and high-current regimes. In the high-current regime, the contribution of the Joule heating to the total heat flux appears to be much smaller, than in the low-current mode. A simple equation for the heat flux from the CL valid in the whole range of current densities is suggested. The analytical results are illustrated with the numerical examples for the catalyst layers of PEMFC and HT-PEMFC.
Acknowledgement This work has been partly supported by the German Federal Ministry of Economics and Technology, Project No. 0327853A.
5 To be precise, the convective term for the transport of liquid water should be added to the left side of Eq. (2). This would lower the conductive heat flux by the value qw and the same value would appear as the convective heat flux. Note that the total heat flux of Eq. (32) would not change, since this flux is determined by the rate of heat production in the CL only.
1.5
Evaporation heat flux / kW m-2
9178
1
0.5
0
20
40
60
80
100
Temperature / oC Fig. 6. Absolute value of the heat flux due to liquid water evaporation in the catalyst layer, Eq. (35). The curve is an upper estimate with the data from Table 1 and pw = 0, which corresponds to zero pressure of water vapor in the CL.
Appendix A. Heat flux due to water evaporation The heat flux due to liquid water evaporation in the CL is obtained if we multiply the molar transfer rate due to evaporation [15] by the respective enthalpy change Hevap . This yields
qevap = −Kevap εCL slt
Hevap w Mw
(psat w − pw )
(35)
where Kevap is the evaporation rate constant, εCL is the CL porosity, s is the CL liquid saturation, Hevap is the enthalpy of water evaporation, w is the liquid water density, Mw is the molecular weight of water, psat w is the pressure of saturated water vapor and pw is the partial pressure of water vapor in the CL. Physically, the product εCL s is the fraction of CL volume occupied by liquid water, Hevap w /Mw is the energy consumed for evaporation of unit volume of water, Kevap (psat w − pw ) is the “driving force” for evaporation, which is proportional to the difference of saturated and available water vapor pressures. The constant Kevap in an efficient manner takes into account the kinetics of water phase transformation in the porous media of the CL. The pressure of saturated water vapor is given by ˜ (T ) psat w = p0 p where p0 = 1 atm and
p˜ (T ) = exp a0 + a1 (T − 273) + a2 (T − 273)2 + a3 (T − 273)3
(36)
(37)
describes the temperature dependence (the coefficients a0 , . . . a3 are given in Table 3). In models which assume that evaporation occurs in the CL, the (negative) heat flux qevap should be added to the right side of Eq. (32). Written in the form of Eq. (35), the flux qevap is determined by the CL liquid saturation, temperature and density of water vapor in the CL and it does not explicitly depend on the cell current density. However, the dependence on j0 is “hidden” in the liquid saturation s, which is a function of current-dependent water transport in the MEA. The absolute value of heat flux due to water evaporation in the CL calculated with the data from Table 1 is depicted in Fig. 6. As can be seen, at typical working temperatures below 70◦ C, this heat flux does not exceed 1 kW m−2 (cf. the data in Table 2). In PEMFC, the GDL is more than ten times thicker than the CL, and the GDL porosity typically is twice that of the CL. Thus, the heat flux due to evaporation in the GDL is at least an order of magnitude larger and we may conclude that water evaporation occurs predominantly in the GDL.
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
9179
Table 3 The coefficients in Eq. (37). a0
a1
a2
a3
−2.1794 ln (10)
0.02953 ln (10)
−9.1837 × 10−5 ln (10)
1.4454 × 10−7 ln (10)
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
M. Eikerling, J. Electrochem. Soc. 153 (2006) E58. D. Harvey, J.G. Pharoah, K. Karan, J. Power Sources 179 (2008) 209. A. Weber, J. Newman, Chem. Rev. 104 (2004) 4679. A.A. Kulikovsky, Electrochim. Acta 55 (2010) 6391. C.-Y. Wang, Chem. Rev. 104 (2004) 4727. A.A. Kulikovsky, J. Power Sources 162 (2006) 1236. M. Eikerling, A.A. Kornyshev, J. Electroanal. Chem. 453 (1998) 89. A.Z. Weber, J. Newman, J. Electrochem. Soc. 153 (2006) A2205. A.A. Kulikovsky, Electrochem. Commun. 9 (2007) 6. A.A. Kulikovsky, in: Analytical Modelling of Fuel Cells, Elsevier, Amsterdam, 2010.
[11] M.L. Perry, J. Newman, E.J. Cairns, J. Electrochem. Soc. 145 (1998) 5. [12] J. Newman, in: Electrochemical Systems, Prentice Hall, Inc., Englewood Cliffs, NJ 07632, 1991. [13] K.C. Neyerlin, W. Gu, J. Jorne, H. Gasteiger, J. Electrochem. Soc. 153 (2006) A1955. [14] V.A. Danilov, M.O. Tade, Int. J. Hydrogen Energy 34 (2009) 6876. [15] D. Natarajan, T.V. Nguyen, J. Electrochem. Soc. 148 (2001) A1324. [16] A.A. Kulikovsky, H.F. Oetjen, C. Wannek, Fuel Cells 10 (2010) 363, doi:10.1002/fuce.200900135. [17] Z. Xie, T. Navessin, K. Shi, R. Chow, Q. Wang, D. Song, B. Andreaus, M. Eikerling, Z. Liu, S. Holdcroft, J. Electrochem. Soc. 152 (2005) A1171. [18] M.J. Lampinen, M. Fomino, J. Electrochem. Soc. 140 (1993) 3537. [19] M. Khandelwal, M. Mench, J. Power Sources 161 (2006) 1106.
Contents lists available at ScienceDirect
Electrochimica Acta journal homepage: www.elsevier.com/locate/electacta
Heat flux from the catalyst layer of a fuel cell A.A. Kulikovsky ∗,1,2 , J. McIntyre Institute of Energy and Climate Research – Fuel Cells (IEK-3), Research Centre “Jülich”, D-52425 Jülich, Germany
a r t i c l e
i n f o
Article history: Received 11 January 2011 Received in revised form 26 July 2011 Accepted 26 July 2011 Available online 5 August 2011 Keywords: Fuel cell Catalyst layer Model Heat transport
a b s t r a c t We report exact solutions to the problem of heat transport in the catalyst layer (CL) of a fuel cell. The solutions are obtained for the low- and high-current regimes of CL operation. The approximate equation for the heat flux from the CL valid for the whole range of current densities is suggested. This equation is suitable for CFD calculations of heat transport in cells and stacks. Heat fluxes from the catalyst layers of PEMFC, HT-PEMFC and DMFC are discussed. © 2011 Elsevier Ltd. All rights reserved.
1. Introduction Key components of fuel cells are catalyst layers, where the half-cell electrochemical reactions produce or neutralize charged particles. Even in the best cells, these reactions generate quite a substantial amount of heat. This heat affects cell performance, since most of the transport and kinetic processes in a cell depend exponentially on temperature. Thus, calculation of the temperature field in cells and stacks is of great importance. Perhaps, the most advanced models for CL function have been developed for PEM fuel cells (see [1–3] and the references cited therein). However, the vast majority of CL models are isothermal and only few works have been devoted to the heat transport in CLs. The reason is that the temperature variation in the CL is usually not large. The temperature variation across the CL is, indeed, small (see below). However, the CL is thin and the heat flux from the CL (temperature gradient) is not small. This flux is obtained from the solution of the heat transport problem coupled to the problem for proton current j and reaction rate Rreac distributions across the CL (performance problem). Fig. 1 shows a schematic of the cathode catalyst layer environment. Proton current j0 enters the CL from the membrane. In the CL, the half-cell reaction converts the proton current into the electron current, so that at the CL/GDL interface, the proton current is
∗ Corresponding author. E-mail address: [email protected] (A.A. Kulikovsky). 1 Also with: Moscow State University, Research Computing Center, 119991 Moscow, Russia. 2 ISE member. 0013-4686/$ – see front matter © 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.electacta.2011.07.113
zero and the electron current reaches the value of the cell current density j0 (Fig. 1). The conversion can occur in one of the two regimes (Fig. 2) [4]. In the low-current (LC) regime, the rate of conversion is distributed uniformly along x and the proton current linearly decays with the distance (dashed curves in Fig. 2). In the high-current (HC) regime, the conversion runs mainly in a small conversion domain close to the membrane interface, where the proton current and the rate of conversion rapidly decay (solid curves in Fig. 2). This situation arises in CLs with poor proton conductivity; the reaction rate has peak at the membrane interface due to lower expenditure for proton transport there. Importantly, in the HC regime, the polarization voltage is much larger than in the LC mode. An energy conservation equation is a necessary part of numerous CFD models of cells and stacks [3,5]. Most of CFD models consider catalyst layer as a thin interface generating heat and current. The usually used expression for the heat flux q emitted by the CL is
q=
TS nF
+ 0 j0 +
j02 lt ion
(1)
Here T is the temperature, S is the entropy change in the half-cell reaction, n is the number of electrons transferred, 0 is the halfcell polarization voltage, j0 is the cell current density, lt is the CL thickness and ion is the CL proton conductivity. The first term on the right side of Eq. (1) describes the thermodynamic (the term with S) and irreversible (the term with 0 ) heat generated in the electrochemical reaction. The last term in Eq.
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
F i* j jc j0 Kevap lt l* n q qtot qlow tot high
qtot
Rreac Rion S T T x
Tafel slope molar concentration of feed molecules (mol m−3 ) reference concentration of feed molecules (mol m−3 ) Faraday constant volumetric exchange current density (A m−3 ) local proton current density (A m−2 ) critical current density (A m−2 ), Eq. (10) cell current density (A m−2 ) evaporation rate constant (atm−1 s−1 ) catalyst layer thickness (m) thickness of the conversion domain (m) number of electrons transferred heat flux (W m−2 ) total heat flux from the CL (W m−2 ) total heat flux from the CL in the low-current regime (W m−2 ) total heat flux from the CL in the high-current regime (W m−2 ) rate of the half-cell reaction, A m−3 CL proton resistivity ( cm2 ) entropy change in the half-cell reaction (J mol−1 K−1 ) cell temperature (K) temperature drop across the CL (K) coordinate across the catalyst layer (m)
Electrolyte
Gas-diffusion layer
transfer coefficient dimensionless parameter (21) enthalpy of water evaporation (J mol−1 ) half-cell overpotential (V) overpotential due to proton resistivity (V) CL thermal conductivity (W m−1 K−1 ) dimensionless reaction penetration depth, Eq. (8) CL porosity density of liquid water (kg m−3 ) CL ionic conductivity (−1 m−1 )
Catalyst layer 0 0
500
0.5
0
0
0.2
0.4
0.6
0.8
1
0
Fig. 2. The two regimes of CL operation. In the low-current regime, the proton current decays along x linearly and the rate of conversion Rreac is constant (dashed lines). In the high-current regime, due to poor proton conductivity of the CL, the proton current and the rate of conversion rapidly decay near the membrane interface (solid lines). The conversion occurs in a narrow conversion domain of a thickness l* .
(1) represents the Joule heat due to the proton current in the CL.3 This expression implicitly assumes, that the proton current is uniform across the CL. This is rather crude approximation: in fact, in the LC regime, j linearly decays along x and the Joule term in Eq. (1) is overestimated by a factor of 3 [6]. What is the heat flux from the CL working in the HC regime, when j(x) is strongly non-linear (Fig. 2)? Below we show that in the HC regime • The Joule (second) term in Eq. (1) is j0 b, where b is the Tafel slope of the half-cell reaction. Note that this changes the dependence of heat flux on cell current from quadratic to linear. • The temperature drop in the CL due to the reaction heat is twice as large as in the LC regime.
Subscripts 0 at the membrane/CL interface 1 at the CL/GDL interface evap evaporation of water w water Greek ˛ ˇ Hevap ion ε εCL w ion
1000 1
ORR rate / A cm -3
b c1 cref
Current density / A cm-2
Nomenclature
9173
1
Fig. 1. Schematic of the generic catalyst layer and the shapes of the proton current density j, the electron current density je , the feed molecules concentration c and the overpotential .
Finally, we derive a general expression for the heat flux from the CL, which interpolates between the LC and HC regimes. This relation can be used in CFD calculations of cells and stacks. In this work, our main goal is to provide cell and stack modelers with the correct equation for the heat flux emitted by the catalyst layer. Other sources and sinks of heat could be added to the cell/stack energy balance equation, taking into account specific types of cells and regimes of operation. 2. Model 2.1. The basic equation In this work, the transport of feed molecules through the CL is assumed to be ideal. The case of strong transport limitations takes place in flooded cathodes of low-temperature cells; CL performance equations for this case are considered in [7,4]. Usually, in typical operating conditions, CLs provide good transport of feed molecules to the catalyst sites. The heat transport in the CL is governed by equation 2
−
∂ T = ∂x2
TS nF
+ Rreac +
j2 ion
(2)
Here is the CL thermal conductivity, x is the coordinate across the CL, is the local polarization voltage (overpotential), Rreac is the volumetric rate of electrochemical conversion (A m−3 ) and j is the local ionic (proton) current density (Fig. 1). Eq. (2) says that the
3 Note that the CL electronic conductivity is usually much larger than the proton one, and the respective Joule heat can be neglected.
9174
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
variation of conductive heat flux (the left side) equals the rate of heating due to reaction and Joule dissipation (the right side). Rreac determines also the rate of proton current decay along x:
∂j = −Rreac ∂x
(3)
The right side of Eq. (2) includes the dominant sources of heat in catalyst layers of low- and high-temperature fuel cells. On the cathode side of low-T cells, liquid water is formed in the oxygen reduction reaction (ORR). Evaporation of liquid water is an important mechanism of fuel cell cooling. However, the question of whether evaporation predominantly occurs in the CL or in the GDL is debated (cf. [1,8]). In the first case, the evaporation term should appear on the right side of Eq. (2). In the second case, this term appears in the GDL problem, while Eq. (2) remains unchanged. Here we assume that evaporation occurs mainly in the GDL; this assumption is supported by the estimate given in Appendix A. The exact solution to Eq. (2) is rather cumbersome [9]; however, this equation can be simplified. Temperature variation across the CL is not large; thus, on the right side we may safely replace T(x) by the temperature at the CL/GDL interface T1 ≡ T(lt ) (Fig. 1). In the LC regime, the variation of with x is small and we may set (x) 0 , where 0 ≡ (0) is the overpotential at the membrane interface (Fig. 1). In the HC regime, the reaction runs close to the membrane and thus only a small domain close to the membrane contributes to the heat production (see below). In this domain 0 and hence in all cases, in Eq. (2) we may replace (x) by 0 . Taking into account these remarks and Eq. (3), we come to 2
−
∂ T T1 S =− + 0 nF ∂x2
∂j ∂x
+
j2 ion
(4)
The heat crossover through the membrane is usually small ([10], Sections 3.4.6 and 3.4.7). Thus, to a good approximation the membrane can be considered as an ideal thermal insulator separating the anode and the cathode sides. The boundary conditions for Eq. (4), therefore, are
∂T ∂x
= 0,
T (lt ) = T1
(5)
x=0
where the second condition fixes the temperature at the CL/GDL interface. It should be emphasized that though Eq. (5) determine the temperature shape in the CL, they do not affect the equation for the total heat flux from the CL. Physically, the total heat flux is determined by the rate of heat energy production inside the CL and it does not depend on the boundary conditions. The adiabatic boundary condition at x = 0 means that the total heat flux qtot leaves the CL at the CL/GDL interface. This condition is used below to calculate qtot . The general solutions to the problem (4) depend on the shape of proton current density j(x); this shape results from the electrochemical behavior of the CL [10]. So far, the general solution to this problem has not been found. However, in the limiting cases of LC and HC regimes, the explicit expressions for the shapes j(x) have been derived [4]. These shapes allow us to solve the heat transport problem (4) and (5). Since the temperature variation across the CL is small, this variation does not affect Rreac . This allows us to decouple the thermal and electrochemical (performance) problems. 2.2. The regimes of CL operation The CL performance model used in this work is based on the standard macrohomogeneous model [11,7]. The flooded agglomerate model adds the pre-exponential effectiveness factor in the Butler–Volmer equation for the half-cell reaction rate [2]. This factor takes into account transport loss due to oxygen diffusion
within the agglomerate and for the reasonable set of parameters this factor weakly depends on the local overpotential. Therefore, the results of our work can be used for the catalyst layers with the flooded agglomerates, taking into account the effective lowering of the exchange current density due to oxygen transport in the agglomerates. Note also that the performance model used below is isothermal. This is justified since the temperature variation in the CL is small. The inequalities determining the regime of CL operation have been derived assuming that the rate of electrochemical conversion follows the Butler–Volmer relation [10]. The LC regime takes place when the cell current density j0 is below the critical value jc : j0 < jc
(6)
where
jc = and
2i∗ ion b
ε=
c1 cref
coth
1 ε
ion bcref 2i∗ lt2 c1
(7)
(8)
is the dimensionless Newman’s reaction penetration depth [12]. Here i* is the volumetric exchange current density (A cm−3 ), b is the Tafel slope, c1 is the feed molecules concentration in the CL and cref is their reference concentration. For air-fed cathode, cref is the oxygen molar concentration in air. Typically, the exchange current density is small (e.g., for oxygen reduction reaction) and ε is large. Physically it means, that the Newman’s reaction penetration depth greatly exceeds the CL thickness.4 For ε 1 we have coth (1/ε) ε and Eqs.(6) and (7) simplify to j0 <
ion b lt
(9)
jc =
ion b lt
(10)
Importantly, Eq. (9) does not contain the exchange current density i* . Physically, at large ε (small i* ), the shape of Rreac is controlled by the proton transport. Regardless of the value of i* , the reaction could be either uniformly distributed in the CL (small current, good proton transport), or “concentrated” at the membrane interface (large current, poor proton transport). The critical current density separates these two regimes. The high-current regime arises if the CL proton conductivity is low, or the CL is thick. For the case of Butler–Volmer kinetics, this regime takes place if j0 > 4 ion b/lt [10]. Our calculations show that for the heat transport problem this inequality should be strengthened to j0 >
8ion b lt
(11)
If the cell current density falls into the range ion b 8ion b < j0 < , lt lt
(12)
the CL operates in the transition region between the LC and HC regimes. In this region, the exact solution of performance problem does not exist. Fortunately, in some cases the LC and HC temperature shapes are close to each other, so that the transition-region shape is “trapped” within a narrow domain between the two analytical solutions (see below). For the catalyst layers considered in
4 The inverse relation takes place on the anode side of hydrogen cells, where the exchange current density is very large.
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
9175
Table 1 Parameters for calculations. (am) – assumed, (ct) – calculated. HT-PEMFC Transfer coefficient ˛ Exchange current density i* (A m−3 ) Catalyst layer thickness lt (m) Overpotential 0 (V) Critical current density jc (A m−2 ) Thickness of the conversion domain l* (m) Parameter ˇ Proton conductivity of the electrolyte phase, p (S m−1 ) ORR entropy change, S (J mol−1 K−1 ) Thermal conductivity of CL, (W m−1 K−1 ) Current density j0 (A m−2 ) Liquid water density w (kg m−3 ) Molecular weight of water Mw (kg mol−1 ) Heat capacity of liquid water cpw (J kg−1 K−1 ) Enthalpy of water evaporation Hevap (J mol−1 ) Evaporation rate constant Kevap (atm−1 s−1 ) CL porosity εCL Liquid saturation in the CL s
PEMFC
0.777 [16] 103 [16] 1.5 × 10−4 [16] 0.665 (ct) 320 5 × 10−6 (ct) 2.858 (ct)
1.0 [13] 103 (am) 10−5 (am) 0.427 (ct) 2955 – – 1 [17] 326.36 [18] 0.27 [19] 104 (am) 103 0.018 4190 41.7 × 103 1 [15] 0.3 (am) 0.3 (am)
this work, the values of critical current density (10) are depicted in Table 1.
Setting here x = 0 we obtain the temperature drop T ≡ T0 − T1 across the CL in the LC regime:
2.3. Low-current regime: the temperature shapes
T =
In the LC regime, the proton current linearly decays with x (Fig. 2):
2.4. High-current regime: the temperature shapes
j(x) = j0 1 −
x lt
(13)
With this relation, Eq. (4) transforms to 2
∂ T − = ∂x2
T S 1 nF
+ 0
j
0
+
lt
j02 ion
x 1− lt
−
∂T = ∂x
T S 1 nF
+ 0
j x 0
+
lt
j02 lt
2
3ion
1− 1−
T S 1 nF
+ 0 j0 +
qlow tot =
nF
x lt
3
j02 lt
(15)
(16)
3ion
(17)
is the proton resistive overpotential and lt ion
is the CL proton resistivity. Integrating Eq. (15) from x to lt we find the temperature profile across the CL:
T (x) = T1 +
+
ˇ 2
ˇ
1−
x lt
(20)
j0 2jc + j0
(21)
and jc is given by (7). Note
that Eq. (21) is an approximate solution to equation j0 = jc ˇ tan ˇ/2 , which results from Eq. (20) with x = 0 [4]. Fig. 2 shows that the electrochemical conversion occurs in a narrow conversion domain at the membrane interface. The characteristic width of this domain is [4] ion b j0
(22)
In this domain, the square of proton current density can be approximated by exponent [4] j2 j02 exp −
j0 Rion = 3
Rion =
(19)
2
x
where ion
j l 0 t
In the HC regime, the proton current decays along x according to [4]
l∗ =
+ 0 + ion j0
3ion 2
where
It is advisable to rewrite this equation in the form
T S 1
+ 0 +
(14)
Setting here x = lt we get the total heat flux leaving the CL [6]: qlow tot =
nF
j = jc ˇ tan
Integrating this equation from 0 to x and taking into account the first of Eq. (5), we obtain
T S 1
T S 1 nF
j02 lt2 3ion
+ 0
1−
j l 0 t
2
1−
x
1
lt
4
1−
x2 lt2
1−
x 3 lt
Note that Eq. (23) approximates the square of the exact expression (20) only close to the membrane (Fig. 2). However, this is exactly what is needed for the heat transport problem: outside the conversion domain, the proton current and the reaction rate are small and this region does not contribute to heat production (Fig. 2). Using Eq. (23) in (4) we get an equation for CL temperature in the HC regime: 2
(23)
l∗
−
∂ T T1 S =− + 0 nF ∂x2
∂j ∂x
+
j02 ion
x
exp −
(24)
l∗
Integrating this equation from 0 to x and taking into account (5), we obtain (18)
−
∂T = ∂x
T S 1 nF
+ 0 (j0 − j) +
j02 l∗ ion
x
1 − exp −
l∗
(25)
9176
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
which describes the linear temperature shape across the CL. Though this linear shape does not satisfy the adiabatic boundary condition at x = 0, it is close to the exact result, Eq. (28) (Fig. 3). Setting in (29) x = 0, we get the temperature drop across the CL:
166 Thigh
164 Tlow, no Joule
T
GDL
Temperature / oC
Thigh, no Joule
Tlow
162
160
0
0.2
0.4
0.6
0.8
1
Fig. 3. Temperature distribution across the 150-m HT-PEMFC cathode catalyst layer calculated with the low- and high-current equations. The low-current curves are denoted as Tlow , the high-current as Thigh . Solid Thigh curve: the exact solution (28), long-dashed curve: the approximate linear solution (29). Parameters for all curves are listed in Table 1.
Setting x = lt here we find high
qtot =
T S 1 nF
+ 0 j0 +
j02 l∗
ion
l t
1 − exp −
(26’)
l∗
In the HC regime, lt /l* 1 and the exponent exp (− lt /l* ) is vanishingly small. Neglecting this exponent, we find the heat flux leaving the CL: high
qtot
T S 1 nF
+ 0 j0 +
j02 l∗
(26)
ion
Rearranging terms and taking into account (22), we come to high
qtot =
T S 1 nF
+ 0 + b j0
(27)
We see that in the HC regime, the Joule heating is equivalent to increase in the half-cell overpotential 0 by the Tafel slope b. In low-temperature fuel cells, b is in the order of 50–100 mV [13], whereas 0 is in the order of 300–500 mV. Thus, in the HC regime, the contribution of the Joule heating to the overall heat flux does not exceed 30% (see also below). In high-temperature solid oxide fuel cells (SOFCs), b 150 − 200 mV [14] and the Joule contribution to the heat flux may exceed 50%. Integrating (25) from x to lt , we get the temperature shape across the CL:
T (x) = T1 +
T S 1 nF
− exp −
+ 0
x
2l∗
x
− exp −
j l 0 t +
1−
x 2l∗ + lt lt
j0 lt b x l∗ 1− + lt lt
exp −
lt 2l∗
l t
exp −
l∗
The exponential terms in Eq. (28) contain a small pre-factor l* /lt 1. These exponents “flatten” the shape T(x) near the membrane to provide the adiabatic boundary condition at the membrane interface (Fig. 3). However, since l* lt , this flattening only occurs in a small conversion domain (Fig. 3). Therefore, to a good approximation exponents in Eq. (28) can be ignored and we come to a simple result T (x) T1 +
T S 1 nF
j l 0 t
+ 0 + b
x 1− lt
nF
j l 0 t
+ 0 + b
(30)
Comparing this to Eq. (19) we see that in the HC regime, the temperature drop due to the reaction is twice as large as in the LC regime. This result could have been foreseen. Ignoring for a moment the Joule heating we note, that the total heat flux from the CL does not depend on the distribution of reaction rate in the CL. Indeed, the total amount of heat energy released in the reaction events per unit time is independent of their spatial distribution. Further, in the HC regime, the temperature shape is close to linear, since in most of the CL thickness, Rreac and j(x) are close to zero. We, therefore, come to a simple geometric problem. The LC temperature shape is parabolic (the Joule term is ignored): Tlow (x) = T1 + a2 (1 − (x2 /lt )2 ), and the HC shape is linear: Thigh (x) =
T1 + c 2 1 − lx , where a and c are constants. Both these shapes t must have the same slope at x = lt ; elementary calculation shows that c2 = 2a2 and hence T|high = 2T|low . 2.5. The general expression for the heat flux Substituting (22) into (26’) and rearranging terms, we come to high
qtot =
T S 1 nF
+ 0 + b 1 − exp −
j0 lt ion b
j0
(31)
Eq. (31) is the exact expression for the HC heat flux from the CL. As discussed above, this expression is simplified to Eq. (27), since in the HC regime, the parameter j0 lt /( ion b) is small. Eq. (31) suggests the general expression for the heat flux from the CL, valid in the whole range of cell current density: qtot =
T S 1 nF
+ 0 + b 1 − exp −
j0 lt 3ion b
j0
(32)
Eq. (32) differs from Eq. (31) by a factor of 3 in the power of exponent. If the LC condition Eq. (9) is fulfilled, the exponent in (32) can be expanded in series; retaining the linear term, we come to the LC Eq. (15). In the HC limit of Eq. (11), the exponent in Eq. (32) can be neglected, and we come to Eq. (27). Thus, in the LC and HC limits, Eq. (32) reduces to the respective exact expression and hence this equation interpolates between the limiting solutions. Note that in Nafion-based cells, ion in Eq. (32) depends on the water content of the membrane phase in the CL. Note also that Eq. (32) does not take into account the situations with the gradient of water content through the CL. 2.6. DMFC cathode: heat flux due to fuel crossover
(28)
l∗
T S 1
(29)
A feature of a direct methanol fuel cell is large fuel crossover through the membrane. Permeated methanol is catalytically burned on the cathode side and gives a significant (sometimes dominating) contribution to the heat flux from the cathode catalyst layer (CCL). Note that methanol comes to the CCL at a temperature close to T1 and hence the methanol flux itself does not change the heat balance in the CCL. The main effect of methanol is heat released in the methanol–oxygen catalytic combustion. The respective heat flux is given by equation [10] qcross
T1 Scross = 6F
ˇ 1+ˇ
(jlim − j0 )
(33)
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
Table 2 The heat flux (kW m−2 ) from the catalyst layer calculated with the low-, high- and general relations (the rows LC, HC and General, respectively). The row “Regime” indicates the regime of CL operation. LC and HC stand for the low- and high-current regimes, respectively, T indicates the transition region between the LC and HC. The row “Reaction” shows the heat flux from the half-cell reaction without the Joule term. The row “Joule” shows the Joule heat flux. The last row indicates the Joule heat flux calculated using the conventional low-current approximation, Eq. (17). The parameters for calculation are listed in Table 1.
Low-current
10
High-current Fitting
j0 (A cm−2 ) Real regime LC equation, Eq. (17) HC equation, Eq.(27) General Eq. (32) Reaction Joule LC Joule
0 0
0.2
0.4
0.6
Current density / A
0.8
1
cm-2
Fig. 4. Heat flux from the HT-PEFC cathode calculated with the low-current Eq. (16) (top curve), high-current Eq. (27) and approximate fitting Eq. (32). The latter two curves are indistinguishable. 0 /298 is the entropy change in the Here Scross = HM 0 is the respective methanol–oxygen combustion reaction, HM enthalpy change, and
ˇ=
Dm lb Db lm
is the crossover parameter. Here Dm , Db are the diffusion coefficients of methanol in the membrane and anode backing layer (ABL), respectively, lm and lb are the membrane and ABL thickness, respectively. The limiting current density jlim is given by jlim =
6FDb cM lb
(34)
where cM is methanol molar concentration in the feed channel. In DMFC cathode, the flux (33) should be added to the right side of Eq. (32). Note that at small cell currents, the heat flux due to methanol crossover, Eq. (33) greatly exceeds all the other heat sources in DMFC. In the opposite limit of j0 jlim this flux is small and the heat sources due to ORR and Joule heating should be taken into account. 3. Numerical examples 3.1. HT-PEMFC cathode Parameters for the calculations in this section are listed in Table 1. For the CCL of HT-PEMFC, the working current density greatly exceeds 8jc (Table 1), i.e., the CCL operates in the HC regime. Nonetheless, for comparison, we also display the LC curves calculated with the same set of parameters. These curves show the temperature shapes which follow from the misconception about the regime of CL operation. Fig. 3 shows the LC and HC shapes T(x). The temperature shapes without the Joule term (the last term in Eqs. (18) and (28) are also displayed. The LC curves exhibit quite significant contribution of Joule heating to the temperature shape (Fig. 3). However, the HC curves show that in fact this contribution is small (Fig. 3). Physically, in the HC regime, the proton current rapidly disappears close to the membrane and the contribution of Joule heat is small. If the CL were operated in the LC regime, the proton current would linearly decay with the distance and fraction of the Joule heating in the overall heat flux would be much larger. This effect also is illustrated in Fig. 4, which shows the heat fluxes from the HT-PEFC cathode calculated with the low-, high- and fitting Eqs. (16), (27) and (32), respectively. As can be seen, as the CCL
HT-PEMFC
PEMFC
1.0 HC 15.31 10.79 10.79 10.31 0.48 5.0
1.0 T 7.50 7.47 7.37 7.17 0.20 0.33
operates in the HC regime, the LC Eq. (15) strongly overestimates the Joule heat flux, whereas Eqs. (27) and (32) give correct and practically indistinguishable results. The numerical data are listed in the first column of Table 2. Note that at small current densities, the contribution of Joule heating is small and all the curves in Fig. 4 coincide. Note also that the curves in Fig. 4 are plotted taking into account the dependence 0 (j0 ) (Eq. (74) in [4]). 3.2. PEMFC cathode The curves for PEMFC cathode are shown in Fig. 5. The pair j0 , jc shows that this cathode works in the transition region (jc < j0 < 8jc , Table 1). However, the gap between the LC and HC curves is small and hence both these solutions are close to the real curve. The PEMFC working current density is large (1 A cm−2 , Table 1); however, the temperature drop in the CCL does not exceed 0.2 K. The main reason is that the CCL of PEMFC is only 10 m thick, 10 times smaller, than in the HT-PEMFC (Table 1). In spite of the small temperature drop, the heat flux from the CCL is about 8 kW m−2 , which is comparable to that flux from the CCL of HTPEMFC (11 kW m−2 , Table 2). 4. Discussion Figs. 3 and 5 demonstrate the difference between the temperature profiles when the CL works in the low and high-current regimes. These shapes do not directly affect the electrochemical performance of the CL; however, they determine the heat flux leaving the CL. If the heat contact between the GDL and current collector
70.2
Temperature / oC
Heat flux / kW m-2
HT-PEFC
70.15
Thigh
Thigh, no Joule
Tlow
70.1
GDL
15
5
9177
Tlow, no Joule
70.05
70
0
0.2
0.4
0.6
0.8
1
Fig. 5. (Color online) Same as in Fig. 3 for the PEMFC cathode. Parameters for all curves are listed in Table 1.
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
is poor, part of this heat flux would be reflected and returned to the CL, thereby heating up the CL and accelerating the reaction. It should be emphasized that the total heat flux leaving the CL does not depend on boundary conditions: this flux is determined solely by the rate of heat production in the CL. Eq. (5) only affect the temperature shape across the CL. Note that Eq. (32) does not specify the direction of the heat flux. On the cathode side of low-T cells, there is a heat flux transported with the liquid water produced in the ORR. This heat flux is implicitly taken into account on the right side of Eq. (32). Physically, ORR generates water molecules in liquid state and part of the reaction heat is spent to warm up these molecules to the working temperature T. The heat flux transported with the flow of liquid water produced in the ORR is qw = Mw cpw Tj0 /(2F). With the data from Table 1 we get qw 1.34 kW m−2 . Comparing this to the total heat flux of 7.5 kW m−2 (Table 2), we see that the liquid water transports about 20% of the total heat produced in the CL.5 Note also that the contribution of water flux through the membrane to the CL heat balance is negligible, since the MEA sandwich is almost isothermal. Fig. 4 and Table 2 indicate that the general equation for the heat flux (32) has accuracy better than 2%. Note that in the case of HTPEMFC, the standard LC Eq. (17) overestimates the total heat flux by 50% (Table 2) due to incorrect (10 times higher) value of the Joule flux. Eq. (32) is recommended for use in CFD calculations. If local current over the cell surface is non-uniform, Eq. (32) allows the user not to care about the regime of CL operation. This equation returns the correct value of the heat flux from domains with low, intermediate and high local current densities. The model discussed ignores the voltage loss due to poor feed molecules transport in the CL. If this loss is high, the reaction rate Rreac has maximum at the CL/GDL interface [4]. In that case, the conversion occurs in a narrow conversion domain at the CL/GDL interface; the thickness of this domain depends on the oxygen diffusion coefficient in the CL. Qualitatively, in this regime, the drop of temperature over the CL thickness would be much smaller, since heat would be effectively removed to the GDL.
5. Conclusions Solutions to the CL performance problem are used to derive expressions for the temperature shape and for the heat flux from the catalyst layer operating in the low- and high-current regimes. In the high-current regime, the contribution of the Joule heating to the total heat flux appears to be much smaller, than in the low-current mode. A simple equation for the heat flux from the CL valid in the whole range of current densities is suggested. The analytical results are illustrated with the numerical examples for the catalyst layers of PEMFC and HT-PEMFC.
Acknowledgement This work has been partly supported by the German Federal Ministry of Economics and Technology, Project No. 0327853A.
5 To be precise, the convective term for the transport of liquid water should be added to the left side of Eq. (2). This would lower the conductive heat flux by the value qw and the same value would appear as the convective heat flux. Note that the total heat flux of Eq. (32) would not change, since this flux is determined by the rate of heat production in the CL only.
1.5
Evaporation heat flux / kW m-2
9178
1
0.5
0
20
40
60
80
100
Temperature / oC Fig. 6. Absolute value of the heat flux due to liquid water evaporation in the catalyst layer, Eq. (35). The curve is an upper estimate with the data from Table 1 and pw = 0, which corresponds to zero pressure of water vapor in the CL.
Appendix A. Heat flux due to water evaporation The heat flux due to liquid water evaporation in the CL is obtained if we multiply the molar transfer rate due to evaporation [15] by the respective enthalpy change Hevap . This yields
qevap = −Kevap εCL slt
Hevap w Mw
(psat w − pw )
(35)
where Kevap is the evaporation rate constant, εCL is the CL porosity, s is the CL liquid saturation, Hevap is the enthalpy of water evaporation, w is the liquid water density, Mw is the molecular weight of water, psat w is the pressure of saturated water vapor and pw is the partial pressure of water vapor in the CL. Physically, the product εCL s is the fraction of CL volume occupied by liquid water, Hevap w /Mw is the energy consumed for evaporation of unit volume of water, Kevap (psat w − pw ) is the “driving force” for evaporation, which is proportional to the difference of saturated and available water vapor pressures. The constant Kevap in an efficient manner takes into account the kinetics of water phase transformation in the porous media of the CL. The pressure of saturated water vapor is given by ˜ (T ) psat w = p0 p where p0 = 1 atm and
p˜ (T ) = exp a0 + a1 (T − 273) + a2 (T − 273)2 + a3 (T − 273)3
(36)
(37)
describes the temperature dependence (the coefficients a0 , . . . a3 are given in Table 3). In models which assume that evaporation occurs in the CL, the (negative) heat flux qevap should be added to the right side of Eq. (32). Written in the form of Eq. (35), the flux qevap is determined by the CL liquid saturation, temperature and density of water vapor in the CL and it does not explicitly depend on the cell current density. However, the dependence on j0 is “hidden” in the liquid saturation s, which is a function of current-dependent water transport in the MEA. The absolute value of heat flux due to water evaporation in the CL calculated with the data from Table 1 is depicted in Fig. 6. As can be seen, at typical working temperatures below 70◦ C, this heat flux does not exceed 1 kW m−2 (cf. the data in Table 2). In PEMFC, the GDL is more than ten times thicker than the CL, and the GDL porosity typically is twice that of the CL. Thus, the heat flux due to evaporation in the GDL is at least an order of magnitude larger and we may conclude that water evaporation occurs predominantly in the GDL.
A.A. Kulikovsky, J. McIntyre / Electrochimica Acta 56 (2011) 9172–9179
9179
Table 3 The coefficients in Eq. (37). a0
a1
a2
a3
−2.1794 ln (10)
0.02953 ln (10)
−9.1837 × 10−5 ln (10)
1.4454 × 10−7 ln (10)
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10]
M. Eikerling, J. Electrochem. Soc. 153 (2006) E58. D. Harvey, J.G. Pharoah, K. Karan, J. Power Sources 179 (2008) 209. A. Weber, J. Newman, Chem. Rev. 104 (2004) 4679. A.A. Kulikovsky, Electrochim. Acta 55 (2010) 6391. C.-Y. Wang, Chem. Rev. 104 (2004) 4727. A.A. Kulikovsky, J. Power Sources 162 (2006) 1236. M. Eikerling, A.A. Kornyshev, J. Electroanal. Chem. 453 (1998) 89. A.Z. Weber, J. Newman, J. Electrochem. Soc. 153 (2006) A2205. A.A. Kulikovsky, Electrochem. Commun. 9 (2007) 6. A.A. Kulikovsky, in: Analytical Modelling of Fuel Cells, Elsevier, Amsterdam, 2010.
[11] M.L. Perry, J. Newman, E.J. Cairns, J. Electrochem. Soc. 145 (1998) 5. [12] J. Newman, in: Electrochemical Systems, Prentice Hall, Inc., Englewood Cliffs, NJ 07632, 1991. [13] K.C. Neyerlin, W. Gu, J. Jorne, H. Gasteiger, J. Electrochem. Soc. 153 (2006) A1955. [14] V.A. Danilov, M.O. Tade, Int. J. Hydrogen Energy 34 (2009) 6876. [15] D. Natarajan, T.V. Nguyen, J. Electrochem. Soc. 148 (2001) A1324. [16] A.A. Kulikovsky, H.F. Oetjen, C. Wannek, Fuel Cells 10 (2010) 363, doi:10.1002/fuce.200900135. [17] Z. Xie, T. Navessin, K. Shi, R. Chow, Q. Wang, D. Song, B. Andreaus, M. Eikerling, Z. Liu, S. Holdcroft, J. Electrochem. Soc. 152 (2005) A1171. [18] M.J. Lampinen, M. Fomino, J. Electrochem. Soc. 140 (1993) 3537. [19] M. Khandelwal, M. Mench, J. Power Sources 161 (2006) 1106.