Effects of temperature and humidification levels on the performance of a proton exchange membrane fuel cell

http://pia.sagepub.com/

Energy

Engineers, Part A: Journal of Power and

http://pia.sagepub.com/content/220/5/435

The online version of this article can be found at:

DOI: 10.1243/09576509JPE180

2006 220: 435

Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy

M-S Chiang and H-S Chu

Fuel Cell

Effects of Temperature and Humidification Levels on the Performance of a Proton Exchange Membrane

Published by:

http://www.sagepublications.com

On behalf of:

Institution of Mechanical Engineers

be found at:

can

Proceedings of the Institution of Mechanical Engineers, Part A: Journal of Power and Energy

Additional services and information for

http://pia.sagepub.com/cgi/alerts Email Alerts: http://pia.sagepub.com/subscriptions Subscriptions: http://www.sagepub.com/journalsReprints.nav Reprints: http://www.sagepub.com/journalsPermissions.nav Permissions: http://pia.sagepub.com/content/220/5/435.refs.html Citations:

What is This?

- Aug 1, 2006

Version of Record

>>

Effects of temperature and humidification levels on the

performance of a proton exchange membrane fuel cell

M-S Chiang1,2
and H-S Chu1

1_{Department of Mechanical Engineering, National Chiao Tung University, Hsinchu, Taiwan, Republic of China} 2

Department of Mechanical Engineering, Nan Kai Institute of Technology, Nantou, Taiwan, Republic of China The manuscript was received on 4 April 2005 and was accepted after revision for publication on 7 March 2006. DOI: 10.1243/09576509JPE180

Abstract: Numerical investigation of a proton exchange membrane fuel cell performance sub-jected to various humidification and thermal conditions is the focus of this study. Governing equations describing species, mass, momentum, and enthalpy conservation are employed and solved by a computational fluid dynamic algorithm to obtain domain’s physical properties and cell performance. The model accounts for electrochemical kinetic as well as two-phase flow with phase change and water transport. Numerical prediction results use polarization curves and contour plots to illustrate the effects of various humidification schemes and temperature gradient scenarios. Findings show that humidification-level perturbation on the anode or cath-ode side creates different effects. Mechanisms influencing performance-variation tendencies are interpreted. In addition modelling results with existing temperature gradient exhibit differ-ent trends on the overpotdiffer-entials depending on the slope and magnitude. At higher cathode temperature, it is shown that polarization curves are dominated by reaction kinetics and membrane water content at medium and high reaction rates, respectively.

Keywords: proton exchange membrane fuel cell, water and thermal management, humidifi-cation level, temperature gradient

1 INTRODUCTION

Fuel cell systems convert fuel chemical energy directly into electrical energy. When a traditional thermal engine combustion process is replaced by an electro-chemical reaction, it has advantages of quietness, high efficiency, low emission, and simplicity. Given these characteristics, the fuel cell is a revolutionary 21st century power conversion device with appli-cations to vehicles, electronic equipments, and other stationary power generation systems. The proton exchange membrane fuel cell (PEMFC), using solid polymers as working ions’ passages, has recently received much attention and intensive study as a promising candidate of these applications [1].

During the fuel cell power conversion process, chemical energy stored in the bounded atom

electronic orbits is liberated and converged into elec-trical energy through electrocatalytic activity of cata-lysts in the electrodes. Catacata-lysts macroscopically liberate electrons of hydrogen molecules which are converted into protons simultaneously at the anode. These two charged species are transported through individual gateways by electric field func-tions and other mechanisms. The electrons move along the external circuit and release electrical energy through potential drop on the loading. Proton transport is driven by a set of more compli-cated factors. The proton exchange membrane (PEM) must be hydrated for ion movement facili-tation. It uses polytetrafluoroethylene as the struc-ture backbone and sulphonic acid side chains for the ion conduction purpose. Although proton move-ment in the membrane is not well understood, suit-able membrane water content is beneficial [2]. Additionally, with polar molecular property, water is carried by protons and moved from the anode to the cathode during cell operation. Concentration and pressure gradients through the membrane also have a degree of influence on water transport.

_{Corresponding author: Department of Mechanical Engineering,}

National Chiao Tung University, Room 342, Engineering Build-ing No. 5, 1001 Ta Hsuch Road, Hsinchu, Taiwan 300, Republic of China. email: sam@nkc.edu.tw

Water is produced in the cathode catalyst layer (CL) because of cell reaction and, with the aforemen-tioned effects, may lead to a flooding problem and dramatically reduce cell performance at high current density. Several works performed on this issue pro-pose a better understanding and fulfilling of fuel cell system water management [3 –6]. The tempera-ture factor is inevitable in water management inves-tigation because water phase change closely relates to corresponding saturation pressure and relative humidity. The energy fraction released by heat during cell operation is as high as 50 per cent of the total energy used at high current densities, such that there is an equal amount of thermal energy produced accompanying the power gener-ated. The extra high temperature might lead to membrane dehydration and poor ion transport without imposing special effort. On the contrary, when the electrode temperature is high, vapour satu-ration pressure is elevated with reduced relative humidity, preventing water formation and a further flooding problem.

Transport and energy conversion processes in the PEMFC involve multi-phase, multi-component flow and heat and mass transfer in porous media with electrochemical reactions. Reactants are deliv-ered and distributed evenly by grooved flow channels on bipolar plates (BPs). Despite the simple geometry of flow channels, the actual flow field is intrinsically three-dimensional that contains longitudinal, trans-verse, and through plane directions. The cathode CL serves as oxygen sink and water source. The most challenging aspect of modelling the fuel cell transport phenomena is that these processes are coupled with complex electrochemical reactions.

Bernardi and [7], Bernardi and Verbrugge [8, 9], Springer et al. [10] pioneered the study of PEMFC modelling, which utilized the lump approach and the conservation theory to develop unidirectional models through electrodes and membrane assum-ing isothermal condition. The Bernardi model [7] suggested reactant gas humidification tuning to accommodate changing cell demands for water. Springer et al. [10] utilized their model to study the water transport coefficient and cell performance for different anode humidifier temperatures. Fuller and Neuman [11] considered a two-dimensional membrane – electrode assembly and examined the water thermal management in PEMFC. Membrane water concentration and gas channel tempera-ture were calculated for different heat transfer coeffi-cients. Findings showed that a lower heat transfer coefficient is responsible for membrane dehydration and substantial cell performance decrease. Yi and Nguyen [12] developed an along-the-channel, two-dimensional model including various thermal man-agement schemes. Their model results showed that

the heat exchanger design with counter flow is more effective than coflow and constant bulk temp-erature schemes. Besides these works, the model of Rowe and Li [13] included variable membrane hydration and water phase change in the electrodes. Baschuk and Li [14] investigated variable degree effect of water flooding in the CL layer and/or cath-ode electrcath-ode backing region on cell performance. Simpalee and Dutta [15] formulated a comprehen-sive three-dimensional model to predict temperature distribution inside a straight channel PEMFC and presented a group of source terms facilitating com-mercial software application. Djilali and Lu [16] ana-lysed transport phenomena in a PEMFC assessing non-isothermal and non-isobaric effects. Further-more, numerical results of Ju et al. [17] with a single-phase and non-isothermal model revealed that temperature rise within the cell could be as high as 13.75 K at medium reaction rate under certain conditions.

It is clear from the above literature reviews that water and thermal management plays an essential role on the PEMFC performance. However, no studies have ever investigated the gradient influences of cell temperature and humidification level arising from transient cold start or inadequate cooling con-ditions. The purpose of this study is to present a sys-tematic investigation of these two factor impacts on PEMFC polarization characteristics. Relationships between temperature and electrochemical kinetics, as well as membrane conductivity, are modelled using empirical formulas. Local temperature and relative humidity factors on water phase change and transport are harnessed in the modelling to evaluate their influences on cell outputs.

2 MATHEMATICAL FORMULATIONS

Various transport phenomena such as mass, momentum, and energy, as well as electrochemical reactions and phase change are coupled during the PEMFC power generation. Governing equations con-cerning these problems are introduced in the follow-ing. Figure 1 shows the physical and computational domain considered in this study consisting of BPs, gas diffusion layers (GDLs), CLs, and the PEM. Species contained in the anode region are hydrogen and water, whereas oxygen, nitrogen, and water are those in the cathode region.

The following assumptions are made for this study. 1. The calculated mixture compressibility factor is close to unity. Ideal gas behaviour is assumed for anode and cathode gas mixtures.

2. Owing to low reactant velocity, the Reynolds number is usually below 100 and the laminar

flow assumption is made for the fluid in the domain.

3. Transport properties in the porous media are isotropic and phase compositions in mixture are smoothly varying.

4. The electric potentials on top and bottom of the cell are constant.

5. The pressure diffusion, thermal diffusion, and Knudsen effect are not considered in species mass transport.

6. Working fluid and the porous media solid phase are in local thermal equilibrium.

7. There is no charge accumulation in the electrodes and the electroneutrality assumption is made in the domain.

2.1 Mass conservation of species and bulk mixture

The mass fraction of individual species in the domain can be described as

r (1rwaUlg) ¼ r  r Xn1 b¼1 Dabrwb ! þSs,a (1)

whereaandbdenote the individual species that rep-resent oxygen, hydrogen, nitrogen, and water, 1 is the effective porosity corrected by the saturation level,r

the bulk density of the mixture, and wathe mass frac-tion of speciesa. The equation source term connects the mass transport and electrochemical reaction in CLs. The consumption or creation of each species depends on the local volume-based current density divided by the Faraday constant and electron num-bers transferred. Note that the gas phase relative mobility lg is employed in the convection term of

equation (1) to consider the two-phase flow phenomenon.

The multi-component diffusion coefficient Dabis a function of species concentration and binary diffusion

coefficient, the latter can be written as [18, 19]

Dgd¼ 0:0018583 T3=2_(1=M gþ1=Md)1=2
 Pc_gdV_gd 1 d ₍₂₎

where T is the temperature in kelvin, P the pressure in atm, and Mgthe molecular weight of species g. The parameters in the denominatorcgdand Vgdare the col-lision diameter and the colcol-lision integral between speciesgandd.

The densities of species and bulk fluid have the relation

r¼X

n

a¼1

r_a(1  s) þrls (3)

where s is the saturation level of the liquid phase and the following continuity equation describes the mass conservation of mixture

r (1rU) ¼ 0 (4)

2.2 Momentum equation

The general form of Navier –Stokes equation is employed with source terms to account for drag forces arising in the porous region. The equation is given as

r (r1UU) ¼ 1rP þ r  (1meffr U) þ Sm (5) wheremeffis the effective viscosity of reaction mixtures in anode and cathode and Smrepresents the Darcy and

Forchheimer terms that describe the momentum dis-sipation in the porous media. Their expressions are

meff ¼ rls þrg(1  s) (krl=nl) þ (krg=ng) (6) Sm¼  12U ffiffiffiffi K p m eff ffiffiffiffi K p þ1CFrjUj   (7) Fig. 1 Physical and computational domain considered in this study

where K is the permeability and CFis the form drag

constant. Equation (5) is restored to the original Navier – Stokes Equation by setting porosity 1 equal to unity and K approach infinity in flow channel regions. At the other place, the convection term and diffusion term drop out because the per-meability in the denominator is extremely small. Employing a unified formulation in different trans-port property regions is advantageous, as there is no need of dealing with the interface condition of physical property continuity used in early works [17, 20].

2.3 Energy equation

Energy equation is important in this study, as it describes the thermal energy conservation in the domain and the temperature field can be extracted from it. Except for inlet and outlet boundaries, enthalpy generation is considered from three sources: joule heating due to current transport in the medium, electrochemical reaction irreversibility due to overpotential at electrodes, and latent heat of water due to phase change. Many physical proper-ties are closely related to temperature in the cell domain. First, mixture velocities at channel inlet ports for a fixed reactant stoichiometry and a given cell geometry are primarily decided by their thermo-dynamic states and subsequently dictated by mix-ture temperatures. Secondly, electrochemical reaction kinetics in the CL and membrane conduc-tivity are also influenced by local temperatures. Thirdly, the temperature field is also responsible for water formation or evaporation. Consequently, it has certain effects on reactant concentration and membrane water content. The generalized steady-state energy equation can be expressed as

r (1rhU) ¼ r  keffrT þX n a¼1 X n1 b¼1 rDabrwbha ! þ1t:rU Sjhþ ji  ij s þSh (8)

where the first two terms on the right-hand side rep-resent the conduction energy and reactant enthalpy flux and the third term is the irreversible viscous dis-sipation. The fourth and the fifth terms describe the electrical-related thermal effects, where Sj is the

volumetric current density in the CL that can be expressed by the Bulter –Volmer equation shown in section 2.5. The last term results from the consider-ation of phase change when water evaporconsider-ation or condensation occurs in the cell domain and is pre-sented in the next section following the phase change formulation.

2.4 Formation and transport of water

Water vapour partial pressure varies from channel inlet to outlet, as cell reactions alter its concen-tration. Water forms as a result, when local relative humidity is greater than 1. On the contrary, the water evaporation process activates as soon as rela-tive humidity is less than 1 and water exists. The gen-eralized Richards equation developed by Wang and Cheng [21] is employed to mimic water transport in the cell domain. This equation essentially orig-inates from water mass conservation in a control volume. It can be given as

r (1rUll) þ r  Nl ¼ r  1rlDcrs  lllgK (rlrg)g n   þSl (9)

In equation (a), the electro-osmotic drag effect is included using a divergent expression of water mass flux with a drag coefficient of 2.5 s [22]. The dif-fusion coefficient of water Dcis subjected to capillary pressure and, in turn, controlled by surface tension and water saturation level [23]

Dc ¼ s 3_{K (1 l} l)(dPc=ds) ml (10) Pc ¼6cosuc 1 K  1=2 ½1:417(1  s) 2:12(1  s)2þ1:262(1  s)3 (11) wherellis the liquid phase mobility given by relative

permeability and kinetic viscosity

ll¼

krl=nl krl=nlþkrg=ng

(12) Water condensation or evaporation is considered by introducing the source term Sl in equation (9)

as [24] Sl¼

Mlrcon1x_RTw(xwP  Psat), if xwP . Psat reva1srl(xwP  Psat), if xwP , Psat 

(13) where rconand revarepresent phase change rate

con-stants. Owing to this process, the source term Shin

equation (8) can be obtained by multiplying Slwith

water latent heat.

2.5 Cell potential and electrochemical reaction When current is drawn from the cell, two charged species move along individual paths by the potential gradient effect. Conducting materials for electrons are the solid phases (carbon cloth or paper, carbon support, and catalyst particle) in CLs, GDLs, and

BPs (flow field plate). Proton transport media are the ionomer phases in membrane and CLs. Potential fields in these two media can be described as

r (ssolrfsol) ¼ Sj (14) r (smemrfmem) ¼ Sj (15) wheressolandsmemare the electrical conductivities

of these two transport media andfsolandfmemare

the aforementioned potential fields. The term Sj

rep-resents the source or sink of the two currents. In the anode CL, they are created, and in the CL, they are consumed. The membrane conductivity smem is a

vital property that exerts great influence on the cell performance. The model of Springer et al. [10] is used in this study

smem¼(0:005139g0:00326)  exp 1268 1 303 1 T     (16) whereg is the membrane water content defined as the number of water molecules per sulphonic acid group and T is the temperature.

The CL overpotential is the electrochemical reac-tion driving force, which is the difference between

fsol and fmem. The volumetric-based transfer

cur-rent density considering electrode structure, reac-tant concentration, temperature, and overpotential is given by the Bulter – Volmer equation. It rep-resents the heterogeneous catalyst reaction rate among reactants and solid platinum particles in an electrode. At open circuit voltage, anodic and catho-dic currents occurred simultaneously. However, the overpotential biases the net current towards a favourite direction. Consequently, an oxygen reduction reaction in cathode or a hydrogen oxi-dation reaction in anode is activated. The equation is written as Sj,k¼Aviref,k ca ca,ref  wa  exp bk,aFhk RT   exp bk,cFhk RT   (17) where Av is the effective catalyst surface area per

unit volume of CL, iref,kthe reference exchange

cur-rent density, ca the molar concentration of species

a, and bk,a and bk,c are the anodic and cathodic

transfer coefficients at anode and cathode, respect-ively. Note that the CL contains nanoscale platinum particles distributed on carbon supports; its effective surface area is several orders of magnitude greater than the geometric area. The following expression provides the relation among Av and parameters

describing CL structure such as catalyst surface area per unit mass, catalyst loading, and layer thick-ness

Av¼

Amasmplt Hcat

(18) The experimental results from Parthasarathy et al. [25] provide the relation between the cathode refer-ence exchange current density and temperature

iref,o2¼104exp 3:507  4001 T   (19) 2.6 Boundary condition

As this study investigates temperature gradient and humidification level impacts on cell performance, specific wall temperatures are prescribed at the top and bottom of the computational domain. Thus, the uneven temperature scenario during system start up, as well as in a cell stack unit, could be simu-lated. At the channel outlet port, a fixed back pressure is prescribed. The reactant temperature at channel inlet is made equal to the humidification level, and the mixture species concentration can readily be calculated from Dalton’s law. With these data, anode or cathode mixture inlet velocity is decided by the following relation

vk¼

jIWcelLchaRTk,in

qFWchaHcha(1  Psat,w,k=Pk)xa,inPk

(20)

where the subscript a refers to hydrogen at anode and oxygen at cathode, respectively; and q equals 2 in anode flow and 4 in cathode flow.

Neumann boundary conditions are used on the symmetry walls and CL interfaces in contact with membrane and GDLs to ascertain thermal energy and electrical current flux impermeability. Fixed value of total cell overpotential is specified on the cathode outer boundary. This value includes ohmic overpotentials in electrodes and membrane, as well as activation overpotentials in CLs. With this method, membrane and solid-phase potential fields are obtained, and local current densities can be cal-culated from the ohmic law. Consequently, cell potential is obtained by subtracting total cell over-potential from open circuit voltage

fcel ¼fochtot (21)

foc¼0:2329 þ 0:0025Tcel (22) The feature of such voltage-to-current algorithm is that it can be used to resolve CL local activation

overpotential [26], needed in the calculation of the Bulter – Volmer equation.

2.7 Numerical procedure

The governing equations were converted to the con-servative convection – diffusion form and solved by the computational fluid dynamic algorithm with control volume approach [27 – 29], in which a stag-gered grid system was used for the velocity field. Momentum equations were solved first, and then the pressure and velocity fields were corrected to satisfy the continuity equation. Following the vel-ocity vector, temperature field, species concen-trations, potential fields, and saturation level were solved. Throughout the calculation procedure, proper under-relaxation factors were employed for numerical stability consideration. Main parameters and dimensions employed in this study are shown in Tables 1 and 2. The calculation process was iter-ated until the relative residual of each variable became smaller than 1  1024. Various grid systems with total cell numbers of 23 010, 53 600, 76 380, and 100 320 were constructed to evaluate the cell number dependence on output current density. The calculated data indicated that a cell number of 53 600 gave a satisfactory result, and the compu-tational time was largely reduced when compared with the other two higher cell number systems. Simulations were performed on a personal computer with a Pentium IV 2.4 GHz CPU. It generally required

48 h to accomplish a typical polarization curve with

13 data points.

Model validation certification was conducted by comparing results with Wang et al.’s report [30], using their model dimensions and operating con-ditions as shown in Fig. 2. It discloses a satisfactory agreement between experimental and current mod-elling results at different operating temperatures.

Table 2 Main physical parameters and properties utilized in this study

Parameters and properties Value Units Porosity of the diffusion and

catalyst layer

0.4 Porosity of the membrane 0.28 Permeability of the diffusion and

catalyst layer

2.3  10211 m2 Permeability of the membrane 1.0  10218 m2 Tortuosity of the diffusion and

catalyst layer

1.5 Tortuosity of the membrane 3 Molar ratio of cathode side

dry air (N2/O2)

79/21 Collision diameter between

H2and H2O

2.734 A˚

Collision diameter between O2and H2O

3.054 A˚

Collision diameter between O2and N2

3.633 A˚

Collision diameter between N2and H2O

3.220 A˚

Surface tension of water 0.06253 N/m

Contact angle 0 8

Condensation rate constant 100 s21

Evaporation rate constant 100 1/atm s Concentration dependence of H2 0.5

Concentration dependence of O2 1

Transfer coefficient at anode 0.5 Transfer coefficient at cathode 1.5

Electrical conductivity of electrode 114 1/V m Thermal conductivity of electrode 1.6 1/Wm K Thermal conductivity of membrane 0.3 1/Wm K Catalyst loading 4.0  1023 _{kg m}2

Catalyst surface area per unit mass 1.0  105 _m2_kg

Fig. 2 Comparison of current model results with Wang et al.’s report [30]. The humidification temperatures are 343 K and backpressures are 3 atm on anode and cathode sides

Table 1 Model dimensions and operating conditions of this study

Parameters and properties Value Units

Domain length 5.0  1022 _m

Domain width 1.6  1023 m

Gas channel width 8.0  1024 _m

Domain height 4.303  1023 m

Gas channel height 8.0  1024 _m

Diffusion layer thickness 2.54  1024 m Catalyst layer thickness 1.0  1026 _m

Membrane thickness 1.75  1024 _m

Molar ratio of cathode side dry air (N2/O2)

79/21 Anode H2stoichiometry 3

Cathode O2stoichiometry 3

Anode pressure 2 atm

3 RESULTS AND DISCUSSION

A series of numerical simulations of underlying fac-tors addressed in this study followed the modelling validation assessment. Prediction results presented are two-fold: those based on equal cell boundary temperatures on both sides while reactant humidifi-cation temperatures are altered and those based on equal temperatures of boundary and humidification on each side with a temperature gradient imposed between both sides. To examine the interrelation between the two main factors considered, humidifi-cation level and cell boundary temperatures varied from higher values (353 – 363 K) to lower values (333 to 323 K). This differs from previous studies with only one parameter changed [30, 31], and hence, a more vivid situation that simulates cell operation is expected.

Note that from previous section comparisons, cur-rent modelling captures satisfactory polarization curve characteristics at low or high temperature range. However, this study focuses the prediction in higher temperature region with perturbed values within 20 K because researchers are more interested in cell performance in this region [13 –17]. Also, owing to extreme thinness of the domain, introdu-cing excessive temperature discrepancy on both sides seems unusual and impractical.

3.1 Effects of humidification scheme

In regard to promoting proton transport through the membrane, proper water vapour concentration in the cell region is crucial. Enhanced or reduced humidification scheme results on one or both sides of channel inlet ports are first reported. Subscripts B, H, A, and C in the legend of Fig. 3 and the following figures represent boundary, humidification, anode, and cathode, respectively. A temperature difference of 10 K between cell and humidification level, either larger or smaller, is posted on designated positions.

Results between the anode-enhanced humidifi-cation scheme and the base cases with these four temperatures being equal are compared in Fig. 3. Three cell boundary temperatures 333, 343, and 353 K are considered. The data in Fig. 3 reveal a con-centration polarization effect for this humidification scheme in higher current density. In addition, this decay onset depends on cell temperature, and at higher cell temperature this phenomenon appears earlier. As inlet gas is fully saturated depending on to the humidification level, the anode mixture temp-erature decreases from heat transfer with the cell solid phase, and saturation level in the anode is elev-ated because of water condensation. This causes

effective transport passage decrease of anode gas, promoting concentration overpotential. Another contributor to this is hydrogen depletion from the oxidation reaction in the anode CL, such that water partial pressure increases as well as condensation probability.

Cathode-enhanced humidification scheme influ-ence with cathode mixture remained at 10 K greater than cell temperature as shown in Fig. 4. The results reveal a more serious concentration overpotential than the previous scheme; and the limiting current density decreases at higher cell temperature, i.e. cell performance dependence on temperature reversed at a high reaction rate for the present humidification scheme. Higher humidification temperature and corresponding low oxygen concen-tration in the inlet port may be considered to cause extra concentration overpotential at high cell temp-erature for this humidification scheme. However, reactant gas inlet velocity adjusts to alleviate dilution caused by high vapour concentration and satisfies the required stoichiometry as indicated in equation (20). In this context, high humidification level and water generation combined effect in the cathode reaction lead to a stronger effective porosity decay in the cathode electrode and decrease limiting current density at high cell temperature.

The simulation results for the anode-reduced humidification scheme are presented in Fig. 5. In this case, effect is found on the larger performance decay, compared with the base cases. This cell Fig. 3 Effect of anode-enhanced humidification scheme with 10 K of anode humidification level greater than cell temperature for cell temperatures 333, 343, and 353 K

voltage reduction attributes to water swell deficiency in the membrane because of lower inlet water con-centration. As the membrane conductivity is dictated by water content – the number of water molecules per fixed charged site – the supplied water vapour

from anode inlet port is not enough to keep the membrane fully hydrated, such that the membrane proton transport is more difficult and activates a larger ohmic overpotential.

The cathode-reduced humidification scheme effect on cell output is shown in Fig. 6 for various cell temperatures. Comparing with base cases, it reveals the merit of applying the current scheme at a higher current density, causing cell voltage and performance increase at each cell temperature. Gas passage clogging prevention in the electrode is responsible for this gain because the produced water in the cathode reaction balances the lower inlet water concentration. However, an additional study (not shown in the figure) revealed that it reaches a plateau above cell temperatures of 348 K as large amounts of water generation exceed reduced humidity. For further improved performance, a larger humidification temperature reduction is needed to influence the membrane water content. In addition, at medium reaction rate and high cell temperature, cell resistance shows a slight variation. This is because the variation of membrane water content and conductivity arise from lower cathode water concentration.

To examine the impact of greater cell humidifi-cation gradient, the combined results of anode-reduced and cathode-enhanced humidification schemes are depicted in Fig. 7. The ohmic and con-centration overpotentials clearly are all altered, com-pared with base cases. Cell performance in this

Fig. 5 Effect of anode-reduced humidification scheme with 10 K of anode humidification level smaller than cell temperature for cell temperatures 333, 343, and 353 K

Fig. 6 Effect of cathode-reduced humidification scheme with 10 K of cathode humidification level smaller than cell temperature for cell temperatures 333, 343, and 353 K

Fig. 4 Effect of cathode-enhanced humidification scheme with 10 K of cathode humidification level greater than cell temperature for cell temperatures 333, 343, and 353 K

operating condition can be regarded as the individ-ual scheme superposition. Nevertheless, the plot for 353 K exhibits a smaller decay at high reaction rate than the result in Fig. 4, as the reaction rate at the conduction-dominated region is suppressed from low anode humidity and the accompanying small membrane conductivity. Consequently, the series concentration overpotential at low cell voltage is slightly released. In a further study that reverses the humidification condition of the previous case, the same expected individual scheme superposition is shown. However, this effect is intrinsically non-linear because of complicated mechanisms of the coupled transport phenomena.

3.2 Effects of cell temperature gradient

Increasing cell bulk temperature reportedly creates a positive influence on cell performance [30]. Predic-tion results containing a series of boundary and humidification temperature variation on one elec-trode with others fixed are reported in this section. The reactants are also fully saturated based on boundary temperature. Polarization curves are plotted as total cell overpotential versus current den-sity as cell temperature used in equation (22) cannot be obtained. Results with anode temperatures increased from 333 to 353 K at 353 K cathode temp-eratures are presented in Fig. 8. Cell current density does not change significantly until a higher reaction rate according to findings. Membrane conductivity

along central channel direction on the membrane middle plane is drawn in Fig. 9 at total cell overpo-tential of 0.98 V for further result information. It shows two influential factors for conductivity; the first is cell temperature with a positive effect. How-ever, owing to the electro-osmotic drag transport process, water concentration decreases along the Fig. 7 Effect of humidification gradient that combines

anode-reduced and cathode-enhanced schemes for cell temperatures 333, 343, and 353 K

Fig. 8 Performance curves with 353 K cathode temperatures and anode temperatures varied from 333 to 353 K

Fig. 9 Membrane conductivity along central channel direction on the membrane middle section for cases in Fig. 8

anode channel, leading to a quicker membrane water content decay at 333 K anode temperature. Conse-quently, its current density exhibits a sensible decrease when compared with other higher anode temperature situations. A similar scenario with 333 K cathode temperature was also calculated. However, no appreciable cell performance variation is revealed. This can be explained by the anode kin-etic independence with temperature, which leads to trivial cell output change.

The prediction result at 353 K anode temperature and continuous cathode temperature variation from 333 to 353 K is represented in Fig. 10. Obviously, the output current density increased monotonously with cathode temperature change, indicating that cell behaviour is dominated by the electrochemical reaction kinetics rather than conduction or con-centration overpotentials. With sluggish cathode oxygen reduction, cathode exchange current density is several orders of magnitude smaller than the anode value, and a higher cathode temperature is beneficial for cell reaction and electricity generation. Temperature gradient effect at low anode temp-eratures of 333 K and cathode temperature variation from 333 to 353 K are shown in Fig. 11. It indicates that current density dependence on cathode temp-erature variation has an entirely different tendency between medium and high cell reaction rates. Mem-brane conductivities along the central channel direc-tion on middle secdirec-tion for total cell overpotentials of 0.5 V and 0.98 V are plotted in Fig. 12 for further information. When total cell overpotential is at

0.5 V, membrane conductivities are shown to be smaller than those of 0.98 V for various cathode temperatures. This can be explained by the smaller membrane water content for the former case with a weaker cell reaction and electro-osmotic drag. In addition, owing stronger cell reaction and water transport, membrane conductivity along the cell Fig. 11 Performance curves with 333 K anode temperatures and cathode temperatures varied from 333 to 353 K

Fig. 10 Performance curves with 353 K anode temperatures and cathode temperatures varied from 333 to 353 K

Fig. 12 Membrane conductivity along central channel direction on the membrane middle section for cases in Fig. 11 at total cell overpotentials of 0.5 and 0.98 V

channel has an earlier and more serious decay for overpotential at 0.98 V and 353 K cathode tempera-ture. For this reason, it reverses the intrinsic relation between cell temperature and performance at high current density.

3.3 Local distributions of physical properties Various temperature and humidification gradient effects on cell performance are previously presented. Displaying local physical property variations such as reactant concentration, temperature, and saturation level offers a more insightful understanding of the coupled transport phenomena in the PEMFC.

Temperature contours in cathode GDL at the x ¼ 0.025 m section with 353 K anode temperatures and 343 and 333 K cathode temperatures are shown in Fig. 13. The total cell overpotential is at 0.8 V. The upper region of the plot is adjacent to the cath-ode channel and shoulders, whereas the lower boundary of the plot is connected to the CL. Obviously, as CL is the main source of heat gener-ation, the highest temperature appears at the

bottom of GDL. Meanwhile, on this CL interface, the central region exhibits greater temperature as the cathode BP and shoulder have more capability to dissipate heat than the channel.

The expected effect on cell physical properties is justifiable, given the temperature distribution. Local oxygen mass fractions for the previous cases are shown in Fig. 14. The data in this figure indicate that a higher GDL and CL temperature corresponds to lower oxygen concentration. This can be attribu-ted to the higher reaction rate from equation (19) and greater membrane conductivity. More oxygen is consumed in the cell reaction as a result, and the local concentration is reduced.

Another temperature field influence can be observed in Fig. 15, which is the water saturation level in the GDL for the same scenarios. There are two important points that can be observed from this plot. First, as the transport direction is opposite to that of oxygen, the GDL shoulder region tends to accumulate more water vapour than the GDL chan-nel region, leading to a higher saturation level at this location. Secondly, water vapour saturation pressure decreases at lower temperature, such that local water vapour partial pressure is apt to exceed saturation pressure. Consequently, the local satur-ation level increases with a 333 K cathode tempera-ture. The further influence of this phenomenon is that reactant and product transport passages are hindered and the concentration overpotential is elevated.

4 CONCLUSIONS

A three-dimensional, non-isothermal numerical model has been developed to explore temperature gradient and humidification perturbation effects on Fig. 13 Local temperature contours in cathode GDL at

section of x ¼ 0.025 m with 353 K anode temperatures and cathode temperatures (a) 343 K and (b) 333 K

Fig. 14 Local oxygen mass fractions in cathode GDL at section of x ¼ 0.025 m with 353 K anode temperatures and cathode temperatures (a) 343 K and (b) 333 K

Fig. 15 Local water saturation in cathode GDL at section of x ¼ 0.025 m with 353 K anode temperatures and cathode temperatures (a) 343 K and (b) 333 K

PEMFC. This model accounts for multi-component species transport and the two-phase flow. The CL local activation overpotential resolution is also achieved. According to the model prediction results and discussion, the following conclusions can be drawn.

1. The anode humidification level bears the impact of concentration and ohmic overpotentials on cell performance at elevated and reduced situa-tions, respectively.

2. Proper cathode humidification level-reduction benefits cathode reactant transport and improves concentration polarization at a high reaction rate. 3. Despite its non-linear behaviour, a suitable PEMFC humidification scheme can be reached by carefully using the individual humidification level effect on either side of the cell.

4. Temperature gradient effect on cell performance exhibits entirely different behaviours depending on its direction and magnitude. When the cathode temperature is fixed at a lower level, the anode temperature variation has no noticeable change in the output current. However, when it is speci-fied at a higher value, a larger temperature gradi-ent is harmful at high currgradi-ent density.

5. Positive cell temperature gradient effect on its performance is found for medium reaction rate at a lower anode temperature. However, this trend is reversed at a higher reaction rate because of membrane dehydration.

6. Local physical property contours demonstrate interrelations among the temperature field, local oxygen concentration, and water saturation level. As a result, cell performance has a close con-nection with these property variations, and the essential role of existing PEMFC temperature and humidification gradients is elucidated.

REFERENCES

1 Spiegel, R. J., Gilchrist, T., and House, D. E. Fuel cell bus operation at high altitude. Proc. Instn Mech. Engrs, Part A: J. Power and Energy, 1999, 213, 57 – 68. 2 Nguyen, T. and White, R. E. A water and heat

manage-ment model for proton-exchange-membrane fuel cells. J. Electrochem. Soc., 1993, 140(8), 2178– 2186.

3 Zawodzinski, T. A., Derouin, C., Radzinski, S., Sherman, R. J., Smith, V. T., Springer, T. E., and Gottesfeld, S. Water uptake by and transport through Nafion 117 membrane. J. Electrochem. Soc., 1993, 140(3), 1041 –1047.

4 Siegel, N. P., Ellis, M. W., Nelson, D. J., and von Spakovsky, M. R. A two-dimensional computational model of a PEMFC with liquid water transport. J. Power Sources, 2004, 128, 173 –184.

5 Voss, H. H., Wilkinson, D. P., Pickup, P. G., Johnson, M. C., and Basura, V. Anode water removal: a water

management and diagnostic technique for solid poly-mer fuel cells. Electrochim. Acta, 1995, 40(3), 321 – 328. 6 Wood, D. L., III, Yi, J. S., and Nguyen, T. V. Effect of direct liquid water injection and interdigitated flow field on the performance of proton exchange mem-brane fuel cells. Electrochim. Acta, 1998, 43(24), 3795 – 3809.

7 Bernardi, D. M. Water-balance calculation for solid-polymer-electrolyte fuel cells. J. Electrochem. Soc., 1990, 137(11), 3344 – 3350.

8 Bernardi, D. M. and Verbrugge, M. V. Mathematical model of a gas diffusion electrode bounded to a poly-mer electrolyte. AIChE J., 1991, 37(8), 1151 – 1163. 9 Bernardi, D. M. and Verbrugge, M. W. A mathematical

model of the solid-polymer-electrolyte fuel cell. J. Electrochem. Soc., 1992, 139(9), 2477– 2491.

10 Springer, T. E., Zawodzinski, T. A., and Gottesfeld, S. Polymer electrolyte fuel cell model. J. Electrochem. Soc., 1991, 138(8), 2334 –2342.

11 Fuller, T. F. and Neuman, J. Water and thermal man-agement in solid-polymer-electrolyte fuel cells. J. Elec-trochem. Soc., 1993, 140(5), 1218 – 1225.

12 Yi, J. S. and Nguyen, T. V. An along-the-channel model for proton exchange membrane fuel cells. J. Electro-chem. Soc., 1998, 145(4), 1149 – 1159.

13 Rowe, A. and Li, X. Mathematical modeling of proton exchange membrane fuel cells. J. Power Sources, 2001, 102, 82 – 96.

14 Baschuk, J. J. and Li, X. Modeling of polymer electro-lyte membrane fuel cells with variable degree of water flooding. J. Power Sources, 2000, 86, 181 – 196.

15 Simpalee, S. and Dutta, S. Numerical prediction of temperature distribution in PEM fuel cells. Numer. Heat Trans. Part A, 2000, 38, 111 –128.

16 Djilali, N. and Lu, D. Influence of heat transfer on gas and water transport in fuel cells. Int. J. Therm. Sci., 2002, 41, 29 – 40.

17 Ju, H., Meng, H., and Wang, C. Y. A single-phase, non-isothremal model for PEM fuel cells. Int. J. Heat Mass Trans., 2005, 48, 1303 – 1315.

18 Wangard, W., III, Dandy, D. S., and Miller, B. J. A numerical stable method for integration of the multi-component species diffusion equations. J. Comp. Phys., 2001, 174, 460 –472.

19 Cussler, E. L. Diffusion mass transfer in fluid system, 2nd edition, 1997 (Cambridge University Press, Cambridge).

20 Gurau, V., Barbier, F., and Liu, H. T. An analytical sol-ution of a half-cell model for PEM fuel cell. J. Electro-chem Soc., 2000, 147(7), 2468 –2477.

21 Wang, C. Y. and Cheng, P. Multiphase flow and heat transfer in porous media. In Advance in heat transfer (Eds J. P. Hartnett, T. F. Irvine, Jr, Y. I. Cho, and G. A. Greene), 1997 (Acadmic Press, San Diego). 22 Mazumder, S. and Cole, J. S. Rigorous 3-D

mathemat-ical modeling of PEM fuel cells II. Model predictions with liquid water transport. J. Electrochem. Soc., 2003, 150(11), A1510 – A1517.

23 Wang, C. Y. and Beckermann, C. A two-phase mixture model of liquid – gas flow and heat transfer in capillary porous media-I. Formulation. Int. J. Heat Mass Trans. 1993, 36(11), 2747 – 2758.

24 He, W., Yi, J. S., and Nguyen, T. V. Two-phase flow model of the cathode of PEM fuel cells using interdigi-tated flow fields. AIChE J., 2000, 46(10), 2053 – 2064. 25 Parthasarathy, A., Srinivasan, S., and Appleby, A. J.

Temperature dependence of the electrode kinetics of oxygen reduction at the platinum/nafion interface-a microelectrode investigation. J. Electrochem. Soc., 1992, 139(9), 2530 – 2537.

26 Nguyen, P. T., Berning, T., and Djilali, N. Compu-tational model of a PEM fuel cell with serpentine gas flow channels. J. Power Sources, 2004, 130, 149 –157. 27 Hwang, J. J., Chen, C. K., Savinell, R. F., Liu, C. C., and

Wainright, J. A three-dimensional numerical simu-lation of the transport phenomena in the cathodic side of a PEMFC. J. Appl. Electrochem., 2004, 34, 217 – 224.

28 Hu, G., Fan, J., Chen, S., Liu, Y., and Cen, K. Three-dimensional numerical analysis of proton exchange membrane fuel cells (PEMFCs) with conventional and interdigitated flow fields. J. Power Sources, 2004, 136, 1 – 9.

29 Patanker, S. V. Numerical heat transfer and fluid flow, 1st edition, 1980 (Hemisphere, New York).

30 Wang, L., Husar, A., Zhou, T., and Liu, H. A parametric study of PEM fuel cell performances. Int. J. Hydrogen Energy, 2003, 28, 1263 – 1272.

31 Berning, T. and Djilali, N. Three-dimensional compu-tational analysis of transport phenomena in a PEM fuel cell – a parametric study. J. Power Sources, 2003, 124, 440 – 452.

APPENDIX Notation

A surface area of the catalyst particle (m2) c molar concentration (mol/m3)

C form drag constant

D multi-component diffusion coefficient (m2/s)

D binary diffusion coefficient (m2/s) F Faraday constant (96 500 C/mol) g gravitational acceleration (m/s2) h enthalpy (J/kg)

H height (m)

i vector of current density (A/m2) i current density (A/m2)

I reference current density (A/m2) k relative permeability

K permeability (m2) L length (m)

m catalyst loading (kg/m3) M molecular weight (kg/mol) n number of species

N flux (kg/m2/s) P pressure (Pa)

q number of electron transferred r rate constant (s21or 1/atm s) R universal constant (8.314 J/mol/K)

s saturation level S source term T temperature (K) U fluid velocity (m/s) v inlet velocity (m/s) w mass fraction W width (m) x molar fraction b transfer coefficient

g water content in membrane

d tortuosity 1 porosity 6 surface tension (N/m) h overpotential (V) u contact angle (8) k thermal conductivity (W/m/K) l relative mobility m dynamic viscosity (N s/m) n kinematic viscosity (m2/s)

j stoichiometry of hydrogen or oxygen

r density (kg/m3)

s electrical conductivity (1/V m) t shear stress tensor (N/m2)

f potential (V) w concentration dependence c collision diameter (A˚ ) V collision integral Superscripts c capillary eff effective value

Subscripts

a anodic

BA boundary of anode side BC boundary of cathode side c cathodic cat catalyst cel cell cha channel con condensation eva evaporation F Forchheimer term g gas h enthalpy

HA humidification of anode side HC humidification of cathode side in inlet

j volumetric transfer current density k anode or cathode

m momentum

mas per unit mass of catalyst carbon mem membrane phase of the cell n hydrogen or oxygen

oc open circuit value plt platinum

r relative

ref reference s species sat saturation

sol solid phase of the cell tot total

v per unit volume