A Three-Dimensional, Multicomponent, Two-Phase Model for a Proton

738

Energy & Fuels 2006, 20, 738-747

A Three-Dimensional, Multicomponent, Two-Phase Model for a Proton Exchange Membrane Fuel Cell with Straight Channels Guilin Hu†,‡ and Jianren Fan*,† State Key Laboratory of Clear Energy Utilization, Zhejiang UniVersity, Hangzhou, 310027, P. R. China, and School of Light Industry, Zhejiang UniVersity of Science and Technology, Hangzhou, 310012, P. R. China ReceiVed August 10, 2005. ReVised Manuscript ReceiVed December 16, 2005

A three-dimensional, steady-state, nonisothermal gas-liquid two-phase mathematical model for the proton exchange membrane fuel cell (PEMFC) is developed to better investigate transport phenomena and electrochemical kinetics in the whole fuel cell with straight channel and to predict the performance of fuel cell. Almost all transport phenomena occurring in the fuel cell such as fluid flows, heat transfer, mass transfer, electrochemical kinetics, and charge transport are accounted for in this comprehensive model. Using a multiphase mixture model to describe two-phase flow and transport phenomena, this model can handle the situation where a single-phase region coexists with a two-phase zone in PEMFC easily. An important feature of this model is that it can model and predict electrochemical kinetics in detail (i.e., transport of electrons in the carbon phase), and protons in the membrane phase are accounted for spatially. Three-dimensional spatial distributions of temperature, species concentrations, and carbon-phase and membrane-phase potentials are illustrated and discussed in detail by numerical simulation. Furthermore, comparison of predicted performance and experimental data published in the literature is conducted to verify the mathematical model. Last, effects of velocity and humidification temperature of the cathode inlet stream and shoulder-to-channel ratio on the cell performance have been analyzed.

1. Introduction Proton exchange membrane fuel cells (PEMFCs), or polymer electrolyte fuel cells (PEFCs), are a type of lower-temperature fuel cells and have emerged as strong candidates as future power sources for both stationary and automotive applications. Especially for automotive applications, this fuel cell is one of the most promising types of fuel cells to replace conventional power conversion devices such as internal combustion engines because of its many distinctive attributes such as high energy efficiency, low noise, low emission, and especially low operating temperature (i.e., quick startup under environmental temperature). However, its performance must be improved and costs reduced before this system becomes commercially viable. Compared to an experimental study, a numerical study can reduce the cycles of cell design greatly and get many important results that cannot be obtained by experiment or are too expensive for experiment. In the past decade, many computational models containing variable complex degrees have been reported in the literature.1-11 * Corresponding author. Fax: +86-0571-87991863. E-mail: fanjr@ zju.edu.cn. † State Key Laboratory of Clear Energy Utilization, Zhejiang University. ‡ School of Light Industry, Zhejiang University of Science and Technology. (1) Springer, T. E.; Zawodinski, T. A.; Gottesfeld, S. J. Electrochem. Soc. 1991, 138, 2334-2342. (2) Bernardi, D. M.; Verbrugge, M. W. J. Electrochem. Soc. 1992, 139, 2477-2491. (3) Eikerling, M.; Kharkats, Yu. I.; Kornyshev, A. A.; Volfkovich, Yu. M. J. Electrochem. Soc. 1998, 145, 2684-2699. (4) Fuller, T. F.; Newman, J. J. Electrochem. Soc. 1993, 140, 12181225. (5) Singh, D.; Lu, D. M.; Djilali, N. Int. J. Eng. Sci. 1999, 37, 431452.

According to the development process, the model is becoming more and more complex, and the researched phenomena are becoming more and more comprehensive. At the first stage, many one- and two-dimensional models were first presented, focusing on flow with transport of the reactants and products in the flow channels and across the membrane. The most prominent earlier work focuses on developing one-dimensional, isothermal models of the MEA,1-3 which is a kernel component of a PEMFC. The other one-dimensional models4,5 include a full hydrodynamic description of water and gas transport in the electrodes and the gas channels as well. The two-dimensional model can be categorized as along the channel or across the membrane. The latter is widely used to investigate transport phenomena in the fuel cell using interdigitated gas distributors.6,7 The pseudo two-dimensional models accounting for composition changes along the fluid channel were also developed by many research groups.8,9 Although the one- and two-dimensional models have presented good understanding on many physical phenomena and obtained many very interesting results that can be employed to optimize cell design and operation, the threedimensional model is required to investigate spatial distributions of characteristics parameters in detail. More recently, three(6) Kazim, A.; Liu, H. T.; Forges, P. J. Appl. Electrochem. 1999, 29, 1409-1416. (7) Yi, J. S.; Nguyen, T. V. J. Electrochem. Soc. 1999, 146, 38-45. (8) Nguyen, T. V.; White, R. E. J. Electrochem. Soc. 1993, 140, 21782196. (9) Gurau, V.; Liu, H. T.; Kakac, S. AIChE J. 1998, 44, 2410-2422. (10) Berning, T.; Lu, D. M.; Djilali, N. J. Power Sources 2002, 106, 284-294. (11) Dutta, S.; Shimpalee, S.; Van Zee, J. W. J. Appl. Electrochem. 2000, 30, 135-146.

10.1021/ef050254b CCC: $33.50 © 2006 American Chemical Society Published on Web 02/10/2006

Model for a Proton Exchange Membrane Fuel Cell

dimensional models10-13 based on the computational fluid dynamics approaches have been presented. In these models, not only the three-dimensional characteristics have been described, but also more transport phenomena, including heat and mass transfer, multicomponent, and even electron transport13 have been studied. Water management of PEMFC plays an important role for operation and optimization of cell performance of PEMFC. If too much liquid water exists in cell, mass transport resistance will be enhanced and some active reaction sites will be covered. On the other hand, when PEMFC operates at low humidity, the PEM will lose water, which causes a rapid increase of the ohmic resistance. Both these phenomena will decrease cell performance. The first mathematical model of PEMFC is a single-phase flow model, which cannot deal with the flooding phenomena, and occurred in the fuel cell especially when the cell operates at high current densities. Therefore, efforts to model two-phase flow and transport in the fuel cell are becoming a more critical issue. Yi and Nguyen7 developed one model, which considers liquid water, for interdigitated flow fields. However, this inconvenient model is based on mass conservation, and the liquid water flow induced by capillary forces and gas flow induced drag was addressed by semi-heuristic equations. Natarajan and Nguyen14 presented a 2D multiphase transport model of the cathode gas-diffusion layer. Later, this model was modified to a molar based equation by He et al.,15 and the equation of motion for liquid water was also modified to account for water movement by both capillary and gas shear induced force. More recently, Natarajan and Nguyen16 developed a pseudo-three-dimensional model by extending their twodimensional isothermal model.14 This model can be used to obtain qualitative insight into the distribution of liquid water in the backing layer and its effect on the distribution of the gasphase reactant species. Then, operating parameters such as temperature, stoichiometric flow rate, and inlet gas stream humidity were analyzed by numerical simulation. Berning and Djilali17 reported a 3D computational fluid dynamics multiphase model for a PEMFC. Transport of liquid water inside the gasdiffusion layers is modeled using viscous forces and capillary pressure terms. The physics of phase change is accounted for by prescribing local evaporation as a function of the undersaturation and liquid water concentration. To locate condensing regions in the electrode, Kermani et al.18 presented a threespecies condensing flow model for the cathode of a PEMFC. Shimpalee and Dutta19 developed a complete three-dimensional model to include the dimension parallel to the reactive surface to account for the electrode regions hidden from the gas channel. They treat the liquid water as a component of the gas mixture. Furthermore, the effect of liquid water and its distribution in the porous gas diffusion layer on the gas transport and electrochemical reactions were not accounted for. Modeling of (12) Hu, G. L.; Fan, J. R.; Chen, S.; Liu, Y. J.; Cen, K. F. J. Power Sources 2004, 136, 1-9. (13) Meng, H.; Wang, C. Y. J. Electrochem. Soc. 2004, 151, A358A367. (14) Natarajan, D.; Nguyen, T. V. J. Electrochem. Soc. 2001, 148, A1324-A1335. (15) He, W.; Yi, J. S.; Nguyen, T. V. AIChE J. 2000, 46, 2053-2064. (16) Natarajan, D.; Nguyen, T. V. J. Power Sources 2003, 115, 66-80. (17) Berning, T.; Djilali, N. J. Electrochem. Soc. 2003, 150, A1589A1598. (18) Kermani, M. J.; Stockie, J. M.; Gerber, A. G. Condensation in the Cathode of a PEM Fuel Cell. Proceeding of the 11th Annual Conference of the CFD Society of Canada, Vancouver, Canada, May 28-30, 2003; pp 292-299. (19) Shimpalee, S.; Dutta, S. Numer. Heat Transfer, Part A 2003, 8, 111-128.

Energy & Fuels, Vol. 20, No. 2, 2006 739

two-phase flow in the porous air cathode of a PEMFC was attempted by Wang et al.20 The transport coefficients in the stationary numerical two-phase model are parametrized as functions of the liquid water saturation. However, the half cell including the gas channel, gas diffuser, and catalyst layer in the cathode side of a PEMFC was considered in the article. Later, this model was further developed to include a whole cell by You and Liu.21 Recently, Pasaogullari and Wang22 used the multiphase mixture model to investigate two-phase transport in the porous layers of a PEMFC cathode, and effects of porosity, thickness, and wettability of a microporous layer on the two-phase transport in PEMFC are elucidated. In most of above model, the catalyst layer is accounted for as an interface between the diffusion layer and the membrane, which cannot really model the detailed special distribution of species concentration and activation overpotential in the catalyst layer. Furthermore, the temperature distribution was not considered in most articles. The modeling approach of our present work is based on the multiphase mixture model23 for analyzing multiphase flow and transport in the fuel cell. The major objective of this work is to develop a three-dimensional model for the simulation of PEMFCs, which can be applied to investigate the multidimensional distributions of current densities, reactants, and concentrations in the fuel cell. This comprehensive model will account for electrochemical and transport phenomena in the whole cell and the coupling between them. In particular, this model considers transport of the electron in the backing layer and catalyst layer, as well as transport of the proton in the membrane and catalyst layer. The effects of velocity and humidification temperature of cathode inlet stream and shoulder-to-channel ratio on the performance are discussed in this article. We consider the catalyst layer as volume not as an interface, and thus a detailed special distribution in the catalyst layer can be obtained, which is important for the design and optimization of the catalyst layer. 2. Description of the Mathematical Model Generally, a PEM fuel cell is composed of nine segments across the fuel cell; they are current collectors, gas channels, diffusion layers and catalyst layers on both anode and cathode sides and polymer membrane sandwiched between them, as illustrated in Figure 1. The two outer segments are graphiteplate current collectors, which are curved with flow channels for reactants and products to enter and exit the cell, respectively. The adjacent segments are electronically conducting porous backing layers that allow for even distribution of reactants to the anode and cathode. The centric polymer membrane spatially separates the fuel cell. It not only avoids species and electrons directly through the cell, but also provides channel for proton transport. The catalyst layers between the membrane and the backing layers are the domain for electrochemical reactions, which is considered as real width, not as simple as an interface as in many previous models. The catalyst layers provide electron, proton, and reactants channels. In this work, we assume that the fuel cell is repeated periodically along the y-direction. To (20) Wang, Z. H.; Wang, C. Y.; Chen, K. S. J. Power Sources 2001, 94, 40-50. (21) You, L.; Liu, H. Int. J. Heat Mass Transfer 2002, 45, 22772287. (22) Pasaogullari, U.; Wang, C. Y. Electrochim. Acta 2004, 49, 43594369. (23) Wang, C. Y.; Cheng, P. Int. J. Heat Mass Transfer 1996, 39, 36073618.

740 Energy & Fuels, Vol. 20, No. 2, 2006

Hu and Fan

sponding to the summary of every species source term, can be defined as

Sm )

∑k Sk

(2)

where k denotes chemical species that include hydrogen, oxygen, nitrogen, and water. Sk stands for the source/sink term of the kth species. 2.2. Momentum Conservation. The fluid flow in the channel, backing layer, and catalyst layer of the fuel cell can be described by the general Darcy equation as

Figure 1. Schematic structure of a PEMFC and computation domain. CH, channel; BP, bipolar plate; BL, backing layer; CL, catalyst layer; MEM, membrane.

save the computation source and computation time, we take one current collector, which is half of two flow channels around which is the computation domain and which represents one part of a three-dimensional whole PEMFC, shown in Figure 1. The following simplifying assumptions have been considered in this model development: (1) The flow in the fuel cell is laminar everywhere. This is reasonable for the low velocity and low Reynolds number. (2) Only the steady state is considered. (3) The proton membrane is impermeable to gas species. (4) Isothermal boundary conditions are used for the external wall and inlet stream. (5) The porous media such as backing layer, catalyst layer, and membrane are isotropic and homogeneous, characterized by effective permeability and uniform porosity. (6) Butler-Volmer kinetics governs the electrochemical reaction. (7) Electron conductivity in the current collectors is very big, that is, the potential in the current collector is constant. (8) The reactants dissolve in water before they are subjected to the electrochemical reaction. At the stage of our model, we will ignore this phenomenon and will assume that they are transported as components of the gas mixture. (9) Water exists as gas state in the fuel cell; condensation occurs after vapor partial pressure exceeds saturation vapor pressure. This is the weakest flaw in this nonisothermal model, which presumes that water is generated in the catalyst layer in the form of vapor, and will be improved in our future work. 2.1. Continuity Equation. In multiphase mixture model, the gas and liquid occupy the control volume together. The conservation equation of mass for the two-phase mixture (i.e., the continuity equation) can be written as

∇‚(Fu) ) Sm

∇‚(Fuu) ) -∇p + ∇‚(µeff∇u) + Su

where p denotes the pressure, µeff is the effective viscous coefficient, and Su is the momentum source. For the fluid flow in the gas channel, the equation becomes the common NavierStokes equation for zero momentum source term. For the fluid flow in the porous media, the momentum source Su is used to describe the Darcy’s drag force exerted by the fluid on the solid surface in the porous backing layer and catalyst layer

Su ) -

µeff 2 u K

(4)

where K is permeability of the porous media. In this case, the momentum equation can be simplified to the well-known Darcy’s law, because the inertia and viscous terms are so small compared with the source term and can be neglected. However, the macroscopic velocity varies rapidly near the interface between the channel and backing layer, and thus the macroscopic viscous force and inertial force cannot be neglected in the momentum equation. Therefore, a general equation with appropriate source terms is more suitable for describing the flow in the coupling channel and porous backing and catalyst layers. 2.3. Species Conservation Equations. The effective and rapid supply of reactants and emission of products have great impact on the performance of a fuel cell. The species conservation of a two-phase mixture can be defined as

∇‚(γkFuCk) ) ∇‚{

∑R [FRsRDk,R(∇Ck,R - ∇Ck)]} + ∇‚(FD∇Ck) - ∇(∑Ck,RJR) + Sk R

(5)

where R implies gas phase and liquid phase in our liquid-gas two-phase system, C stands for the mass fraction of species, D is the diffusion coefficient of species, s is phase saturation, and J is the interphase diffusion mass flux. γk is the convection correction factor of species and can be defined as

(1)

where is the porosity of the porous media, which is equal to the unit for the fluid channels. F is the density of mixture, and u is the intrinsic fluid velocity vector. u b reflects the superficial velocity in the porous media. The source term Sm denotes mass source/sinks to the multiphase mixture for each control volume in the computational domain. Because the mass increase of liquid water due to vapor condensation equals the mass decrease of water vapor in the gas phase, it has no effect on the final source term. Thus, this source term is zero in most of the computation domain except for the activated reaction sites as the anode and cathode catalyst layers. This term, corre-

(3)

γk )

FλlCk,l + λgCk,g FlslCk,l + FgsgCk,g

(6)

Mass source term Sk of a species induced by electrochemical reaction in the catalyst layer can be determined as

Sk )

r j nF

(7)

where n is number of transfer electron for the kth species, r is the stoichiometric number, and j is the transfer current density. Thus, the mass source term of reacting species


{


hydrogen, oxygen, and water reads as

ja - for H2 2F jc Sk ) for O2 4F jc - for H2O 2F

(8)

ηact ) φa,c - φm - Voc

For the water, the electroosmotic drag from the anode to the cathode and the electrochemical reaction must be accounted for. Thus, the mass source term of water in the cathode becomes

jc β SH2O ) ∇‚ - bı c F 2F

(

ζa and ζc stand for so-called transfer coefficients. The coefficient 1 - sl is used to consider the influence of liquid water on the area of catalyst reaction linearly. jo,ref is the reference transfer current depending on the operation temperature and catalyst loading and so on. ηact is the surface overpotential, defined as

)

(9)

where β is the electroosmotic drag coefficient and bıc stands for the local current density vector. Water vapor will condense when the vapor partial pressure is greater than the saturation vapor pressure corresponding to local temperature (i.e., the relative humidity of water vapor exceeds 100%). The relationship between saturation vapor pressure pv and temperature reads as:1

log10 pv ) 2.1794 + 0.02953T - 9.1837 × 10-5T2 + 1.4454 × 10-7T3 (10)

(17)

where φa,c and φm denote the potentials for the carbon phase and membrane phase, respectively. Voc is the reference opencircuit potential of the electrode. It equals zero at the anode but depends on temperature at the cathode.24 2.4. Potentials. The movements of electrons and protons are controlled by carbon-phase potential (solid part) and membranephase potential (polymer electrolyte), respectively. Generally, the backing layer only has electron transport (the proton can only be transported in the membrane), and the catalyst layer is an overlapping region of solid base and electrolyte membrane, and thus it can transport both electrons and protons. The electrical potentials of the membrane phase φm and carbon phase φa, φc are introduced following the concept of mean parameters. Taking local electroneutrality into account, we find that the electronic and ionic currents produced or consumed in the catalyst layers lead to a voltage drop via Ohm’s law according to:

i ) (σ∇φ

(18)

and the corresponding vapor density is

FHv 2O

)

pvMH2O

(11)

RT

where MH2O stands for the molecular weight, R is the general gas constant, and T is the temperature. Thus, concentration of water vapor in the gas mixture is determined as

CH2O,g )

FHv 2O Fg

(12)

where Fg is the density of gas mixture. Gas-phase saturation and liquid-phase saturation for a gasliquid system can be defined as

sl )

Fg(CH2O - CH2O,g) Fl(1 - CH2O) + Fg(CH2O - CH2O,g) sg ) 1 - sl

(13) (14)

Because species transport and electrochemical reaction are coupled with each other, the local equation of transfer current density distribution must be determined to obtain reactive species concentration distribution. In this work, the Butler-Volmer equation is introduced to describe transfer current density:

ja ) jref o,a

( )[ C H2

CHref2

jc ) -(1 -

γa

(

exp

sl)jref o,c

ζaa

RT

F‚ηact - exp - F‚ηact RT

( )[ C O2

COref2

γc

) ( )] ( ) ( )] ζca

The governing equation for each sort of current follows the current conservation:

∇i ) j

(19)

Thus, the governing equations of carbon-phase potential and membrane-phase potential become

∇‚(σe,eff∇φe) + Sφe ) 0

(20)

∇‚(σs,eff∇φs) + Sφs ) 0

(21)

where σe,eff is proton conductivity of electrolyte and σs,eff is electron conductivity of carbon phase of solid base. Sφe and Sφs are product and consumption rates of charge in the electrochemical reaction of catalyst layer, respectively, defined as:

j , At anode Sφe ) -ja, At cathode c

(22)

Sφs ) -ja At anode

(23)

Sφs ) jc

(24)

At cathode

The local membrane proton conductivity is a function of temperature and water content and is given by1

(15)

ζccF ζacF η - exp - ηact exp RT act RT

[ (3031 - T1)]

σe(T) ) (0.005139λ - 0.00326 ) × exp 1268

(25) and water content of the membrane surface is dependent

(16)

where subscripts a and c denote anode and cathode, respectively.

(24) Um, S.; Wang, C. Y.; Chen, K. S. J. Electrochem. Soc. 2000, 147, 4485-4493.


Hu and Fan

on water vapor activity

λ)

{

0.043 + 17.81a - 39.85a2 + 36.0a3, 0 < a e 1 14 + 1.4(a - 1), 1 < a (26)

where ie is the proton current density. However, in the catalyst layers, the energy equation includes heat production due to reversible and irreversible mechanisms, and therefore the heat source term can be determined as

ST )

where a is water activity at the membrane surface, which can be defined as

a)

X wp

(27)

psat

where Xw is the molar fraction of water vapor at the surface of membrane, p is gas pressure, and psat stands for saturation vapor pressure corresponding to local temperature. After obtaining the phase potentials, the operation current density in the external circuit I can be calculated by averaging local protonic current density at the interface between the membrane and catalyst layer:

I)

1 1 × × H L

∫0H ∫0L

( |) σe

∂φm ∂x

IF

dy dz

∇‚(λeff∇T) ) ∇‚(γhFCpuh) + ∇‚[

∑R (hRJR)] + ST

(29)

where λeff stands for an effective thermal conductivity; it is a combination of a multiphase mixture with the solid matrix phase. The second term on the right-hand side describes the energy flux due to relative phase motions, including both sensible and latent heat transport:

hl ) hg + hfg +

∫0T (cl - cg) dT

(30)

γh is the energy convection correction factor, defined as:

γh )

F[λlhl + λghg] Flslhl + Fgsghg

is2 σs,eff

where is is local current density. In the membrane, the corresponding heat effect reads as 2

ST )

ie σe,eff

(35)

Velocity: Fu ) Flul + Fgug

(36)

Concentration: FCk ) FlslCk,l + FgsgCk,g

(37)

Diffusion coefficient: FDk ) FlslDk,l + FgsgDk,g (38) Viscosity: µ )

F (KrlFl/µl) + (KrgFg/µg)

(33)

(39)

In a partially saturated porous media, the flow is strongly dependent on the local saturation, as effective diffusivity. In this work, the relative permeabilities for gas and liquid phases are represented by cubed phase saturation:21

KrR ) s3R

(40)

Relative mobilities are defined as

λR )

KrR/νR KrRFR/µR + KrRFR/µR

(41)

As in a traditional mixture theory, a diffusive mass flux of phase R within the two-phase mixture can be defined as

JR ) FRuR - λRFu

(42)

and

∑R JR ) 0 From eq 42, the individual phase velocities can be determined from the mixture flow field after phase fluxes are obtained. Liquid-phase mass flux Jl in the two-phase mixture can be determined as

Fλl Fλlλg (∇p - ∇pl) ) K ∇pc µ µ

Jl ) K

(32)

(34)

Density: F ) Flsl + Fgsg

(31)

ST is a volumetric heat source. In the backing layer, ohm heating by electron current is considered as:

ST )

)

It is notable that this equation is suitable for the current collector and membrane (the velocity is zero and becomes the pure heat conduction problem), fluid domain as fluid channel, and porous media such as the backing layer and catalyst layer. 2.6. Mixture Variables and Properties. The mixture variables and properties in the above formulation are defined as

(28)

where H is the height of the catalyst layer, L is the length, and IF implies the interface between the membrane and catalyst layer. Under the assumption of a perfectly conductive solid for current collectors, φ0a is equal to zero on the anode side at the anodic current collector and φ0c equals the cell voltage on the cathode side at the cathodic current collector. Solution of the problem for a given φ0c provides the mean current density in the fuel cell (i.e., a point on the I-V curve). 2.5. Conservation of Energy. Resolving the temperature distribution involves the current collector, channel and porous backing layer, catalyst layer, and membrane, and thus it is a coupled convection-conduction heat transfer problem. A generalized conservation equation for energy in the fuel cell is:

(

dU0 i2 + j ηact + T σeff dT

(43)

where pc is the capillary pressure between two phases. The capillary pressure is a function of liquid saturation, usually defined by experimental expression20 as

pc ) pg - pl ) ϑ

1/2

(K )

J(sl)

(44)

and JSl is the Leverett function:

J(sl) ) 1.417(1 - sl) - 2.120(1 - sl)2 + 1.263(1 - sl)3 (45)


2.7. Boundary Conditions. The periodic boundary conditions for all variables are imposed on the bottom and top ends (y ) 0 and y ) H): f|y)0 ) f|y)H, because the computational domain covers a part of the cell assembly repeating periodically along the y-direction. The von Neumann boundary condition is applied to membranephase potential for all boundaries of the computational domain equation (i.e., zero flux: ∂φm/∂n ) 0). The thermal boundary conditions of outer surfaces of both the cathode and anode are considered in the fixed temperature way, that is to say, the temperature on these surfaces equals the cell temperature. This can be employed by control flux of coolant in practice. At the inlets, the Dirac boundary conditions are adopted to all variables, including the fluid velocity, temperature and the species concentrations, and can be prescribed as

U ) Uin, V ) W ) 0, T ) Tk,in, Ck ) Ck,in At the outlets, the fully developed conditions are applied to velocity, temperature, and species concentration fields. At the interfaces between the gas channels and the plate collectors, boundary conditions of the first kind are prescribed for the gas mixture velocity component (no slip condition). 2.8. Numerical Procedures. All model equations containing mass, momentum, species, temperature, and charges were converted to a finite difference form by the control volume method. In this work, the employed approach for obtaining the polarization curve is the so-called potentiostatic, which is setting the cell potential and calculating current density. Thus, the numerical procedure is described as follows: first, we solved the anode carbon-phase potential equation, membrane-phase potential equation, and cathode carbon-phase potential equation in sequence, and we obtain convergence in inter-iteration. Then, the moment equation and hydrogen for the anode are iteratively solved for handling the variable parameters with composition variance. The solutions of the velocity component equations are obtained in a staggered control volume. The semi-implicit method for pressure-linked equations consistent (SIMPLEC) scheme was adopted for the pressure correction. In the same way, the moment equation and species for the cathode are solved. Last, the temperature of the whole cell is solved after the velocity and electrochemical rate are obtained. The source terms of the species concentration equation are dependent on the transfer current densities, which are relative to the phase potential. Thus, the transport equations for the gas components and phase potential are solved iteratively. All the coupled equations were solved iteratively until the relative error in each field achieved the specific convergent standard (usually 10-5). Once the converged solutions are obtained, the cell average current density is obtained as the arithmetic average of all local current density along the flow direction. Changing carbon-phase potential of the cathode current collector, different current densities, and polarization curves are obtained. 3. Results and Discussion The results presented below were obtained for the parameters and base conditions listed in Table 1, which are obtained mainly from Bernardi and Verbrugge.2 It is difficult to completely visualize the concentration, velocity, temperature, and so on in the fuel cell for the three-dimensional characteristics, and thus two-dimensional figures are adopted to describe these phenomena.

Energy & Fuels, Vol. 20, No. 2, 2006 743 Table 1. Physical Parameters and Base Conditions quantity

value

operating temperature/K length of the channel/cm channel height/cm shoulder height/cm channel width/cm backing layer width/cm catalyst layer width/cm membrane width/cm anodic transfer coefficient at anode, ζaa cathodic transfer coefficient at anode, ζca anodic transfer coefficient at cathode, ζac cathodic transfer coefficient at cathode, ζcc concentration parameter, γa concentration parameter, γc permeability of backing layer, K (cm2) porosity of the backing layer, 0 carbon-phase conductivity, σa,c (Ω-1‚cm-1) operating voltage/V pressure of anode inlet/kPa pressure of cathode inlet/kPa reference exchange current density at anode jo,ref/A‚cm-3 reference exchange current density at cathode jo,ref/A‚cm-3 velocity at the anode inlet/m‚s-1 velocity at the cathode inlet/m‚s-1 humidity temperature of anode inlet stream humidity temperature of cathode inlet stream

353.15 6.0 0.08 0.08 0.08 0.024 0.0028 0.025 0.5 0.5 2 2 0.5 1.0 1.76 × 10-7 0.5 40 0.6 3 × 101.3 5 × 101.3 1.0 × 103 1.35 × 10-4 0.36 0.6 60 °C 60 °C

Figure 2 shows concentration contours of hydrogen, oxygen, and water in the x-z plane at half-height channel (y ) 0.08 cm) and characterizes the distributions of species in the channel, backing, and catalyst layers. It can be seen from the figure that hydrogen and oxygen appear to follow the same trend (i.e., reduce along the channel), and water concentration increases along the channel. The reason for this phenomenon is that the electrochemical reaction occurring in the catalyst layer renders the reactants concentrations to decrease and products concentration to increase. Liquid water saturation is illustrated in Figure 3. The distribution is symmetric around the centerline in the z-direction for the same flow in the channel and increases along the channel for the electrochemical reaction in the catalyst layer. As a result of the electrochemical reaction and accumulation effect along with flow, liquid water may first appear in the vicinity of the interface between the membrane and cathode catalyst layer near the channel outlet. Water removal in front of the shoulders must be transported to the channel first (i.e., walks a longer path in comparison to the species in the region in front of the channel). Thus, the saturation in this region is greater than that in the region in front of the channel. Also, liquid water saturation in the channel is very small; comparatively, the value in the catalyst layer is much larger, which may bring side effects to the transport and reaction of species. Figure 4 plots representative 2D velocity vectors in the x-y plane at half-length of the fuel cell channel (z ) 3.0 cm). Figure 4a shows the velocity vectors in the cathode side; these velocities are pointing outward from the membrane to the channel. This implies that species have a bulk flow out from reaction sites, which is opposite to diffusion, because that oxygen is less consumed than water vapor produced in the catalyst layer of the cathode. Figure 4b displays the velocity field in the anode side. Unlike the cathode, the velocity points from the channel to the reaction surface, that is, species have a bulk flow from the channel to the reaction surface; this is a benefit to transport species. The reason is that hydrogen is consumed in the catalyst


Hu and Fan

Figure 3. Three-dimensional plots of the liquid saturation of water in the cathode (× 103).

Figure 4. Two-dimensional velocity fields in the x-y plane at halflength of the channel (z ) 3.0 cm). Figure 2. Species concentration distributions in the x-z plane at halfheight of the channel (y ) 0.08 cm).

layer, and water is transported from the anode to the cathode for electroosmotic drag. Figure 5 illustrates the distributions of oxygen concentration and local current density in the y-z plane in the middle of the catalyst layer. The concentration has maximal value lying in the middle of the inlet channel and decreases along the flow direction induced by the reaction consumption. The value in the region in front of the shoulder is smaller than that in front of the channel for the limitation of the diffusion, as showed in Figure 5a. The local current density is shown in Figure 5b. We can find that the maximal current density exists in the region in front of the shoulder, which is different from the results

obtained in the literature.10,20 This is because the spatial distribution of overpotential is considered in this model. The electrochemical reaction rate is influenced not only by reactant concentration but also by overpotential. It is highly nonlinear since the activation overpotential term goes into the exponential operator of the Butler-Volmer equations. The low concentration in front of the shoulder will decrease the current density, but the low overpotential causes larger current density in this area. That is, ohm potential loss has a more important impact on the current density than concentration loss under this case. It is also found that the distribution is similar to previous data when the fuel cell operates at high current density, because oxygen transport limitation becomes the main influence factor for current density distribution.



Figure 5. Distributions of oxygen concentration and local current density in the y-z plane in the middle of the catalyst layer.

Figure 6 shows the distributions of potential in the x-y plane at half-length channel (z ) 3.0 cm). Carbon-phase potential at the anode is plotted in Figure 6a. It can be seen that the change is very small (i.e., ohm loss caused by electron transport is very low with the order of l mV). This is because electron conductivity is very big in the carbon phase. Similar phenomena can be found in Figure 6c, which illustrates the distribution of carbon phase at the cathode. Furthermore, it should be noted that the gradient of potential has maximum value near the intersect position of channel and shoulder; this is because of the accumulation effect caused by the electron in the backing layer in front of the channel flowing to the shortest position. From Figure 6b, we can find that in the membrane the gradient of the protonic potential is constant due to the ohmic drop. Membrane potential has a large change from anode to cathode for the big membrane resistance, about 0.24 V. Thus, potential loss induced by membrane resistance occupies a large ratio in the total polarization. Figure 7 shows temperature distribution in the x-y plane at half-length channel (z ) 3.0 cm). It illustrates the nonuniform nature of the thermal field, with the higher temperatures

Figure 6. Potential distributions in the x-y plane at half-length of the channel (z ) 3.0 cm).

achieved inside the cell. Particularly, the maximal temperature appears in the catalyst layer of the cathode, because the electrochemical reaction occurring here can produce much heat as well as ohmic heating. Furthermore, temperature in the region in front of the channel is larger than that in the region in front of the shoulders. This is because that heat generated inside a cell is primarily removed through the backing layer to the current collector land by lateral conduction.25 The temperature change at the cathode is much greater than that at the anode, because much reactive heat is generated in the catalyst layer of cathode, and phase change also releases heat. (25) Ju, H.; Meng, H.; Wang, C. Y. Int. J. Heat Mass Transfer 2005, 48, 1303-1315.


Figure 7. Temperature distribution in the x-y plane at half-length of the channel (z ) 3.0 cm).

Hu and Fan

Figure 9. Effect of gas humidification temperature of the cathode on cell performance.

Figure 10. Effect of inlet velocity of the cathode gas on cell performance. Figure 8. Comparison of computed cell performance and experimental data (b: this model; 4: experimental data).

The results predicted with this model are compared with experimental data in public literature26 and displayed in Figure 8. The numerical model has been run under the conditions similar to those of the experimental case (i.e., temperature is 353.15 K, pressure is 5 and 3 atm at cathode and anode inlet, respectively). It should be pointed out that the reference transfer current density was adjusted24 to achieve the good agreement between predictions of the present model and the experimental data. However, it just indicates that our model can reflect fuel cell work relatively well, because the details on the membrane and electrode assembly (morphological and wetting characteristics) used in obtaining the experimental data are not known, and only average current densities and not current density distributions could be compared. The humidity of reactive gas has an important effect on the performance of the cell for the PEMFC operation. Computation was conducted for different cathode humidification temperatures of 40, 60, and 80 °C, keeping other parameters unchanged. The results are plotted in Figure 9. Under low current operation, performance will decrease for higher humidification temperature. This is because reactant concentrations are reduced, especially along the flow process for the reaction, rendering the reactants insufficient (i.e., dilution of reactants makes the performance decrease). However, under high current operation, performance will be enhanced at higher humidification temperatures. The reason is that water capacity at the cathode increases, which reduces water loss at the anode by increasing back diffusion. (26) Ticianelli, E. A.; Berry, J. G.; Srinivasan, S. J. Electroanal. Chem. 1988, 251, 275-295.

Figure 11. Cell performance varies with shoulder-to-channel ratio.

Thus, membrane resistance can be reduced, and conductivity of the membrane is improved effectively. The limiting current density is also higher for a higher humidification temperature. Figure 10 displays the performance of a fuel cell for a different cathode inlet velocity. It can be seen that the performance increases as inlet velocity increases, especially for high current operation. This can be explained by two side effects. On the one hand, high velocity can provide reactants faster by enhancing convection transport and then avoid starvation of reactants at the reactive area. On the other hand, the products can be removed in time for high gas velocity, and especially after condensation front, effective liquid water removal can improve fuel cell performance greatly. This model can predict the effect of geometry on the performance of a fuel cell, because it can describe spatial distributions of phase potential and activation overpotential in detail. Figure 11 shows the performance of a fuel cell for


different shoulder-to-channel ratios of 1:1, 1:2, and 1:3, keeping the total width of shoulder and channel constant. It can be seen that a small shoulder-to-channel ratio (i.e., fewer shoulders and more channels over a given electrode area) is favorable for the performance of a fuel cell, especially for high current operation. The reason for this phenomenon is that a wider fluid channel can effectively transport reactants and products to and from the reactive sites by decreasing the transport distance above the shoulder, and then the utilization of the reactive area is enhanced. That is, current density for a wider fluid channel case is higher than a narrower fluid channel in most areas, so as a whole the average current density is higher. Thus, cell performance is improved. 4. Conclusions A steady-state, three-dimensional, two-phase mathematical model of the PEMFC was developed in the framework of a computational fluid dynamics code. This model comprehensively considers almost all physical and chemical processes occurring in the fuel cell, such as flow, heat and mass transfer, charge transport, and electrochemical kinetics coupling with transport phenomena. Electron transports in the backing and catalyst layers, as well as protons in the membrane, are investigated by electrode electrochemical kinetics. The predicted distribution of current densities is radically different from those obtained with previous 2D and 3D models, where the variation in potential distribution was not fully accounted for. The multidimensional flow structure, distributions of reactants concentrations, activation overpotential, and current densities for the fuel cell have been discussed in detail. The polarization curve predicted by this model is close to experimental data in the published literature. The effects of inlet velocity and humidification temperatures of cathode stream and shoulderto-channel ratio on the cell performance are discussed in this work. The results indicate that an increase in velocity can enhance performance of the fuel cell. Inlet humidity conditions have a complex effect on the performance: under low current density operation, low humidity has a better performance; however, under high current density operation, high humidity is favorable for performance. Keeping the total width of shoulder and fluid channel constant, the large channel-to-shoulder ratio can increase the performance of the fuel cell effectively, especially for high current operation. This study can provide a theoretical reference for the design and operation of the fuel cell.

Energy & Fuels, Vol. 20, No. 2, 2006 747 Acknowledgment. We thank the Natural Science Foundation of Zhejiang Province (P. R. China, Grant No. 501140).

Nomenclature J‚kg-1‚k-1

c ) Specific heat, C ) Molar concentration, mol‚cm-3 D ) Diffusion coefficient, cm2‚s-1 F ) Faraday constant, 96 487 C/mol h ) Enthalpy H ) Total height of shoulder and channel, cm I ) Mean current density, A‚cm-2 j ) Transfer current density, A‚cm-3 J ) Mass flux K ) Hydraulic permeability, cm2 KrR ) Relative permeability of phase R, cm2 L ) Length of channel, cm n ) Electron number participating in the electrochemical reaction p ) Pressure, atm pc ) Capillary pressure r ) Stoichiometric coefficient or relative value R ) Universal gas constant, 8.314 J/(mol‚K) s ) Saturation T ) Cell temperature, K u ) Velocity vector, cm‚s-1 Greek letters β ) Electroosmotic drag coefficient φ ) Potential, V γ ) Concentration parameter µ ) Fluid viscosity, kg‚cm-1‚s-1 ) Porosity of the porous media material ϑ ) Interfacial tension, N/cm F ) Density of the fluid, mol‚cm-3 λ ) Relative mobility of phase and thermal conductivity σ ) Conductivity, Ω-1‚cm-1 ζ ) Transfer coefficient Superscripts eff ) Effective value ref ) Reference Subscripts a ) Anode c ) Cathode g ) Gas phase l ) Liquid phase k ) Chemical species R ) Phase (gas phase or liquid phase) EF050254B

A Three-Dimensional, Multicomponent, Two-Phase Model for a Proton

Recommend Documents