Two-dimensional modeling of electrochemical and transport phenomena in the porous structures of a PEMFC

international journal of hydrogen energy 34 (2009) 3091–3103

Available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/he

Two-dimensional modeling of electrochemical and transport phenomena in the porous structures of a PEMFC Melik Sahraouia, Chafik Kharratb, Kamel Halouanib,* a

Institut Pre´paratoire aux Etudes d’Inge´nieurs de Tunis (IPEIT), Tunisia UR: Micro-Electro-Thermal Systems (METS-ENIS), Industrial Energy Systems Group, Institut Preparatoire aux Etudes d’Ingenieurs de Sfax (IPEIS), University of Sfax, B.P: 1172, 3018 Sfax, Tunisia

b

article info

abstract

Article history:

A two-dimensional CFD model of PEM fuel cell is developed by taking into account the

Received 1 October 2007

electrochemical, mass and heat transfer phenomena occurring in all of its regions simul-

Received in revised form

taneously. The catalyst layers and membrane are each considered as distinct regions with

3 November 2008

finite thickness and calculated properties such as permeability, local protonic conductivity,

Accepted 4 November 2008

and local dissolved water diffusion. This finite thickness model enables to model accu-

Available online 18 February 2009

rately the protonic current in these regions with higher accuracy than using an infinitesimal interface. In addition, this model takes into account the effect of osmotic drag in the

Keywords:

membrane and catalyst layers. General boundary conditions are implemented in a way

PEM fuel cell

taking into consideration any given species concentration at the fuel cell inlet, such as

Numerical modeling

water vapor which is a very important parameter in determining the efficiency of fuel cells.

Finite volume

Other operating parameters such as temperature, pressure and porosity of the porous

Heat

structure are also investigated to characterize their effect on the fuel cell efficiency.

Mass and charge transfer

ª 2008 International Association for Hydrogen Energy. Published by Elsevier Ltd. All rights reserved.

1.

Introduction

As the human needs for energy increase and with them the prices of fossil fuels, it is becoming vitally essential to develop alternative energy sources with minimum negative impact on the environment. The fuel cell technology has been identified as one of the alternatives that offers better efficiency than the combustion of hydrocarbons while emitting less pollution. The fuel cell is an electrochemical device that converts chemical energy directly into electricity with heat and water as byproducts. The Proton Exchange Membrane Fuel Cell (PEMFC) is one of the fuel cell technologies that is widely used because of its lower operating temperature with applications ranging from powering vehicles to mini and micro devices such as laptops or cell phones.

In recent years a lot of progress has been achieved in the development of cost effective fuel cell technologies, which enabled more widespread use of this technology. However, more research in design optimization and improved materials needs to be done in order to make the fuel cell technology ubiquitous. For the design optimization, modeling of the transport phenomena and chemical reactions for two- or three-dimensional representations of the fuel cells are very critical in the development of efficient fuel cells. Recently research on PEMFC technology has received a great deal of attention focusing on materials, size optimization, water and thermal management, and reliability. In studying transport and electro-kinetics, many researchers have focused on different aspects of the PEMFC and the approaches in dealing with the problem are quite varied. Fuel

* Corresponding author. Tel.: þ216 74 241 223; fax: þ216 74246347. E-mail address: [email protected] (K. Halouani). 0360-3199/$ – see front matter ª 2008 International Association for Hydrogen Energy. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.ijhydene.2008.11.012

3092

international journal of hydrogen energy 34 (2009) 3091–3103

Nomenclature

yi

mass fraction

A c cp D e F H I j ha hc K k L M p R S T u u v Vth Vcell

Greek l 3 4 s m r

water content porosity potential protonic conductivity dynamic viscosity density

area concentration specific heat mass diffusivity Thickness Faraday’s constant fuel cell overall height current density volumetric current density anode height cathode height permeability thermal conductivity fuel cell length molar mass pressure universal gas constant source term temperature x-direction velocity average inlet velocity y-direction velocity theoretical voltage working cell voltage

cell models can be categorized as analytical, semi-empirical or numerical as summarized by Chedie and Munroe [1]. Simple analytical models limited to one-dimensional representations, such as the ones developed by Standaert et al. [2,3] and by Charradi and Halouani [4], are able to predict the cell polarization curve for a given operating current density using ohm’s law and Butler–Volmer equation. However, analytical models are unable to predict accurately transfer processes inside the cell but remain useful for quick calculations of simple systems. More detailed analysis of the multiple transport and electrochemical mechanisms requires the use of numerical models (CFD). Recently numerical modeling from one- to three-dimensional representations of fuel cells has received the most attention. Bernardi and Verbrugge [5] presented a steady state, one dimensional, isothermal model for a fully hydrated membrane. This model studied the species transport in the gas, liquid and solid phases, the influence of the porosity of the electrodes and the effects of membrane properties. They identified inefficiencies due to unused species and low catalyst utilization at practical operating current densities. You and Liu [6] presented an isothermal one-dimensional model of the cathode catalyst layer. Their model predicts the current limit condition at high current densities. Their model shows that the porosity and polymer content in the catalyst layer has to be optimized in order to produce the highest current density. Fuller and Newman [7] presented a transient, two dimensional, non-isothermal membrane electrode assembly model. The model showed that the performance of these cells improves when gas streams are saturated with water. Gurau et al. [8] used a steady state, two dimensional, non-isothermal

Subscripts 0 reference value Amb ambient a anode c cathode CL catalyst layer drag frictional drag DW dissolved water e electronic i,j species (H2, O2, H2O, N2) hs heatsink source term m membrane ohm ohmic eff effective rev reversible th theoretical WV water vapor

model describing regions of the gas channels, gas diffusers, catalyst layers and membrane. Their model shows that a nonuniform reactant distribution has an important impact on the current density. But, this model is unable to predict the concentration overpotential due to mass transport limitation in the catalyst layer at high current densities. Eaton [9] used a transient, one-dimensional model for heat, mass and charge transfer in a PEM (Proton Exchange Membrane). Through this model, which accounts for a pressure gradient across the fuel cell, the authors concluded that the water flux from anode to cathode increases with current. But the water flux may be offset by the application of a pressure gradient between the anode and cathode. Genevey [10] presented a transient, one-dimensional comprehensive model of the cathode catalyst layer. The model is based on an agglomerate catalyst layer structure and includes the transport of thermal energy, gas species, liquid water, and charge transfer. Nguyen et al. [11] developed a three-dimensional model including the bi-polar plates with serpentine gas flow channels. Their model does not take into account the liquid water phase and the local overpotential on the cathode side is calculated locally from the potential equation. The overpotential on the anode side is neglected. Their model showed the current non-uniformity due to the presence of the solid plates and it is found that local maximal values for current production change with cell usage conditions. Ramousse et al. [12] developed a steady state and onedimensional model for heat and mass transfer in a PEMFC. Their results show that temperature gradient in the membrane can be affected by temperature difference of the inlet gases and this effect can be much bigger than the

international journal of hydrogen energy 34 (2009) 3091–3103

gradient due to resistive heating. Their results also show that the temperature is greatly affected by the thermal conductivity of the membrane. Hmidi and Halouani [13] developed a steady state, one-dimensional numerical model of a PEMFC studying the effect of temperature, pressure and current density. This model showed that the water content increases with current density but decreases with an increase of temperature or pressure whereas the membrane potential increases with both temperature and pressure. Siegel et al. [14] developed a two-dimensional model for PEMFC with a modified Butler–Volmer model in order to take into account the agglomerate geometry and take into account the pore level mass transfer limitation on the performance of the fuel cell. Their model uses the Thiele’s modulus to incorporate the pore level mass transfer limitation. Siegel et al. [15] extended the previous work to include the effect of liquid water and have improved further the electrochemical model. Their results show that the agglomerate geometry and liquid water influence significantly the efficiency of the fuel cell. Sivertsen and Djilali [16] developed three-dimensional model for PEMFC using a commercial CD tool (Fluent) with user defined functions (UDF) to incorporate the electrochemistry in the membrane-electrode assembly. In their model they assumed that the membrane is fully humidified, therefore, limiting the validity of the model to high inlet relative humidity. Phase change for water vapor is not taken into account. The results produced by their model were in good agreement with the measurements. However, the model does not predict, the accurate variation of the slope of the polarization curve when the fuel cell is operating at high current density and near the current limit. Um and Wang [17] examined water transport in PEMFC using a 3D model. Their model assumed that only water vapor is present and they used it to solve for the gaseous phase of water throughout the entire MEA. In the polymer they used a fictive gaseous phase assumed to be in equilibrium with the water content. With their model they examined the effect of humidity ratio of inlet mixtures for different values of membrane thickness, and for co-flow and counter flow conditions. Their results show with the counter flow design and at low humidity ratio, the current density produced is comparable to that of the co-flow condition with fully humidified inlet mixtures. Wu et al. [18] presented an optimization methodology using design of experiments and a few results from a CFD three-dimensional tool in order to find the best solution which maximizes power density. Their results show that many solutions can be found. This approach avoids running a large number of three-dimensional solutions, which may take a long time to get results. Other studies have examined the optimization of PEMFC by proposing different flow configurations in the channel by using flow baffles (Liu et al. [19]), variables cross section (Liu et al. [20] and Yan et al. [21]). These studies have shown that improved performance can be obtained with the above modifications since they help in spreading reacting species from the channel to the catalyst layer. Harvey et al. [22] using a single phase model for a threedimensional representation of a PEMFC examined three modeling approaches for the catalyst layers which are the agglomerate model as discussed by Siegel et al., the discrete

3093

volume where the pores are assumed to be filled with water, and the third is the thin film model where the catalyst layer is used as a boundary condition. Their results indicate that the best model that predicts well the experiments is the agglomerate model, which can predict the sharp drop off in voltage when the current limit is reached. Kamarajugadda and Mazumder [23] used a two-dimensional single phase model to investigate different models for membrane conductivity and they found that the choice of the membrane model can affect the current density results significantly. In this work, a mathematical model of PEMFC is developed in order to build a numerical tool enabling the design optimization of PEMFC for a wide range of parameters. This model introduces a different approach in dealing the pore level mass transfer than previously presented as in Siegel et al. [14,15]. This model uses in the Butler–Volmer equations governing the electrochemical source terms a fictitious pore level concentration instead of the volume averaged concentration. The fictitious concentration on the surface of the catalyst particle allows to take into account the pore level mass transfer around the catalyst particles and predict the current limit condition. The pore level model is based on the local solution of the mass transfer equation in the polymer and around a spherical catalyst particle. Using the presented model, the performance of the fuel cell is presented for many parameters such as species inlet pressure, relative humidity, temperature, membrane thickness, and catalyst layer porosity.

2.

Overview of fuel cells

A proton exchange membrane fuel cell (PEMFC) is an electrochemical device that converts the chemical energy of a hydrogen oxidation reaction directly into electrical energy with production of heat and water. A schematic representation of a PEMFC is shown in Fig. 1 with all its components, namely:  The gas channel (GC), which is machined into the collector plate, serves as channels for the reacting gas species.  The gas diffusion layer (GDL) is made of a porous material such as carbon cloth or carbon paper. The porous nature of the GDL facilitates reactant distribution across the catalyst layer while being a good electrical conductor, providing a low electrical resistance connection between the catalyst

Fig. 1 – Computational domain.

3094

international journal of hydrogen energy 34 (2009) 3091–3103

layer and the collector. Also the GDL facilitates the removal of liquid water produced at the cathode and heat produced in the membrane and catalyst layers.  The catalyst layer (CL) main role is to accelerate the electrochemical reaction. It is a mixture of the three different components, namely membrane polymer, GDL, and catalyst particles supported by the structure. Platinum is the typical catalyst used in PEMFC. The catalyst enables the chemical reaction to occur. The overall reaction in the anode CL is given by:H2 / 2Hþ þ 2eand in the cathode CL by: 2Hþ þ 2e þ 1/2O2 / H2O. The overall cell reaction is given by: H2 þ 1/2O2 / H2O þ Heat.  The proton exchange membrane (PEM), sandwiched between the two fuel cell electrodes, is a polymer network containing covalently bonded negatively charged functional groups capable of exchanging ions. The polymer matrix consists of polymer backbone, which can be hydrocarbonbased polymers such as polystyrene and polyethylene or their fluorinated polymer analogs. For fuel cell applications, sulfonic acid functional group ðSO1 3 Þ is the most widely used.

3.

Numerical model

The two-dimensional representation of a PEMFC used in this study is shown in Fig. 1. In this study, the single domain modeling approach is used. The continuity, momentum, energy, and species transfer equations for all regions of the computational domain of the fuel cells can be written in the same manner. However, the diffusion, convection, and source terms will vary depending on the local properties of the region (e.g., GC, GDL, CL, and PEM). As an example, the source term for hydrogen consumption is activated in the anode CL only, whereas consumption of oxygen and production of H2O is activated only in the cathode CL. The PEMFC model is reduced to the following transport equations for the entire domain with variable diffusion and source terms which will be described below for each region. The continuity equation is given by: vru vrv þ ¼ SH2 þ SO2 þ SH2 O þ Sd vx vy

The protonic potential equation is given by:     v vf v vf sm m þ sm m þ Sp ¼ 0 vx vy vx vy

(7)

The model assumes: Steady state; Compressible and laminar flow; Newtonian fluid; Ohmic losses due to electrons transfer are neglected; Liquid water is not taken into account; Dispersion in the porous media is neglected; Each region is isotropic and homogeneous; Local thermal equilibrium between the porous solid structure and the fluid (GDL, CL, PEM);  H2 is present in the anode only, O2 and N2 are present in the cathode only, and H2O vapor is present in the entire computational domain except in the membrane.        

3.1.

The gas channel

In the gas channels, the continuity Eq. (1) has no source terms in both of these regions since there are no chemical reactions occurring. The momentum Eqs. (2) and (3) in both channels do not have the darcean terms since they are not porous. The source terms for the energy, mass, and potential are also nil in these regions. The properties of the gaseous mixture in all regions of the computational domain are calculated by taking into account the different species present locally. The mixture density, viscosity and mass diffusion are respectively given by: X Mi ci (8) r¼ i

X ri

X Mi ci

mi

(9)

  X yj Di;mix ¼ 1  yi = Dij jsi

(10)

m¼

i

r

mi ¼

i

r

and

where yi is the molar fraction of species i is given by: (3)

The Darcy term is included in order to account for the porous structure of the GDL and CL. The energy equation is given by:      v     v   v vT v vT 3 rcp f uT þ 3 rcp f vT ¼ keff þ keff vx vy vx vx vy vy (4) þ Srev þ Sohm þ Shs The species equation is given by:

(5)

The diffused water equation in the polymer phase is given by:     v vcDW v vcDW þ þ Sdrag þ Sd ¼ 0 (6) DDW DDW vx vy vx vy

(1)

The momentum equations in the x- and y-directions with the Darcy term included are given by:     vruu vrvu vp v vu v vu 3m þ ¼ þ m þ m  u (2) vx vy vx vx vx vy vy K     vruv vrvv vp v vv v vv 3m þ ¼ þ m þ m  v vx vy vy vx vx vy vy K

    vuci vvci v vci v vci þ ¼ þ þ Si Di;eff Di;eff vx vy vx vy vx vy

ci yi ¼ P cj

(11)

j

The diffusion coefficient given in Eq. (10) is based on a simplification of the Stefan–Maxwell equations [8]. The binary diffusion coefficient Dij for two species i and j is scaled using the local temperature and pressure [14]:    p0 T 1:5 Dij ðT; pÞ ¼ D0ij T0 ; p0 p T0

(12)

3095

international journal of hydrogen energy 34 (2009) 3091–3103

The specific heat and thermal conductivity of the fluid mixture are given by: cp ¼

X ri i

kf ¼

X i

r

cp;i ¼

yi ki ¼

X Mi ci r

i

cp;i

X ci P ki cj i

(13)

(14)

j

where the species properties ki, cp,i, and mi are given as a function of temperature [24].

3.2.

Gas diffusion layer (GDL)

As mentioned earlier, the gas diffusion layer uses an electrically conductive porous structure and being porous allows the convective and diffusive migration of the various species. The various transport quantities are volume averaged over a representative elementary volume, which contains many pores (Kaviany [25]). For the sake of model simplicity, we used the same notation as in the plain medium, however any changes due to the presence of the solid matrix is taken into account in the properties as will be discussed below. All the source terms in both the anode and cathode GDL are nil with the exception of the darcean terms of the momentum Eqs. (2) and (3) and the energy Eq. (4) with the heat source ðShs Þ which represents the heat transfer from GDL to bi-polar plates. The Darcy source term accounts for the pressure drop caused by frictional drag in the pore structure and is directly proportional to the averaged local velocity. The species and energy conservation equations keep the same form as those given in the gas channel domain but with averaged variables (velocities, temperature and concentrations) and corrected fluid property (diffusion coefficient) to account for the porous structure. The species effective diffusion coefficients are given by [5]: Di;mix;eff ¼ Di;mix 31:5

(15)

model assumes that no phase change for water occurs (monophase model). The governing equations in the CL layers are similar to those applied in the GDL except for activated source terms to account for chemical reactions. The source terms for the continuity Eq. (1) are only activated in the catalyst layers and are given by SH2 ¼ ja =2F; SO2 ¼ jc =4F; SH2 O ¼ jc =2F. These expressions model the consumption of hydrogen and oxygen respectively in the anode and cathode catalyst layers and production of water vapor in the cathode catalyst layer. The anode and cathode local current densities per unit volume, ja and jc, are expressed in terms of the Butler–Volmer equation for the anode is as follows: !gH2   ceff aa ha ð2FÞ H2 sinh (17) ja ¼ 2Acv i0a;ref cH2;ref RT which is used for modeling the reaction 2H2 / 4Hþ þ 4e and for the cathode: !gO2   ceff ac hc ð4FÞ O2 sinh jc ¼ 2Acv i0c;ref (18) cO2;ref RT which is used to model the reaction 4Hþ þ 4e þ O2 / 2H2O. Acv (m1) is the specific surface reaction area of the catalyst eff layer. ceff O2 and cH2 are respectively the effective oxygen and hydrogen concentrations at the surface of the agglomerate while cH2;ref and cO2;ref are H2 and O2 reference concentrations. In order to model, the concentration of the reactants at the surface of the catalyst particle and incorporate the microscopic mass transfer, the local mass diffusion equation is solved for half a sphere embedded in the polymer as shown in Fig. 2. This allows the calculation of a fictitious surface concentration on the catalyst surface. A representation of the domain around the platinum particle is depicted in Fig. 2. The mass transfer equation in spherical coordinates is given by   d 2 dC r ¼0 dr dr

(19)

where 3 is the porosity and Di,mix is given by Eq. (10). The effective thermal conductivity of the porous medium is given by: keff ¼ 3kf þ ð1  3Þks

(16)

where kf is the conductivity of the fluid mixture given by Eq. (14) and ks is the solid conductivity. The source term Shs in the energy equation represents the heat removal from the system and is activated on the channels/GDL interfaces at y ¼ hc and H  ha. This source term is given by (per unit depth) Shs ¼ kal =hðT  Tamb Þ; h ¼ ha or hc where kal is the thermal conductivity of aluminum.

3.3.

Catalyst layer (CL)

The catalyst layer can be described as a union of three different parts: the porous matrix from the GDL, the polymer and the catalyst particles. Therefore, the catalyst layer is modeled as a volume region and not simply as an interface between GDL and membrane. This description gives the possibility of coexistence of gas species, dissolved water and protons in the catalyst layer region. As mentioned earlier this

GAS (Pore) H+

Polymer

R0

H

ee-

Carbon

O RP H

Pt O

O

Fig. 2 – Microscopic model of mass diffusion around platinum particles.

3096

international journal of hydrogen energy 34 (2009) 3091–3103

with the following boundary condition: CðR0 Þ ¼ CB

(20)

for the bulk concentration away from the reaction surface and represents the volume averaged values. Note that the bulk concentration CB is calculated from the equation of species (Eq. 5) for each reacting species (i.e., H2 and O2). On the reaction surface, the rate of species diffusion is equal to the reaction rate, which is proportional to the surface concentration:  vC D  vr

(21)

¼ aCðRP Þ

r¼RP

where a represents the reaction rate at the surface and is derived from the Butler–Volmer equation assuming gO2 or H2 ¼ 1   i0c;ref ac hc ð4FÞ sinh a¼2 cO2;ref RT

C1 þ C2 r

(22)

Using the boundary conditions (20) and (21) allow the computation of the constants C1 and C2 and Eq. (22) can be used to calculate the fictitious surface concentration and is given by C ¼ Cðr ¼ RP Þ ¼ CB

DRRP0

!

aðR0  RP Þ þ DRRP0

(23)

Eq. (23) is modified using a correction factor to account for the fact that the model introduced above for the region around the platinum particle is ideal and one dimensional and does not take into account the effect of neighboring particles. Then the surface concentration in Eq. (23) becomes: 

C ¼ xC

DRR0P

B

!

aðR0  RP Þ þ DRRP0

(24)

where x is determined experimentally from the current limit. The Butler–Volmer equation is modified with the fictitious surface concentration ic ¼ Acv hdc TaCO2

and ia ¼ Acv hda TaCH2

The coefficients hdc and hda are included to account for the fact that the oxygen and hydrogen species are dissolved in the polymer and the dissolved concentrations are given by Henry’s law relation [14]: cdH2 ¼ ha TcH2 and cdO2 ¼ hc TcO2 . ha in Eq. (17) is the anode activation overpotential given by: ha ¼ fe;a  fm

(27)

with 3m,c the membrane fraction in the catalyst layer. The dissolved water equation takes into account the transport of water within the membrane due to diffusion and osmotic drag caused by the migration of Hþ through the membrane. The osmotic drag is given by the source term Sdrag ¼ Vðndrag ðsm / V fm =FÞÞ [14] where ndrag is the measured drag coefficient: lH2 O=SO3 22

(28)

where lH2 O=SO3 represents the number of water moles per mole of sulfonic acid and is given by lH2 O=SO3 ¼ ðEW=rdry ÞðcDW Þ. The diffusion coefficient of the dissolved water in the membrane is given by DDW ¼ ð1:3  104 ð expð2416ðð1=303Þ  ð1=TÞÞÞ. The dissolved water Eq. (6) takes into account the conversion of dissolved water into water vapor occurring in both catalyst layers though the source term Sd given by Sd ¼ hd(hw  T  Cdw  Cwv).

3.4.

The membrane (PEM)

The membrane must successfully transfer protons from the anode side to the cathode side and block the transfer of gas species H2, O2,N2 and H2O. Therefore in the membrane, the permeability for the Darcy term is set to be a very small number (K ¼ 1030 m) in order to inhibit the migration of the gas species due to pressure differences between the cathode and the anode. Since the approach chosen here is single domain, the momentum equations are solved in the membrane and they result into negligible velocity across it. The equations that describe transport in the membrane are Eq. (7) for protonic Hþ transfer (protonic current conservation), conservation of dissolved water (CDW) Eq. (6) and the energy Eq. (4). The protonic source term in the membrane is nil since no chemical reaction is occurring. The average fuel cell current density is calculated using the local protonic current passing through the membrane and is calculated using ohm’s law at the CL/membrane interface: Z 1 L vfm dx (29) I¼ sm vy L 0

(25)

where fe;a is the anode electrode potential given by fe;a ¼ Vth  Vcell , and fm the local membrane potential. In Eq. (18), hc is the cathode activation overpotential given by: hc ¼ fe;c  fm

seff ¼ sm 3m;c

ndrag ¼ 2:5

Solving Eq. (19) gives the local concentration as a function of the radius CðrÞ ¼

For the energy equation the source terms activated in the catalyst layers are the source terms due to the reversible heat and ohmic resistance heating given respectively by: Srev ¼ ðjd =nd FÞðTDSentropy;d Þ (d ¼ a,c) and Sohm;m ¼ ðsm / Vfm Þ / V fm [15] where sm is the local protonic conductivity given by sm ¼ ð0:5139l  0:326Þ expð1268ðð1=303Þ  ð1=TÞÞÞ [8]. As for the local properties, the protonic conductivity in the catalyst layer is corrected with membrane fraction:

(26)

Similarly, fe;c is the cathode electrode potential set to zero.

3.5.

Boundary conditions

The applied boundary conditions are as follows: x-direction  For the cathode inlet (x ¼ L, 0  y  hc) a parabolic profile for u with the average cathode velocity uc specified and

international journal of hydrogen energy 34 (2009) 3091–3103

v¼0 CH2 ¼ 0; CO2 ¼ C1 ; CN2 ¼ C2 ; Cwv ¼ C3 ; CDW ¼ 0 vfm ¼0 vx T ¼ Tc  for the CL and membrane side boundary conditions at x ¼ L, hc < y < H  ha u¼v¼0

vCi vfm vT ¼ ¼ ¼0 vx vx vx  for the anode inlet (x ¼ L, H  ha  y  H ), a parabolic profile for u where the average anode velocity ua is specified

3097

In order to perform a sensitivity analysis of the grid used in this study, two grids with different mesh densities are used to develop the polarization curve over the entire current density range for the parameters of the validation case discussed below. For the first grid the number of grid points used the xand y-directions are 22 by 68 respectively. While for the second grid the number of grids is nearly doubled in both directions giving a total number of grid points of 42 by 121. The results for the current density given by the two grids changed by an average of 2% over the entire range of the polarization curve. Therefore, we deemed that the accuracy given by the coarser grid is acceptable and we used it to simulate all cases except for the ones where geometries change (e.g., changing membrane thickness) and a different grid has to be constructed in order to accommodate the new boundaries. For the coarse grid, in the GDL layers we used 9 grids for each one (average grid size of w0.02 mm) while for the CL and membrane we used 6 grids for each (average grid size of w0.02 mm).

v¼0 CH2 ¼ C4 ; CO2 ¼ 0; CN2 ¼ 0; Cwv ¼ C5 ; CDW ¼ 0 and

vfm vx

¼0 T ¼ Ta  For the exhaust boundary conditions at x ¼ 0, 0 < y < hc and H  ha < y < H vu vv ¼ ¼0 vx vx

vCi vfm vT ¼ ¼ ¼0 vx vx vx  For the CL and membrane side at x ¼ 0, hc < y < H  ha vCi vfm vT ¼ ¼ u ¼ v ¼ 0 and ¼0 vx vx vx y-direction: At the lower and upper boundaries of the computational domain, i.e. y ¼ 0 and y ¼ H, respectively, the following boundary conditions are applied: u ¼ v ¼ 0 and

vCi vT vfm ¼ ¼0 ¼ vy vy vy

fe;a ¼ Vth  Vcell

and fe;c ¼ 0

3.6.

The numerical approach

The set of differential equations governing the PEMFC model along with the boundary conditions given above are solved using the finite volume method based on the SIMPLE algorithm (Patankar [26]). The equations are solved using a nonuniform staggered grid. Since for a two-dimensional model, the current density varies across the membrane, in this solution the fuel cell potential is specified and the average current density is calculated using Eq. (29).

4.

Results and discussion

4.1.

Model validation

The numerical results for the cell voltage as a function of current density are compared to the experimental and numerical results reported by Siegel et al. [14] for the operating conditions given in Table 1. The model results for the polarization curve are shown in Fig. 3 along with the experimental data. The model results show good agreement with the experimental data over the entire range of fuel cell voltage and the average error is within 7%. For the low current density range the present model, is very accurate and predicts the experimental results within 2.5% and in the high current density range the present model is within 10%. The errors are due to the fact that the transport in PEFMC involves many nonlinear and linked equations and properties which introduce errors especially at high current densities when the gradients are pronounced. At higher current densities, the present model shows as in the experiments and the model of Siegel et al. that the fuel cell reaches a current limit. This current limit effect is mainly caused by the concentration overpotential. This occurs as the rate of the reaction becomes much too large, compared to the speed of species diffusion, which limits the amount of reacting species reaching the reaction sites. Compared to the numerical results of Siegel et al. [14], both models predict the current limit. Their model seems to be more accurate in the transition regime before reaching the current limit. However, the present model seems to work better at low current and high voltage. The effect of the current limit due to species diffusion is clearly shown in Fig. 4, which shows the overpotential losses as a function of the current density. At higher currents (i.e., higher Vth  Vcell) and as the current limit is reached, the ohmic and activation overpotentials becomes constant and the increase in losses is solely due to species mass diffusion speed from the bulk (or void spaces) to the reaction sites. At lower currents, the dominant effect becomes that of the activation overpotential since the current density is too low in

3098

international journal of hydrogen energy 34 (2009) 3091–3103

Table 1 – Properties validation data. Property Gas channel length (mm) Gas channel width (mm) Gas channel height (mm) Gas diffusion layer thickness (mm) Catalyst layer thickness (mm) Membrane thickness (mm) Collector thickness (mm) Gas diffusion layer void fraction Catalyst layer void fraction Membrane volume fraction in the catalyst layer Membrane void fraction Faraday’s constant, F Specific porous area of the catalyst layer (mm1) Reference anode exchange current density (A/mm2) Reference cathode exchange current density (A/mm2) Anodic transfer coefficient Cathodic transfer coefficient Oxygen reference concentration (mol/m3) GDL thermal conductivity W m1 s1 Membrane thermal conductivity W m1 s1 Correction factor z Catalyst particle diameter RP (mm) Characteristic polymer radius R0 (mm)

Values

Sources

66

[14]

2

[14]

2

[14]

0.175

[14]

0.1

[14]

0.127 3 0.5

[14] [14] [14]

0.09

[14]

0.72

[14]

0.28 96,485 8000

[14] [14] [14]

5  106

[14]

6  1010

[14]

1 1 4.55

[27] [27] [27]

150.6

[8]

100

[8]

17 0.001

[14] [14]

0.01

[14]

Fig. 3 – Comparisons of numerical model and experimental results [14].

modeling the CL as a distinct region rather than an interface using two- or three-dimensional CFD model. Fig. 6 shows the oxygen concentration profile in the cathode gas channel, GDL and catalyst layer. The oxygen concentration decreases in the flow direction as it is being consumed. The effect of oxygen consumption can be clearly seen in the catalyst layer region where a significant concentration gradient exists. As in the case of hydrogen, the concentration gradient is highest in the reaction zone where oxygen is combining with the protons and the electrons to form water vapor. Fig. 7 shows the water vapor concentration in the cathode. Water vapor concentration gradient is highest in the cathode reaction zone where it is created by the chemical reaction and dragged by osmotic effects from the anode. From the anode

order to have any significant ohmic or species diffusion losses. For the case presented in Fig. 4, the membrane is well hydrated, therefore the ohmic losses are smaller than the activation and concentration overpotential.

4.2.

Results discussion

All the results discussed below are for a cell voltage of 0.6 V, since at this level of current density and voltage there are more pronounced gas species concentration gradients due to the higher species consumption or production. Fig. 5 shows the hydrogen concentration profile in the anode gas channel, GDL and CL regions. The hydrogen concentration starts from the inlet value and decreases through the gas channel, GDL and CL as hydrogen is consumed. The concentration gradient in the catalyst layer is highest due to the local reaction occurring in this region. The hydrogen distribution also shows that the gradient in the hydrogen concentration is significant in both directions. Hence, this shows the significance of

Fig. 4 – PEMFC potential losses.

3099

international journal of hydrogen energy 34 (2009) 3091–3103

CL GDL

Channel

U

Channel

U

GDL CL

Fig. 5 – Hydrogen concentration in the anode gas channel, GDL and CL.

CL, the water vapor diffuses to the cathode GDL then to the channel where it is being convected out of the computational domain. The highest water vapor concentration is observed in the CL near the inlet region. This is due to the fact that the reaction rate is highest in this region; therefore, producing more water vapor. The humidity content in the inlet gases have a significant effect on the fuel cell efficiency since the polymer protonic conductivity is strongly dependent on water content. The effect of the cathode inlet gas humidity is shown in Fig. 8. The cathode water content does not seem to affect the results greatly due to the fact that in this region water is produced through the chemical reaction. In addition water is transported from the anode to the cathode by osmotic drag and diffusion of the dissolved water. The effect of humidity content at the anode inlet is shown in Fig. 9. The results show a significant drop in current density as moisture content in the anode inlet decreases. This is due to the fact that the polymer in the catalyst layer and membrane contains less moisture; therefore, the protonic conductivity is lower which increases the ohmic loss. The dry condition in the anode is exacerbated by osmotic drag of water molecules from the anode to the cathode and diffusion of dissolved water from the cathode to anode is not enough to humidify the polymer on the anode side. Fig. 10 shows the effect of the cathode inlet pressure on fuel cell performance. As the cathode pressure increases with oxygen concentration, the overall cell performance increases since the concentration overpotential effects become smaller.

Fig. 6 – Oxygen concentration in the cathode gas channel, GDL and CL.

This is due to the fact that local current production is linearly dependent on the local concentration as given by the Butler– Volmer Eq. (18). Also the current limit of the fuel cell increases as the concentration of oxygen increases. For low inlet pressure, the drop in cell performance is drastic as shown for the case of 1 bar inlet pressure. The current limit of the fuel cell is reached even at high cell voltage when the inlet pressure is low; in this case a low current limit of 0.12 A/cm2 is reached at 0.84 V. The performance of PEM fuel cell can be drastically improved by increasing the inlet pressure but in practical applications, this increase has a limit since the membrane cannot sustain a large pressure difference between the anode and cathode sides and may result in internal leaks of the fuel and the oxidant. In addition higher pressure applications require much better sealing. The effect of the catalyst layer porosity on the fuel cell performance is shown in Fig. 11. In this case the fraction of carbon/platinum amount is kept as the same while changing the polymer volume and the void fractions. In all of the cases shown in Fig. 11, the combination of polymer volume and void fraction must be 0.81 (¼ 1  0.19 where 0.19 is the carbon/ platinum volume fraction). The results in Fig. 11 show that a very low catalyst layer porosity can have a significant impact on the fuel cell performance by reducing the overall current limit. A low porosity in the CL layer will increase the resistance to species transfer and therefore increasing the concentration overpotential.

3100

international journal of hydrogen energy 34 (2009) 3091–3103

CL GDL

Channel

U

Fig. 9 – Effect of anode RH on the current density (Cathode RH [ 100%).

Fig. 7 – Water vapor concentration in the cathode.

The effect of GDL porosity is shown in Fig. 12. The results show that as the porosity of the GDL decreases, the cell performance decreases as well as the current limit. As the GDL porosity decreases, the void space used for mass diffusion becomes tighter; therefore, the concentration overpotential becomes significant. As a matter of fact these results show two types of concentration overpotential, namely macroscopic and microscopic overpotentials. The former is caused at the GDL layer (macro effects) and results in the reduction of the current limit level. The latter is caused at the agglomerate side (micro effects) as modeled in the Butler–Volmer equation and this induces the current limit. As mentioned earlier, this

Fig. 8 – Effect of cathode RH on the current density (Anode RH [ 100%).

model assumes that the variation of the electronic potential in the CL and GDL is negligible. This assumption becomes no longer true as the porosity of the GDL increases and the decreased amount of solid will increase the electrical resistance. Therefore for a given operating condition, there will be an optimum GDL porosity for which the species mass transfer and electron flow combine to result into the highest possible current density. The most important component in the operation of the fuel cell is the membrane; therefore, in this section the thickness of the membrane is investigated. Fig. 13 shows the performance curve of the PEMFC for different values of membrane thickness. As the thickness of the membrane decreases, the ohmic losses due to the protonic conductivity decreases and for the same current density a higher cell voltage can be achieved. This is shown by the power delivered by the fuel cell. The results indicate that by reducing the membrane

Fig. 10 – Effect of cathode inlet pressure.

international journal of hydrogen energy 34 (2009) 3091–3103

3101

Fig. 13 – Effect of membrane thickness on cell voltage.

Fig. 11 – Effect of catalyst layer porosity on PEMFC performance.

thickness from 0.1 mm to 0.01, the power density can be doubled. Another advantage that the smaller membrane thickness offers the reduced heat generation due to lower ohmic losses and thereby reducing the need for more cooling which incurs higher costs and occupies more volume due to larger cooling ducts, heatsink fins, etc. However, the minimum membrane thickness that can be used is dictated by the structural integrity of the membrane, manufacturing processes, and cross diffusion of the species. For a membrane thickness of 0.01 mm, the estimated current loss due to cross diffusion, using one-dimensional analysis, is about 0.013 A/

Fig. 12 – Effect of GDL porosity on cell performance.

cm2 and it is 0.0013 A/cm2 for a membrane thickness of 0.1 mm. For the first case, the current loss corresponds to about 1% loss in power when the cell operates at the current limit. This percentage increases as cell current decreases and the voltage increases since current cross diffusion depends on species concentration on either side of the membrane. Fig. 14 shows the effect of the catalyst layer thickness on the current density. As the catalyst layer thickness decreases, the current density also decreases. This effect is mainly due to lower number of reaction sites (Pt particles). For higher CL thickness (>1 mm), the catalyst layer does not have any effect on the polarization curve since the ohmic loss in the membrane becomes the dominant effect and the concentration overpotential is comparatively small. It is well known that the PEMFC performance depends greatly on the cell temperature. The results of the present model confirm this dependence as shown in Fig. 15 where the temperature is varied between 328 K and 368 K while maintaining 100% RH for both the anode and cathode. The water vapor concentration for 100% RH for 318 K is 5.53 mol/m3 and 23.6 mol/m3 for 368 K. The results in Fig. 15 show that as the difference in cell performance can be cut by 25% when the

Fig. 14 – Effect of catalyst layer thickness on cell voltage.

3102

international journal of hydrogen energy 34 (2009) 3091–3103

Fig. 15 – Effect of temperature on cell performance.

temperature decreases from 348 K to 318 K while maintaining a fully humidified inlet gases. Also, the current density limit decreases as the temperature decreases since less reacting species are dissolved in the membrane, therefore affecting the reaction rate by increasing the concentration overpotential. The drop in power when the temperature changes from 368 K to 348 K, is about 15% and when it changes from 368 K to 318 K, the drop in overall power density is about 25%. Also, at higher operating temperature, the cell voltage as a function of current density exhibits a lower slope of variation until it reaches the current limit. This shows that for higher temperatures the fuel cell performance offers an improved stability due to the smaller variation of cell voltage as the current density increases. However in order to achieve the beneficial effects of the higher operating temperature, the humidification of inlet gases on molar basis has to be increased in order to attain a higher hydration level of the membrane. If the water vapor content on molar basis is kept the same as the temperature increased, the cell performance will only improve marginally. The model was also run for a smaller size channel, 1 mm high by 1 mm wide, on both the anode and cathode sides. The results show that the current density remained the same except near the current limit where the current limit decreased by 1.5%. Further decrease of the channel height below 1 mm, is not beneficial since the pumping power will increase drastically especially if the flow channel length is much longer than the length used in the present study. In order to address fully the effect of channel height for a real application, the entire channel length with serpentine flow pattern needs to be examined with a full three-dimensional model.

a role in determining its performance, namely the catalyst layers, the membrane, the GDL, and the channels. The model solves the mass, momentum, energy, gaseous species, potential and the electrochemical equations. The model uses variables properties as function of the local temperature and species. The model includes the effect of mass transfer around the catalyst particles by including a fictitious local concentration, which represents the species concentration on the surface of the particle. The model results for current density are compared to experimental results available in the literature and they show good agreement. The model is then used to conduct a parametric study to show the effect of the main parameters on the performance of the fuel cell. One of the most significant results is that the performance of the PEMFC is affected by the moisture content on the anode side much more than the cathode side. Therefore it is of great importance to maintain a fully humidified hydrogen gas in order to achieve high fuel cell performance. The model shows that there are two concentration overpotential effects, which are the microscopic and macroscopic effects. The first one is due to slow diffusion effect at the agglomerate level. The second one, which becomes important as the porosity of the GDL or CL become smaller, is due to the smaller void space necessary for mass transfer within the pores. This model shows that the optimal catalyst layer thickness is about 0.1 mm as for the membrane thickness, the smaller dimension offers less ohmic resistance however manufacturing limits, reliability, and cross diffusion of species dictate the minimum possible membrane thickness that can be used. Based on the present two-dimensional model, the optimum parameters needed to achieve high performance are as follows:  Catalyst layer porosity: 0.1  Channel height: 1 mm  Membrane thickness: 0.1 mm The optimum GDL thickness and porosity cannot be found with the present model since it does not solve for the electronic potential in the porous solid. In order to conduct a comprehensive optimization, further work is needed to improve this model such as the inclusion of liquid water management and the solution of the electronic potential (electron transfer) in the GDL and CL. This model will also be extended to simulate three-dimensional effects in order to capture geometrical phenomena that are not possible to simulate using two-dimensional models and these effects include interdigitated flow channel, point contact of the bipolar plates with the GDL, heat transfer and cooling the fuel cell from the exterior.

references

5.

Conclusions

In this study, a two-dimensional numerical model for a PEMFC is developed in order to enable optimization of fuel cell design. This model includes all the regions of the fuel cell that play

[1] Cheddie D, Munroe N. Review and comparison of approaches to proton exchange membrane fuel cell modeling. Journal of Power Sources 2005;147:72–84.

international journal of hydrogen energy 34 (2009) 3091–3103

[2] Standaert F, Hemmes K, Woudstra N. Analytical fuel cell modeling. Journal of Power Sources 1996;63:212–34. [3] Standaert F, Hemmes K, Woudstra N. Analytical fuel cell modeling; non-isothermal fuel cells. Journal of Power Sources 1998;70:181–9. [4] Charradi H, Halouani K. PEMFC modeling, graduation project report. Tunisia: ENIS, University of Sfax; 2005. [5] Bernardi DM, Verbrugge MW. Mathematical model of a gas diffusion electrode bonded to a polymer electrolyte. AIChE Journal 1991;37(8):1151–63. [6] You L, Liu H. A parametric study of the cathode catalyst layer of PEM fuel cells using a pseudo-homogeneous model. International Journal of Hydrogen Energy 2001;26:991–9. [7] Fuller TF, Newman J. Water and thermal management in solid polymer electrolyte fuel cells. Journal of the Electrochemical Society 1993;140(5):1218–25. [8] Gurau V, Liu H, Kakac¸ D. Two dimensional model for proton exchange membrane fuel cells. AIChE Journal 1998;44:2410–22. [9] Eaton BM. One dimensional, transient model of heat, mass, and charge transfer in a proton exchange membrane, MS thesis. Virginia Tech; 2001. [10] Genevey DB. Transient model of heat, mass and charge transfer as well as electrochemistry in the cathode catalyst layer of a PEMFC, Ph.D. thesis. Virginia Tech; 2001. [11] Nguyen PT, Berning T, Djilali N. Computational model of a PEM fuel cell with serpentine gas flow channels. Journal of Power Sources 2004;130:149–57. [12] Ramousse J, Deseure J, Lottin O, Didierjean S, Maillet D. Modeling of heat, mass, and charge transfer in a PEMFC single cell. Journal of Power Sources 2005;145:416–27. [13] Hmidi R, Halouani K. Polymer electrolyte fuel cells (PEMFC) technology & mathematical modeling, graduation project report. Tunisia: EPT, University of Carthage: 2004. [14] Siegel NP, Ellis MW, Nelson DJ, Von Spakovsky MR. Single domain PEMFC model based on agglomerate catalyst geometry. Journal of Power Sources 2003;115:81–9. [15] Siegel NP, Ellis MW, Nelson DJ, Von Spakovsky MR. A twodimensional computational model of a PEMFC with liquid water transport. Journal of Power Sources 2004;128:173–84.

3103

[16] Sivertsen BR, Djilali N. CFD-based modeling of proton exchange membrane fuel cells. Journal of Power Sources 2005;141:65–78. [17] Um S, Wang C-Y. Computational study of water transport in proton exchange membrane fuel cells. Journal of Power Sources 2006;156:211–23. [18] Wu J, Liu Q, Fang H. Toward the optimization of operating conditions for hydrogen polymer electrolyte fuel cells. Journal of Power Sources 2006;156:388–99. [19] Liu HC, Yan WM, Soong CY, Chen F. Effects of baffle-blocked flow channel on reactant transport and cell performance of a proton exchange membrane fuel cell. Journal of Power Sources 2005;142:125–33. [20] Liu HC, Yan WM, Soong CY, Chen F, Chu HS. Reactant gas transport and cell performance of proton exchange membrane fuel cells with tapered flow field design. Journal of Power Sources 2006;158:78–87. [21] Yan WM, Liu HC, Soong CY, Chen F, Cheng CH. Numerical study on cell performance and local transport phenomena of PEM fuel cells with novel flow field designs. Journal of Power Sources 2006;161:907–19. [22] Harvey D, Pharoah JG, Karan K. A comparison of different approaches to modeling the PEMFC catalyst layer. Journal of Power Sources 2008;179:209–19. [23] Kamarajugadda S, Mazumder S. On the implementation of membrane models in computational fluid dynamics calculations of polymer electrolyte membrane fuel cells. Computers and Chemical Engineering 2008;32: 1650–60. [24] Todd B, Young JB. Thermodynamic and transport properties of gases for use in solid oxide fuel cell modeling. Journal of Power Sources 2002;110:186–200. [25] Kaviany M. Principles of heat transfer in porous media. 2nd ed. New York: Springer; 1999. [26] Patankar SV. Numerical heat transfer and fluid flow. NY: Hemisphere Publication, Mc Graw-Hill; 1980. [27] Guvelioglu GH, Stenger HG. Computational fluid dynamics modeling of polymer electrolyte membrane fuel cells. Journal of Power Sources 2005;147:95–106.