Water transport through a PEM fuel cell: A one-dimensional model with heat transfer effects

Chemical Engineering Science 64 (2009) 2216 -- 2225

Contents lists available at ScienceDirect

Chemical Engineering Science journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / c e s

Water transport through a PEM fuel cell: A one-dimensional model with heat transfer effects D.S. Falcão a , V.B. Oliveira a , C.M. Rangel b , C. Pinho a , A.M.F.R. Pinto a,∗ a b

CEFT, DEQ, Faculdade de Engenharia da Universidade do Porto, Rua Dr. Roberto Frias, 4200-465 Porto, Portugal Instituto Nacional de Engenharia, Tecnologia e Inovação, Paço do Lumiar, 22, 1649-038, Portugal

A R T I C L E

I N F O

Article history: Received 24 April 2008 Received in revised form 12 January 2009 Accepted 19 January 2009 Available online 7 February 2009 Keywords: PEM fuel cell Mathematical modelling Heat transfer Mass transfer Transport processes Net water transfer coefficient

A B S T R A C T

Models play an important role in fuel cell design/development. The most critical problems to overcome in the proton exchange membrane (PEM) fuel cell technology are the water and thermal management. In this work, a steady-state, one-dimensional model accounting for coupled heat and mass transfer in a single PEM fuel cell is presented. Special attention is devoted to the water transport through the membrane which is assumed to be a combined effect of diffusion and electro-osmotic drag. The transport of heat through the gas diffusion layers is assumed to be a conduction-predominated process and heat generation or consumption is considered in the catalyst layers. The analytical solutions for concentration and net water transport coefficient are compared with recent published experimental data. The operating conditions considered are various cathode and anode relative humidity (RH) values at 80 ◦ C and 2 atm. The studied conditions correspond to relatively low values of RH, conditions of special interest, namely, in the automotive applications. Model predictions were successfully compared to experimental and theoretical I–V polarization curves presented by Hung et al. [2007. Operation-relevant modelling of an experimental proton exchange membrane fuel cell. Journal of Power Sources 171, 728–737] and Ju et al. [2005a. A single-phase, non-isothermal model for PEM fuel cells. International Journal of Heat and Mass Transfer 48, 1303–1315]. The developed easy to implement model using low CPU consumption predicts reasonably well the influence of current density and RH on the net water transport coefficient as well as the oxygen, hydrogen and water vapour concentrations at the anode and cathode. The model can provide suitable operating ranges adequate to different applications (namely low humidity operation) for variable MEA structures. © 2009 Elsevier Ltd. All rights reserved.

1. Introduction Fuel cells offer an innovative alternative to current power sources with potential to achieve higher efficiencies with renewable fuels with minimal environmental impact. In particular, the proton exchange membrane (PEM) fuel cells are today in the focus of interest as one of the most promising developments in power generation with a wide range of applications in transportation and in portable electronics (Biyikoglu, 2005; Hirschenhofer et al., 1994). These fuel cells convert the chemical energy of the fuel and oxidant into electrical energy releasing the balance as heat. Although prototypes of fuel cell vehicles and residential fuel cell systems have already been introduced, their cost must be reduced and their efficiencies enhanced. To improve the system

∗ Corresponding author. E-mail addresses: [email protected] (C.M. Rangel), [email protected] (A.M.F.R. Pinto). 0009-2509/$ - see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2009.01.049

performance, design optimization and analysis of fuel cell systems are important. Mathematical modelling and simulation are needed as tools for design optimization of fuel cells, stacks and fuel cells power systems. Several coupled fluid flow, heat and mass transport processes occur in a fuel cell in conjunction with the electrochemical reactions. These processes have a significant impact on two important operational issues: water and thermal management. To achieve optimal fuel cell performance, it is critical to have an adequate water balance to ensure that the membrane remains hydrated for sufficient proton conductivity while cathode flooding and anode dehydration are avoided (Chang et al., 2002). Water content of the membrane is determined by the balance between water production and three water transport processes: electro-osmotic drag of water (EOD), associated with proton migration through the membrane; back diffusion from the cathode; and diffusion of water to/from the oxidant/fuel gas streams. The net transport of water is therefore determined by this balance and may be quantified in M /N + . Understanding terms of the net transport coefficient  = NH H 2O

D.S. Falcão et al. / Chemical Engineering Science 64 (2009) 2216 -- 2225

the water transport in the PEM (Baschuk and Li, 2000; Wang et al., 2001) is a key issue to avoid cathode flooding and membrane dehydration and can also serve as a guide for materials optimization and development of new MEAs. Different models were developed in the last decade to describe several water transport mechanisms through the membrane such as Springer et al. (1991) using a diffusion model, Bernardi and Verbrugge (1992) considering a hydraulic permeation model and Kulikovsky (2004) developing a semi-analytical model 1D + 1D taking into account oxygen and water transport across the cell and deriving an expression for the limiting current density. More recently, two- and three-dimensional models on computational fluid dynamics (CFD) approach have been presented. A three-dimensional model presented by Natarajan and Nguyen (2003) considers the effects of flowrate, inlet humidity and temperature on the liquid water flooding in the cathode. The works of Meng and Wang (2005) and Wang and Wang (2006) also predict water flooding inside a PEMFC and the liquid water effects in the cell performance. Sukkee and Wang (2006) presented a unified water transport CFCD (computational fuel cell dynamics) model incorporating various modes of water transport (diffusion, convection and electro-osmotic drag). The authors applied the model to elucidate the water management in a three-dimensional fuel cell with dry to low humidified inlet gases. Ju et al. (2007) conducted a two-phase model on PEMFC cathode to address the liquid water saturation. Very recently, Le and Zhou (2008) report a general model specially focused on the liquid water management. The developed three-dimensional unsteady model with detailed thermo-electrochemistry, multispecies and two-phase effects with the interface tracking by using the volume-of-fluid (VOF) method was implemented into the CFD software package FLUENT. Low humidity operation of PEM fuel cells is of increasing interest namely in automotive applications. In these conditions the current density distribution will be strongly dependent on the membrane hydration level which is markedly influenced by temperature. Thermal management is therefore required under these operation conditions. The importance of accounting for temperature gradients in fuel cells modelling was demonstrated in the work presented by Djilali and Lu (2002) and Ramousse et al. (2005). Ju et al. (2005a) developed a three-dimensional, non-isothermal model to account for heat generation mechanisms, including irreversible heat due to electrochemical reactions, entropic heat and Joule heating arising from electrolyte ionic resistance. A model coupling two-phase with temperature effects (solid-phase and fluid-phase temperatures inside a two-layer porous cathode of a polymer electrolyte fuel cell are calculated simultaneously) was presented (Hwang et al., 2007). Nowadays, a general trend has been observed to apply the methods of CFD to fuel cell modelling (Um et al., 2000; Zhou and Liu, 2001; Berning et al., 2002; Ju et al., 2005a; Sukkee and Wang, 2006; Le and Zhou, 2008). While the transport phenomena occurring in the fuel cell are inherently three-dimensional, the practical usefulness of CFD models is relatively low (computing times, for example, are still prohibitive). Some authors (Spakovsky and Olsommer, 2002; Cheddie and Munroe, 2005; Oliveira et al., 2007) refer the need of the smart use of the “old fashioned” analytical or semi-analytical models (using applied mathematical techniques and computing power) to obtain practically useful, reduced models. Validation is also essential when analysing the different available models. There are numerous models “validated” against a polarization curve and using often a significant number of fitting parameters. Min et al. (2006) stated that different groups of parameters can result in almost the same global polarization curve. As pointed out by Neyerlin et al. (2006), many advanced fuel cell performance models have been compromised by using inaccurate or incorrectly defined kinetic parameters. In this work, an easy to implement and using low CPU, steadystate and one-dimensional model accounting for coupled heat and

2217

mass transfer in a single PEM fuel cell is presented. A rationally chosen set of parameters is used. The model predictions of the fuel cell performance (polarization curves) are compared with experimental and one-dimensional model predictions data published by Hung et al. (2007). Special attention is in the present work devoted to the cell water management. The model predictions of the net water transfer coefficient are compared with experimental values recently published by Liu et al. (2007). 2. Model development 2.1. Model assumptions In the development of the model, the fuel cell consists of different layers represented in Fig. 1. The cell consists in an aluminum plate (AL), an acetate sheet (ACE), a cupper current collector (Cu) and a flow channel (C), at the anode and cathode sides and a MEA. The MEA includes the backing layers (BL), the catalytic layers (CL) and the membrane (M). The acetate sheet isolates the end plate (aluminum), from the current collector plate. At the anode side, the fuel H2 (wet or humidified) is oxidized liberating electrons and producing protons while air (wet or humidified) is supplied to the cathode. The free electrons flow to the cathode, via an external circuit and there react with the protons and O2 , producing water and heat. The electrode reactions are: Anodic reaction: H2 (g) → 2H+ (aq) + 2e−

(1)

Cathodic reaction: 1/2O2 (g) + 2H+ (aq) + 2e− → H2 O(g)

(2)

The model presented relies on the following assumptions: • mass and heat transport are steady-state and one-dimensional (direction z in Fig. 1); • ideal gas mixture assumed (operation pressures and pressure gradients are usually small); • heat and mass transport through the gas diffusion and catalyst layers assumed to be a diffusion-predominated process (negligible convection effects); • effective Fick models for the mass transport in the diffusion layers and membrane are considered; • the thermal energy model is based on the differential thermal energy conservation equation (Fourier's law); • the thermal conductivity for all the materials is assumed to be constant; • the anodic and cathodic overpotential is constant through the catalyst layers; • heat generation or consumption is considered in the catalyst layers; • the heat generated by the current flow (Joule effects) through each cell layer is ignored (usually negligible when compared to the other terms in the thermal balance namely from low to moderate current densities); • water transport through the membrane assumed to be a combined effect of diffusion and electro-osmotic drag; • membrane proton conductivity is a function of , the number of water molecules per ionic group; • local equilibrium at interfaces is represented by partition functions; • reactions in the catalyst layers are considered as homogeneous; • kinetics of the cathode is described by a Tafel expression; • kinetics of the anode is described by a Tafel expression; • anode and cathode flow channels are treated as a continuous stirred tank reactor (CSTR), so, the composition and temperature

2218

D.S. Falcão et al. / Chemical Engineering Science 64 (2009) 2216 -- 2225 AC

CC

Tin Air

Tin H2

y Section I Section II

q1

q1

q1

T1

T2

T3

T4

AAL

AACE

ACu

T5

AGraph AC

0

Section II Section I

q3

qACL

q5

qCCL

q4

T6

T7

T8

T9

T10

T11

ABL

ACL

M

CCL

CBL

CC

q2

q2

q2

T12

T13

T14

T15

CGraph

CCu

CACE

CAL

Z

zAC/ABL zABL/ACL zACL/M zCCL/M zCBL/CCL zCC/CBL H2 QAC

T16

Air, H2O QCC

Fig. 1. Schematic representation of a PEM fuel cell.

inside the channels are uniform (this assumption is generally taken in one-dimensional models); • anode and cathode streams act as heat transfer fluids removing heat from the cell at the exit temperatures; • two-phase flows effects are not considered (this assumption usually generates deviations from experimental results for high values of current densities due to the formation of water drops at the cathode side).

In Eqs. (3)–(6), the area As is given by

2.2. Mass transfer model

where j stands for hydrogen (H2 ), oxygen (O2 ) or water (H2 O) and eff ,BL is the effective molecular diffusion coefficient through the BL. Dj For the CL the equations are

The values of reactants and water concentration at the CL were obtained considering mass transfer transport in the flow channels (C) and in BL along the z direction for both cathode and anode. The flow channels are treated as a CSTR and the resulting equations are the following: Q AC

As

AC NH = 2

Icell I AC/ABL AC/ABL in in = s (CH − CH2 ) ⇒ CH2 = CH − cell AC 2 2 2F A 2FQ

(3)

CC NO = 2

Q CC Icell I As CC/CBL CC/CBL = s (COin2 − CO2 ) ⇒ CO2 = COin2 − cell CC 4F A 2FQ

(4)

AC NH = 2O

Icell I As Q AC in AC/ABL AC/ABL in − CH2 O ) ⇒ CH2 O = CH − cellAC = s (CH 2O 2O 2F A 2FQ

CC NH = − (1 + ) 2O CC/CBL

⇒ CH2 O

(5)

( + 1)Icell As 2FQ CC

M NH 2O

Icell /2F

where W and L are the channel width and length, respectively. For the BL the general transport equations are the following: BL eff ,BL dC j

NjBL = −Dj

(9)

dz

CL eff ,CL,CL dC j

NjCL = −Dj

(10)

dz

with j representing hydrogen (H2 ), oxygen (O2 ) or water (H2 O). The species concentrations at the C/BL and BL/CL interfaces (for both cell sides) are giving by assuming local equilibrium with partition coefficients K C/BL and K BL/CL , respectively. The boundary conditions for the set of Eqs. (9) and (10) are At

z = zC/BL → CjBL = K C/BL C C/BL

(11)

At

z = zBL/CL → CjCL = K BL/CL C BL/CL

(12)

(6)

where F is the Faraday constant, Icell the cell current density, Cjin the inlet concentration for species j and the net water transport coefficient, , defined as

=

(8)

Using the set of differential equations (9) and (10) together with the boundary conditions (11) and (12), it is possible to determine the ACL , C CCL , C ACL , C CCL : interface concentrations CH O2 H2 O H2 O 2

Icell Q CC in CC/CBL − CH2 O ) = s (CH 2O 2F A

in = CH + 2O

As = nChannels × W × L

ACL CH 2

⎡ Icell ⎣ K ABL/ACL zAC/ABL zABL/ACL K ABL/ACL zABL/ACL = + − eff ,ABL eff ,ACL eff ,ABL 2F DH2 DH2 DH2 ⎤ K AC/ABL K ABL/ACL As z ⎦ − − eff ,ACL Q AC D H2

(7)

in AC/ABL ABL/ACL + CH K K 2

(13)

D.S. Falcão et al. / Chemical Engineering Science 64 (2009) 2216 -- 2225

COCCL 2

⎡ Icell ⎣ K CBL/CCL zCC/CBL zCBL/CCL K CBL/CCL zCBL/CCL = + − eff ,CBL eff ,CCL eff ,CBL 4F DH2 DO2 DO2 ⎤ K CC/CBL K CBL/CCL As z ⎦ − − eff ,CCL Q CC D

For Gore select membranes (the type of membrane used in the simulations presented in this work) the effective water membrane diffusion coefficient is given (Ye and Wang, 2007) by   1 1 eff − D = 0.5 × 10−6 exp 2416 303 T

O2

+ COin2 K CC/CBL K CBL/CCL ⎡ Icell ⎣ K ABL/ACL zAC/ABL zABL/ACL K ABL/ACL zABL/ACL ACL CH2 O =  + − eff ,ABL eff ,ACL eff ,ABL 2F DH2 O DH2 O DH2 O ⎤ K AC/ABL K ABL/ACL As z ⎦ − − eff ,ACL Q AC D

(14)

CCL CH 2O

−

eff ,CBL

DH2 O

−

K CC/CBL K CBL/CCL As Q AC

−

dry

(z ⎤ z

⎦

DH2 O

in + CH K CC/CBL K CBL/CCL 2O

(16)

The hydrogen and oxygen concentrations at the catalyst layers enable the determination of cell overpotentials as will be explained in Section 2.4. The water concentration at the anode and cathode catalyst layers, ACL and C CCL , is a function of the net water transport coefficient,  CH H2 O 2O (which, in turn, is a function of the membrane water content), and therefore  must be determined. The water content in the membrane, , can be related with the water activity, aH2 O , by O'Hayre et al. (2006):

 = 14aH2 O for 0 < aH2 O  1

(17)

 = 12.6 + 1.4aH2 O for 1 < aH2 O  3

(18)

The water activity can be defined as aH2 O =

CH2 O RT Psat

(19)

where Psat is the water saturation pressure at temperature T and CH2 O is calculated by Eq. (15) or (16). The water formation occurs at the cathode catalyst layer, therefore it was assumed that the water activity is less than 1 for the anode and higher than 1 for the cathode. Accordingly, to determine the water content at the interface ACL/membrane and at the interface CCL/membrane, one must use Eqs. (17) and (18), respectively, with the water activity calculated from Eq. (19). However, as the water activity is also a function of the water concentration, an alternative equation is necessary. The water transport in the membrane can be described by considering electro-osmotic drag and back diffusion: M NH = 2nsat drag 2O

CCL/M

dry eff j Icell  D () − 2F 22 Mm  jz

2.3. Heat model The heat model development is based on the assumptions referred in Section 2.1, knowing the temperatures at the external cell walls. The heat generated in the cell is partly lost to the ambient and partly removed by the cathode and anode fuel streams. Heats involved in cathodic and anodic reactions, qACL and qCCL , were calculated (Shan and Choe, 2005) by   −(S)(T8 + T7 )/2 Icell qACL = a + (25) 2F   −(S)(T10 + T9 )/2 Icell (26) qCCL = c + 4F In these equations the first term represents the heat due to the activation overpotentials at the anode/cathode and the second term represents the entropy change in the anodic/cathodic electrochemical reaction. The heat balances for the anode and cathode streams can be written as •

qAC =

m CpH 2

AC (T6 − Tin )

(27)

CC (T11 − Tin )

(28)

AS •

(20) qCC =

is the electro-osmotic drag coefficient, D is the memwhere nsat drag brane water effective diffusivity, dry is the dry membrane density and Mm the membrane equivalent weight. M Substituting NH by Icell /2F (from Eq. (7)) in Eq. (20) and inte2O grating gives ⎡ ⎤ Icell Mm nsat 11 drag ⎦ ∗ ⎣ (z) = sat  + C exp z (21) eff ndrag 22F dry D

(24)



An iterative process is at this point set-up to determine the  parameter. Equalizing Eq. (17) to Eq. (23) and Eq. (18) to Eq. (24) gives the  parameter and the integration constant as a function of membrane water diffusivity. Setting a value for the membrane water diffusivity enables the calculation of the water contents by Eqs. (17) and (18) or (23) and (24). With the water contents value (an average value of anode and cathode water contents) it is possible to calculate the water membrane diffusivity by Eq. (22). At this point the Solver from Excel is used and the process just stops when the difference between the set and calculated value of the water membrane diffusivity is near zero, by changing the set value.

eff

where C ∗ is the integration constant.

(22)



⎤ Icell Mm nsat 11 drag CCL/M ⎦ ∗ ⎣ ) = sat  + C exp z eff ndrag 22F  D dry

eff ,CCL

3

⎡

(15)

⎡ Icell ⎣ K CBL/CCL zCC/CBL zCBL/CCL = − ( + 1) + eff ,CBL eff ,CCL 2F DH2 O DH2 O K CBL/CCL zCBL/CCL

2

× (2.563 − 0.33 + 0.0264 − 0.000671 )

The values of the water content at the ACL/membrane and CCL/membrane interfaces are given by ⎡ ⎤ Icell Mm nsat 11 drag ACL/M ∗ ACL/M ⎦ (z ) = sat  + C exp ⎣ z (23) eff ndrag 22F  D

H2 O

in + CH K AC/ABL K ABL/ACL 2O

2219

m Cp O2 AS

The overall heat transfer equation is the following: ACL

q

+ qCCL = q1 + q2 + qAC + qCC

where q1 =

T1 − T4 rea

rea = rAAL + rAACE + rACu

(29) (30) (31)

and where rea is the thermal resistance of the anode ending plates.

2220

D.S. Falcão et al. / Chemical Engineering Science 64 (2009) 2216 -- 2225

Also q2 =

T13 − T16 rec

activation and ohmic losses were considered. Thus, the cell voltage can be defined by the following equation: (32)

where rec is the thermal resistance of the cathode ending plates. The following equations result from heat balances to the different cell layers (Fig. 1): q1 = q3 − qAC

(33)

q2 = q4 − qCC

(34)

ACL

q3 = q

− q5

(35)

q4 = qCCL + q5

(36)

2.3.1. Anode and cathode flow channels In a single fuel cell, the graphite plates have channels machined on only one surface, contacting with the adjacent diffusion layer (Fig. 1). To establish the heat transfer equations in the graphite plates and channels, two sections were considered, section I where heat is transferred by conduction and section II where heat is transported by conduction and convection such as in a finned surface, exchanging heat with the channel fluid (Çengel, 1998). The resulting equations are q1 = −kAC

dT dz

(37)

q2 = −kCC

dT dz

(38)

where kACgraph and kCCgraph are the anode and cathode graphite plate thermal conductivities. 2.3.2. Anode and cathode CL The temperature profiles in the anode and cathode CL are obtained from the second Fourier's law and the heat fluxes by the first Fourier's law: d2 T dz

2

d2 T dz

2

= =

qACL ACL

kACL 

qCCL CCL kCCL 

(39)

(40)

q3 = −kACL

jT jz

(41)

q5 = −kCCL

jT jz

(42)

The boundary conditions for the integration of Eqs. (39) and (40) are the temperatures at both sides of the catalyst layers. The heat generated by the electrode reactions (qACL and qCCL ) depends on the activation overpotentials which in turn depend on temperatures in the CL, therefore heat and mass transfer models must be solved simultaneously. The numerical tool used to perform this task is Matlab. 2.4. Cell performance The characteristic shape of the voltage/current density results from three major irreversibilities: activation losses (caused by the slowness of the reactions taking place on the surface of the electrodes), ohmic losses (losses due to the materials resistances to the current transport) and mass transport or concentration losses (losses due to the reduction of reactants concentration). In this model only

Vcell = E0 − a − c − ohm

(43)

Anode overpotential and cathode overpotential are defined by ⎛  ⎞ ref Icell CH2 ⎠ RT (44) a = ln ⎝ ref CCL 2 a F i 0 CH 2 ⎛

c = ln ⎝

ref

Icell CO2



 ref i0 COCCL 2

⎞ ⎠ RT 4c F

(45)

The ohmic overpotential is defined by

ohm = Icell Rm

(46)

The membrane resistance, Rm , must be known to calculate the ohmic overpotential. Since the membrane conductivity can change locally (depending on the amount of water in the membrane) the total resistance was defined as: RM =

m 0

jz [(z)]

(47)

The conductivity and, consequently, the membrane resistance are a function of the amount of water in the membrane. Conductivity of the Gore select membranes (type of membranes used in the model predictions) can be determined (Ju et al., 2005b) from  

1 1 − (48) (T, ) = 0.5 × (0.005193 − 0.00326) exp 1268 303 T The ohmic overpotential can be calculated after integration of Eq. (47). The results shown in the next section were obtained using the set of parameters presented in Tables 1 and 2. All parameters were carefully chosen from recent literature, namely reference exchange current density and transfer coefficients. In numerous published works, there are numerous values for the same parameter, and it seems that some authors use the parameters that better fit their experimental results and not the most adequate to their operating/design conditions. Namely, the exchange current density is linked especially to the cell temperature, so, this value should evaluated at the cell operation temperature. Exchange current densities (see notes attached to Table 1) and transfer coefficients (c as 0.2 and a as 0.5) were chosen based on recommendations from two recent specific works for anode and cathode kinetics. More details can be found in Neyerlin et al. (2006, 2007). 3. Results and discussion The numerical tools used to implement the model were Matlab and Excel. In order to validate the model predictions in terms of polarization curves, the model results are compared with data from Hung et al. (2007). These authors developed a one-dimensional water and thermal management model considering the effects of water transport across the membrane, activation, ohmic and concentration overpotentials, pressure drops, and current density distribution along the channel of a PEM fuel cell. In their work, the authors obtained design and modelling parameters via regression from four sets of experimental data. In the present work, as explained in the previous section, all parameters are carefully chosen according to the operating conditions, following the suggestions of literature (Neyerlin et al., 2006, 2007). The results of the model developed by Hung et al.

D.S. Falcão et al. / Chemical Engineering Science 64 (2009) 2216 -- 2225

2221

Table 1 Kinetic parameter values used for the presented polarization curves simulations. Parameter

Value

Reference ⎡⎛

1.23 − 0.9 × 10−3 (T − 298) +

E0

⎞2 ⎛

PO 2.303RT ⎢ PH log ⎣⎝ ref2 ⎠ ⎝ ref2 4F PH PO 2

ref

ref

i0,H i0,O Liu et al. ref

2

2 ref

i0,H , i0,O Hung et al. 2 2 Ea 0.9V Ec kroa

a c a c

⎞⎤ ⎠⎥ ⎦V

Neyerlin et al. (2006)

2

0.483 A/cm2 a , 4.28 × 10−5 A/cm2 b

Neyerlin et al. (2007)/Neyerlin et al. (2006)

0.242 A/cm2 a , 5.25 × 10−5 A/cm2 b 10 kJ/mol 10 kJ/mol 1.0 0.5 0.2 0.5 1.0

Neyerlin et al. (2007)/Neyerlin et al. (2006) Jiang and Kucernak (2004) Neyerlin et al. (2006) Spiegel (2008) Neyerlin et al. (2007) Neyerlin et al. (2006) Min et al. (2006) Min et al. (2006)

a Calculated for the temperature of 333 K and H2 pressure of 2 atm (Liu et al.) and for the temperature of 343 K and H2 pressure of 1 atm (Hung et al.), based in the following equation, adapted for anode kinetics from Neyerlin et al. (2006):

⎛ i0 (T) = i0 (Tr )⎝ ref

ref

PH2 ref PH 2

⎞kroa ⎠

 exp

Ea

R

1 Tref

 −

 1 T

Value for 353 K and PHydrogen ≈ 1 atm = 0.267 A/cm2 . b Calculated for the temperature of 333 K and O2 pressure of 0.42 atm (Liu et al.) and for the temperature of 343 K and O2 pressure of 0.21 atm (Hung et al.), following equation proposed by Neyerlin et al. (2006). ⎛ i0 (T) = i0 (Tr )⎝ ref

ref

⎞0.79 ⎛ PO2 ref

PO

⎠

⎝

2

PH2 ref

PH

⎞nc /2 ⎠

 exp

2

−Ec0.9V R

1 Tref

 −

 1 T

Value for 353 K and Pair ≈ 1 atm = 1.8 × 10−4 A/cm2 .

Table 2 Parameter values used in simulations. Parameter M

 Hung et al. (2007) M Liu et al. (2007) ABL , CBL Hung et al. (2007) ABL , CBL Liu et al. (2007) ACL , CCL AC , CC ABL , CBL ACL , CCL

Pa · Pc Hung et al. (2007) Pa · Pc Liu et al. (2007) DACL,CCL H 2

DCBL,CCL O 2

DHABL,ACL 2O

Value

Reference

0.0129 cm 0.003 cm 5.55 × 10−6 cm 0.006 cm 0.001 cm 0.05 cm 0.5 0.6 1 atm 2 atm 0.5 10−3 × T 1.75 × (0.5555) ACL,CCL2.5

Hung et al. (2007) Liu et al. (2007) Hung et al. (2007) Ju et al. (2005b) Ju et al. (2005b) Liu et al. (2007) Ju et al. (2005b) Ju et al. (2005b) Hung et al. (2007) Liu et al. (2007)



CBL,CCL

cm2 /s

Reid et al. (1977)

cm2 /s 2 P × (5.17) 0.5 1.75 ×T × (0.5555) cm2 /s 2 P × (4.25) 0.5 1.75 ×T × (0.1013) cm2 /s 2 P × (5.05)

Reid et al. (1977)

2

2.5

ABL,ACL2.5



2.5

10

−3

10

−3

−3

10

P × (4.25) 0.5 × T 1.75 × (0.067)

DHABL,CBL O

ABL,CBL

nsat drag

1.07 0.002 kg/cm3 0.95 kg/mol 1 cm 0.03 cm 0.06 cm 0.35 cm 0.0095 W/cm K 0.015 W/cm K (1 − ACL ) × (0.66 × 71 + 0.34 × 117) × 10−2 + ACL kH2 W/m K (1 − CCL ) × 71 × 10−2 + CCL kair W/m K 343.15 K 353.15 K 333.15 K 353.15 K

2

dry

Mm

AAL , CAL AACE , CACE ACu , CCu AGraph · CGraph kM kABL · kCBL kACL kCCL TinAC · TinCC Hung et al. (2007) T1 · T16 Hung et al. (2007) TinAC · TinCC Liu et al. (2007) T1 · T16 Liu et al. (2007)

(2007) were validated with experimental data for several operating conditions (e.g., fuel cell temperature, anode and cathode pressure, hydrogen and air stoichiometric ratio and hydrogen and air humid-

Reid et al. (1977) Reid et al. (1977) Ye and Wang (2007) Ju et al. (2005b) Liu et al. (2007) Assumed Assumed Assumed Assumed Ju et al. (2005b) Ju et al. (2005b) Schultz (2004) Schultz (2004) Hung et al. (2007) Hung et al. (2007) Liu et al. (2007) Liu et al. (2007)

ification temperature). The cell used had an active area of 25 cm2 with a Gore select membrane and the channels configuration was dual serpentine. More details can be found in the referred work.

2222

D.S. Falcão et al. / Chemical Engineering Science 64 (2009) 2216 -- 2225

1.2

0.0 0.0

0.1

0.2

0.3 0.4 I (A/cm2)

0.5

0.6

0.7

Fig. 2. Comparison between the present model predictions and experimental data and model predictions of Hung et al. (2007) for a cell temperature of 353 K with H2 /air flowrate = 1.5 × /2.5×, H2 /air humidification temperatures = 343 K/343 K and H2 /air pressures = 1 atm/1 atm.

In Fig. 2 the experimental and modelling results from Hung et al. (2007) are plotted together with the predictions from the present model, for an experiment with a cell temperature of 80 ◦ C, gas pressures of 1 atm and humidification temperatures of 70 ◦ C. As can be seen in Fig. 2, the developed model compares well with both experimental and modelling results from the selected work). It must be emphasized that the set of parameters used in the present model predictions was chosen according to operating conditions following recommendations from literature while the authors used a regression based in experimental data to obtain the best fitted model parameters. Due to its importance in PEM fuel cells development, particular attention is devoted to the water management in the cell. Most of the published research on the water quantification in the cell has been the estimation of the net water crossover through the determination of average values through water collection at the anode and cathode exhausts and overall balances. Liu et al. (2007) have recently reported for the first time experimental measurements of the net water coefficient distribution through a Gore select membrane (30 m thick and active area of 6 cm2 ; catalyst load of 0.4 mg/cm2 on the anode and cathode sides) along the flow direction by simultaneously measuring local current density and water concentration. The cell is segmented in 10 sub cells. By measuring the local current at the 10 segments simultaneously, the current distribution in the operating cell was mapped. A microelectric actuator was used in the anode and cathode sides in order to measure the species concentrations. The actuator allows a connection to a microgas chromatograph. More details about the experimental set up can be found in Lu et al. (2007). As the developed model is one-dimensional along z it enables the determination of an unique value of the net water transport coefficient, to compare with the  values reported by Liu et al. (2007), an average of the values presented by the authors for different operated conditions was used in the model validation. The operating conditions considered (the same as those reported) by Liu et al. (2007) are various cathode and anode relative humidity (RH) values at 80 ◦ C and 2 atm. The studied conditions correspond to relatively low values of RH, conditions of special interest, namely, in the automotive applications. The model predictions for the water concentration profile along the cell are shown in Fig. 3. For the same RH=42% on both sides, and a high current density of 0.685 A/cm2 , the water concentration in the membrane is higher at the cathode side because, under these operating conditions, the amount of water transferred by back diffusion is less than the amount resulted from the sum of electro-osmotic drag and water formation. In Fig. 4, mean values of the water content in the membrane (determined at the middle of the membrane) are represented against

CH2O (mol/cm3)

0.2

Model Exp Hung et al. (2007) Mod Hung et al. (2007)

ACL M

CCL

CBL

CC

1.00E-02 1.00E-03 1.00E-04 1.00E-05 1.00E-06 1.00E-07 1.00E-08 0.048

0.053

0.058 z (cm)

0.063

0.068

Fig. 3. Water concentration profile along the cell different layers for RHa = 42% and RHc = 42% at 0.685 A/cm2 and for a flow H2 /air fuel cell with a = 2 at 1 A/cm2 and

c = 2 at 1 A/cm2 .

13.0 12.0 Water Content λ

E (V)

0.4

ABL

1.00E-01

0.8 0.6

AC

1.00E+00

1.0

11.0 10.0 9.0

RHa = 42%, RHc = 42% RHa = 0%, RHc = 42%

8.0 7.0 0.0

RHa = 42%, RHc = 0% 0.2

0.4

0.6

0.8

1.0

1.2

I (A/cm2) Fig. 4. Average membrane water content variation with current density for different operating conditions and for a flow H2 /air fuel cell with a = 2 at 1 A/cm2 and c = 2 at 1 A/cm2 .

the current density for different reactants humidity levels. Higher values of water content were obtained when both reactants are humidified at a level of 42%. For the other two conditions, the most favourable values were obtained when only the cathode stream was humidified. Although no water enters in the anode side, the water transport by back diffusion is enough to maintain the membrane reasonably hydrated. Also evident from Fig. 4 is that the water content decreases with the increasing current density. For high current densities, the water transport from the anode to the cathode by electro-osmotic drag exceeds the water transport by back diffusion from the cathode to the anode and the membrane will eventually dry out. The water back diffusion is limited by the membrane pores shrink due to dehydration. Therefore, the water transport by back diffusion is not enough to prevent membrane dehydration. Humidified anode and cathode streams are needed to avoid the membrane dehydration, particularly at high current densities. Predicted values of the water content across the membrane are shown in Fig. 5. As expected, the amount of water is higher for the cathode side due to the relevance of water transport by electroosmotic drag and water generation by reaction. As can also be seen from the plots, the water content near the anode catalyst layer is lower for the three curves, in particular for the dry anode. Model predictions of the net water transport coefficient are represented in Fig. 6 as a function of the current density. The values of

D.S. Falcão et al. / Chemical Engineering Science 64 (2009) 2216 -- 2225

Table 4 Hydrogen, oxygen and water molar fractions and net water transport coefficients from Liu et al. (2007), model predictions and absolute deviation for RHa = 42%, RHc = 0% at 0.35 A/cm2 .

14.0

Water content λ

2223

12.0

10.0 RHa = 42%, RHc = 42% RHa = 0%, RHc = 42%

8.0

Parameter

Data from Liu et al. (2007)

Model

Absolute deviation

xH2 xO2 xH2 Oa xH2 Oc

0.93 0.19 0.08 0.05 0.16

0.92 0.19 0.05 0.02 −0.22

0.01 0.00 0.03 0.03 0.38



RHa = 42%, RHc = 0% 6.0 0.000

0.001

0.002

0.003

z (cm) Fig. 5. Water content along the membrane for different operating conditions and for a flow H2 /air fuel cell with a = 2 at 1 A/cm2 and c = 2 at 1 A/cm2 .

2.0

Table 5 Hydrogen, oxygen and water molar fractions and net water transport coefficients from Liu et al. (2007) model predictions and absolute deviation for RHa =0%, RHc =42% at 0.35 A/cm2 . Parameter

Data from Liu et al. (2007)

Model

Absolute deviation

xH2 xO2 xH2 Oa xH2 Oc

0.98 0.19 0.02 0.09 −0.15

0.95 0.20 0.07 0.03 −0.61

0.03 0.01 0.05 0.06 0.47



1.0 0.0

-2.0

RHa = 42%, RHc = 42%

357

RHa = 0%, RHc = 42%

356 T (K)

α

358 -1.0

RHa = 42%, RHc = 0% -3.0

TCCL

I=0.35A/cm2 I=0.685 A/cm2

355 354

-4.0 0.0

0.2

0.4

0.6 I (A/cm2)

0.8

1.0

1.2

352 0.0

Fig. 6. Net water transport coefficient value variation with current density for different operating conditions and for a flow H2 /air fuel cell with a = 2 at 1 A/cm2 and c = 2 at 1 A/cm2 .

Table 3 Hydrogen, oxygen and water molar fractions and net water transport coefficients from Liu et al. (2007), model predictions and absolute deviation for RHa = 42%, RHc = 42% at 0.685 A/cm2 . Parameter

Data from Liu et al. (2007)

Model

Absolute deviation

xH2 xO2 xH2 Oa xH2 Oc

0.90 0.17 0.10 0.14 0.03

0.85 0.18 0.09 0.08 0.64

0.05 0.01 0.01 0.06 0.56



353

 increase with increasing current densities, again due to the importance of electro-osmotic drag when compared to water back diffusion. For the studied conditions the values of  become near positive or positive for very high current densities. Negative values of  mean that the net water flux in the membrane is toward the anode, working under negative values of  helps guaranteeing the membrane humidification. Plots like those in Fig. 6 together with information obtained from Figs. 4 and 5 give very useful information to select operating conditions avoiding membrane dehydration. For validation, the model predictions are compared with the results presented by Liu et al. (2007). The values of hydrogen, oxygen and water molar fractions at the gas diffusion layer interface and of the net water transport coefficient are presented in Tables 3–5, together with the absolute deviation between model and experimental

0.5

1.0

1.5 2.0 z (cm)

2.5

3.0

3.5

Fig. 7. Temperature profile along the cell for RHa = 42% and RHc = 42% at two different current densities and for a flow H2 /air fuel cell with a = 2 at 1 A/cm2 and

c = 2 at 1 A/cm2 .

results. For all the studied conditions an average value for the experimental  was extracted from the published data (Liu et al., 2007) along the y direction. The model predictions for the hydrogen, oxygen and water molar fractions show very good agreement with experimental data (low absolute deviations). Relatively to the net water transport coefficient, the absolute deviations are higher. These differences are explained mainly by the fact that a mean value was used from the variable data along the flow channel, published by the authors. The original data show significant variations along the channel. An example of predicted temperature profiles along the cell is presented in Fig. 7. As expected, for identical RH values, higher temperatures are obtained for the higher current density, corresponding to a higher heat of reaction. The maximum temperature reached for both situations occurs at the cathode CL consistently with the exothermal reaction occurring here. In Fig. 8, the polarization curves, predicted by the developed heat and mass transfer model, are presented. The curve is more favourable when both gas streams are humidified due to the higher membrane hydration. For the remaining two situations, there is no significant difference between curves for these operating conditions. The differences between the polarization curves are more evident for the higher values of current density, since the electro-osmotic drag also

2224

D.S. Falcão et al. / Chemical Engineering Science 64 (2009) 2216 -- 2225

1.2

K L

partition coefficient length, cm

1.0

m Mm nchannels nsat drag N P q Q r R RH Rm T V W z

mass flow, mol/s membrane equivalent weight, kg/mol number of channels electro-osmotic drag coefficient molar flux, mol/cm2 s pressure, atm thermal flux, W/cm2 volumic flow, cm3 /s thermal resistance, cm2 K/W universal gas constant, J/mol K relative humidity, % membrane electrical resistance,  cm2 temperature, K voltage, V width, cm coordinate, cm

•

E (V)

0.8 0.6 Rha = 42%, RHc = 42%

0.4

Rha = 0%, RHc = 42% 0.2

Rha = 42%, RHc = 0%

0.0 0.0

0.2

0.4

0.6 I (A/cm2)

0.8

1.0

1.2

Fig. 8. Predicted polarization curves for different operating conditions and for a flow H2 /air fuel cell with a = 2 at 1 A/cm2 and c = 2 at 1 A/cm2 .

increases. For these conditions the anode stream humidification becomes more relevant. 4. Conclusions A simple, using low CPU, steady-state, one-dimensional model accounting for coupled heat and mass transfer occurring in a PEM fuel cell is presented. The model outputs are the temperature and concentration across the cell and the water content in the membrane. Model predictions were successfully compared to experimental and theoretical I–V polarization curves presented by Hung et al. (2007) and Ju et al. (2005a). The set of parameters used in the present model predictions was rationally chosen according to operating conditions following recommendations from literature. The model predicts reasonably well the influence of current density and RH on the net water transport coefficient. Humidified cathode and specially humidified anode streams are needed to avoid the membrane dehydration, particularly at high current densities. The identification of working conditions corresponding to low and even negative values of the net water transport coefficient helps avoiding the membrane dehydration. The model predictions of the oxygen, hydrogen and water concentrations at the anode and cathode show very good agreement with recent experimental data. The model can provide suitable operating ranges adequate to different applications (namely low humidity operation) for variable MEA structures. This easily to implement reduced model is suitable for use in real-time PEMFC simulations. Notation a As C C∗ Cp Deff eff

D E0 F i0 Icell k kro

activity channels area, cm 2 concentration, mol/cm3 integration constant heat capacity, J/mol K effective diffusivity, cm2 /s water membrane effective diffusivity, cm2 /s open circuit voltage, V Faraday's constant, C/mol exchange current density, A/cm2 current density, A/cm2 thermal conductivity, W/cm K kinetic reaction order

Greek letters

 a c a c  E S   dry  

net water transfer coefficient anode transfer coefficient cathode transfer coefficient anode concentration dependence cathode concentration dependence thickness, cm activation energy, kJ/mol entropy, J/kmol overpotential, V membrane water content dry membrane density, kg/cm3 conductivity, S/cm stoichiometry flow ratio

Subscripts A C ea ec H2 H2 O in j O2 Ohm Sat

anode cathode ending anode plates ending cathode plates hydrogen water inlet species j oxygen ohmic saturation value

Superscripts AACE AAL ABL AC ACL ACu AGraph CACE CAL CBL CC CCL CCu CGraph in M Ref

acetate anode sheet anode aluminium plate anode backing layer anode channel anode catalyst layer anode copper plate anode graphite plate acetate cathode sheet cathode aluminium plate cathode backing layer cathode channel cathode catalyst layer cathode copper plate cathode graphite plate inlet membrane reference value

D.S. Falcão et al. / Chemical Engineering Science 64 (2009) 2216 -- 2225

Acknowledgement The partial support of “Fundação para a Ciência e Tecnologia— Portugal” through Project POCI/EME/55497/2004 and scholarship SFRH/BD/28166/2006 is gratefully acknowledged. POCI (FEDER) also supported this work via CEFT. References Baschuk, J.J., Li, X., 2000. Modeling of polymer electrolyte membrane fuel cells with variable degrees of water flooding. Journal of Power Sources 86, 181–195. Bernardi, D.M., Verbrugge, M.W., 1992. A mathematical model of the solid-polymerelectrolyte fuel cell. Journal of Electrochemical Society 139 (9), 2477–2491. Berning, T., Lu, D.M., Djilali, N., 2002. Three-dimensional computational analysis of transport phenomena in a PEM fuel cell. Journal of Power Sources 106, 284–294. Biyikoglu, A., 2005. Review of proton exchange fuel cell models. International Journal of Hydrogen Energy 30, 1185–1212. Çengel, Y.A., 1998. Heat Transfer a Practical Approach. McGraw-Hill, New York. Chang, H., Kim, J.R., Cho, S.Y., Kim, H.K., Choi, K.H., 2002. Materials and processes for small fuel cells. Solid State Ionics 8312. Cheddie, D., Munroe, N., 2005. Review and comparison of approaches to proton exchange membrane fuel cell modelling. Journal of Power Sources 147, 72–84. Djilali, N., Lu, M., 2002. Influence of heat transfer on gas and water transport in fuel cells. International Journal of Thermal Science 41, 29–40. Hirschenhofer, J.H., Stauffer, D.B., Englemen, R.R., 1994. Fuel Cells Handbook. U.S. Department of Fossils Energy, Morgantown Energy Technology Center, WV, USA. Hung, A.J., Sung, L.Y., Chen, Y.H., Yu, C.C., 2007. Operation-relevant modelling of an experimental proton exchange membrane fuel cell. Journal of Power Sources 171, 728–737. Hwang, J.J., Chao, C.H., Chang, C.L., Ho, W.Y., Wang, D.Y., 2007. Modeling of twophase temperatures in a two-layer porous cathode of polymer electrolyte fuel cells. International Journal of Hydrogen Energy 32 (3), 405–414. Jiang, J., Kucernak, A., 2004. Investigations of fuel cell reactions at the composite microelectrode/solid polymer electrolyte interface. I. Hydrogen oxidation at the nanostructured Pt/Nafion membrane interface. Journal of Electroanalytical Chemistry 567, 123–137. Ju, H., Meng, H., Wang, C.Y., 2005a. A single-phase, non-isothermal model for PEM fuel cells. International Journal of Heat and Mass Transfer 48, 1303–1315. Ju, H., Wang, C.Y., Cleghorn, S., Beusher, U., 2005b. Nonisothermal modeling of polymer electrolyte fuel cells I. Experimental validation. Journal of Electrochemical Society 152 (8), A1645–A1653. Ju, H., Luo, C., Wang, C.Y., 2007. Journal of Electrochemical Society 154, B218. Kulikovsky, A.A., 2004. Semi-analytical 1D + 1D model of a PEMFC. Electrochemistry Communications 6, 969–977. Le, A.D., Zhou, B., 2008. A general model of proton exchange membrane fuel cell. Journal of Power Sources 182, 197–222. Liu, F., Lu, G., Wang, C.Y., 2007. Water transport coefficient distribution through the membrane in a polymer electrolyte fuel cell. Journal of Membrane Science 287, 126–131.

2225

Lu, G.Q., Liu, F.Q., Wang, C.Y., 2007. An approach to measuring spatially resolved crossover coefficient in a polymer electrolyte fuel cell. Journal of Power Sources 164, 134–140. Meng, H., Wang, C.Y., 2005. Journal of Electrochemical Society 152, A1733. Min, C.H., He, Y.L., Liu, X.L., Yin, B.H., Jiang, W., Tao, W.Q., 2006. Parameter sensitivity and discussion of PEM fuel cell simulation model validation Part II: results of sensitivity analysis and validation of the model. Journal of Power Sources 160, 374–385. Natarajan, D., Nguyen, T.V., 2003. Three-dimensional effects of liquid water flooding in the cathode of a PEM fuel cell. Journal of Power Sources 115, 66–80. Neyerlin, K.C., Gu, W., Jorne, J., Gasteiger, H.A., 2006. Determination of catalyst unique parameters for the oxygen reduction reaction in a PEMFC. Journal of the Electrochemical Society 153 (10), A1955–A1963. Neyerlin, K.C., Gu, W., Jorne, J., Gasteiger, H.A., 2007. Study of the exchange current density for the hydrogen oxidation and evolution reactions. Journal of the Electrochemical Society 154 (7), B631–B635. O'Hayre, R., Cha, S.W., Colella, W., Prinz, F.B., 2006. Fuel Cells Fundamentals. Wiley, New York. Oliveira, S.B., Falcão, D.S., Rangel, C.M., Pinto, A.M.F.R., 2007. A comparative study of approaches to direct methanol fuel cell modelling. International Journal of Hydrogen Energy 32, 415–424. Ramousse, J., Deseure, J., Lottin, O., Didierjean, S., Maillet, D., 2005. Modelling of heat, mass and charge transfer in a PEMFC single cell. Journal of Power Sources 145, 416–427. Reid, R.C., Prausnitz, J.M., Sherwood, T.K., 1977. The Properties of Gases and Liquids. McGraw-Hill, New York. Schultz, T., 2004. Experimental and model-based analysis of the steady-state and dynamic operating behaviour of the direct methanol fuel cell (DMFC). Ph.D. Thesis. Shan, Y., Choe, S.Y., 2005. A high dynamic PEM fuel cell model with temperature effects. Journal of Power Sources 145, 30–39. Spakovsky, M.R., Olsommer, B., 2002. Fuel cell systems and system modeling and analysis perspectives for fuel cell development. Energy Conversion and Management 43, 1249–1257. Spiegel, C., 2008. PEM Fuel Cell Modeling and Simulation Using MATLAB. Elsevier, Amsterdam. p. 51. Springer, T.E., Zawodzinski, T.A., Gottesfeld, S., 1991. Polymer-electrolyte fuel cell model. Journal of Electrochemical Society. 138 (8), 2334–2342. Sukkee, U., Wang, C.Y., 2006. Computational study of water transport in proton exchange membrane fuel cells. Journal of Power Sources 156, 211–223. Um, S., Wang, C.Y., Chen, K.S., 2000. Computational fluid dynamics modeling of proton exchange membrane fuel cells. Journal of Electrochemical Society 147 (12), 4485–4493. Wang, Y., Wang, C.Y., 2006. Journal of Electrochemical Society 153, A1193. Wang, Z.H., Wang, C.Y., Chen, K.S., 2001. Two-phase flow and transport in the air cathode of proton exchange membrane fuel cells. Journal of Power Sources 94, 40–50. Ye, X., Wang, C.Y., 2007. Measurement of water transport through membraneelectrode assemblies I. Membranes. Journal of Electrochemical Society 154 (7), B676–B682. Zhou, T., Liu, H., 2001. A general three-dimensional model for proton exchange membrane fuel cells. International Journal of Transport Phenomena 3 (3), 177–198.