- Email: [email protected]
The effects of the cathode flooding on the transient responses of a PEM fuel cell
ARTICLE IN PRESS
Renewable Energy 33 (2008) 1824–1831 www.elsevier.com/locate/renene
The effects of the cathode flooding on the transient responses of a PEM fuel cell Mustapha Najjaria,b,, Faycel Khemilib, Sassi Ben Nasrallahb a
Institut Pre´paratoire aux e´tudes d’Inge´nieur, Av. Ibn Eljazzar, Monasir 5019, Tunisia Laboratoire d’Etudes des Syste`mes Thermiques et Energe´tiques, ENIM, Av. Ibn Eljazzar, Monasir 5019, Tunisia
b
Received 10 April 2007; accepted 12 October 2007 Available online 19 November 2007
Abstract This study investigates the effects of the flooding of the gas diffusion layer (GDL), as a result of liquid water accumulation, on the performance of a proton exchange membrane fuel cell (PEMFC). The transient profiles of the current generated by the cell are obtained using the numerical resolution of the transport equation for the oxygen molar concentration in the unsteady state. The dynamics of the system are captured through the reduction of the effective porosity of the GDL by the liquid water which accumulates in the void space of the GDL. The effects of the GDL porosity, GDL thickness and mass transfer at the GDL–gas channel interface on the evolution with time of the averaged current density are reported. The effects of the current collector rib on the evolution of the molar concentration of oxygen are also examined in detail. r 2007 Elsevier Ltd. All rights reserved. Keywords: Proton exchange membrane fuel cell; Gas diffusion layer; Flooding; Averaged current density
1. Introduction The proton exchange membrane fuel cell (PEMFC) is one of the most promising technologies (stationary power generation, small portable applications, powering cars, y) and will compete with many other types of energy devices. PEMFC is an electrochemical energy conversion device that converts hydrogen and oxygen into water and, in the process, it produces heat and electricity. The oxygen required for a fuel cell comes from the air pumped into the cathode. However, hydrogen has some limitations that make it impractical for use in most applications. One of the most critical aspects of a PEMFC is the water management. There is a delicate balance between the membrane hydration and avoiding cathode from flooding. Cathode flooding occurs when the water produced exceeds the water removal rate. Flooding in the cathode reduces the Corresponding author. Institut Pre´paratoire aux e´tudes d’Inge´nieur, Av. Ibn Eljazzar, Monasir 5019, Tunisia. Tel.: +216 97472458; fax: +216 73500512. E-mail address: [email protected] (M. Najjari).
0960-1481/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.renene.2007.10.003
oxygen transport to the reaction sites and decreases the effective catalyst area. In modern fuel cells, the cathode (or the anode) consists of: a catalyst layer, gas diffusion layer (GDL) and bipolar plates. The GDL is a porous medium electrically conductor. Liquid water that builds up at a fuel cell cathode decreases the performances and inhibits robust operation. A number of studies on this problem have been reported in the literature. Experimental work, for the case of a transparent PEMFC under a constant discharge mode, is reported by Tuber et al. [1]. Images of water formed inside the cathode gas channel are presented to explain the phenomena of water flooding. Results concerning the current density as a function of time are also presented. Theoretical analysis of the effects of the cathode flooding on the dynamics of the cell performance is reported by Wang et al. [2–4] and Chang et al. [5]. These models was based on the computational fluid dynamics approach, obtained by solving transport equations governing conservation of mass, momentum, species, energy, charge and taking account on a possible liquid–vapour phase-change. In this case, calculations become difficult to explore and
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
Nomenclature A C C ref Deff D0 Er F hl hm i i0 im JO lm M N el N in Nm N out Np NT p
electrode surface area ðm2 Þ oxygen concentration ðmol m3 Þ reference oxygen concentration ðmol m3 Þ effective diffusion coefficient ðm2 s1 Þ diffusion coefficient in pure gas phase ðm2 s1 Þ reversible cell potential (V) Faraday constant GDL thickness (m) interfacial mass transfer coefficient ðm2 s1 Þ local current density ðA m2 Þ reference current density ðA m2 Þ averaged current density ðA m2 Þ molar diffusion flux of oxygen membrane thickness (m) molar weight of water ðkg mol1 Þ mass of water produced by electrochemical reaction (kg) mass of water fed to the cell by humidification (kg) mass of water fed to the cathode through the membrane (kg) mass of water removed by gas convection (kg) net mass of water produced (kg) mass of liquid water corresponding to a complete flooding (kg) operating pressure (Pa)
numerical resolution is sensitive to different parameters (source terms, boundary conditions,y). However, simplified models are presented by Maggio et al. [6], Pasaogullari et al. [7], and Roshandel et al. [8]. In these models, an effective porosity of the GDL is defined which changes with the amount of liquid water generated by the PEMFC. This approach has been used to describe water transport in proton conductive membranes [6] and to study the effects of the compression on the porosity distribution [8]. Numerous theoretical models have used a diffusional approach for oxygen transport in the porous GDL [9,10]. Jeng [10] developed a pure diffusion steady-state equation that introduced equivalent oxygen diffusivity. The results show the influence of gas channel width and the GDL porosity. However, this study has considered a singlephase-flow and does not account for emergence of liquid water which accumulates in the GDL. The objective of the present paper is to provide a comprehensive model that describes the evolution of the current density under transient conditions in the presence of a partial flooding of the porous cathode. For this purpose, a two dimensional, transient, single-phase oxygen model is developed. Changes of the molar concentration of oxygen gas-phase with respect to the time and the liquid water rate production will be derived in terms of an effective GDL porosity.
psat R RH re Rcell S T t V cell wc wl x O2
1825
saturated pressure of the water (Pa) gas constant relative humidity electrical resistivity ðO m1 Þ ohmic resistance ðOÞ liquid water saturation temperature (K) time (s) cell voltage (V) gas flow channel width (m) GDL height (m) inlet molar fraction of oxygen
Greek symbols ac a dc 0 g Zc rl sm x
the cathode transfer coefficient net water transport coefficient through the membrane the thickness of the catalyst layer (m) initial porosity of the GDL effective gas porosity cathode activation overpotential (V) liquid water density ðkg m3 Þ protonic conductivity of the membrane ðO1 m1 Þ air stoichiometry
This article is organized as follows: first we formulate the problem, next the derived model is validated with experimental results. Results concerning the effects of various parameters such as the GDL thickness, the rib’s width and the mass transfer at the GDL–gas channel interface are presented and discussed.
2. Model development The present study focuses on oxygen mass transport, under transient conditions that takes place in the cathode of a PEMFC. The problem domain under consideration (Fig. 1) is confined to the porous GDL where the major mass transport limitations occur. The oxygen phase diffuses into the GDL from the humidified air supplied at the inlet of the gas flow channel. The dynamics of the system are captured through the dependence of the effective porosity on liquid water production from the electrochemical reaction at the catalyst layer. The model employs the following basic assumptions: (i) Homogenous porous media with uniform porosity. (ii) At the beginning, the vapour pressure in the GDL exceeds its saturated value. The product water condenses in the liquid phase.
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
1826
due to the electro-osmotic drag, diffusion of water from higher to lower concentration and convective flux due to the pressure gradient. The net water transport through the membrane can be written as a function of a net water transport coefficient a [13–15]. The amount of water crossing the membrane is then given by
y wl
Current collector rib
Membrane + catalyst layer
wc
Porous GDL
lm
0
iA Mt. (5) F A positive a means that the net water transport is from the anode side to the cathode side. The net mass of water produced at time t is
Gas flow channel
hl
Nm ¼ a
x
Fig. 1. Schematic of the problem showing coordinates and dimensions of the porous GDL.
N p ¼ N in þ N el þ N m N out .
(6)
2.2. Effective gas porosity (iii) The oxygen model flow rate is related to the gradient of its molar concentration. We assume equivalent oxygen diffusivity. (iv) Isothermal model for the PEMFC cathode.
g ¼ 0 ð1 SÞ,
2.1. Water transport in the cathode Water transport in the cathode of a PEMFC can originate from three sources: the reactant gases (air), electrochemical reaction at the catalyst layer and the net water transport through the membrane. The produced water by electrochemical reaction is related to the current density i and is equal to [11] iA Mt, (1) 2F where F is the Faraday constant. The amount of water that is fed to the cell by humidification is given by [11,12] N el ¼
N in ¼
iA x RHpsat Mt, 4F xO2 pRHpsat
The presence of liquid water in the porous GDL modifies the effective porosity for the oxygen gas transport. The effective gas porosity is determined from the following expression [6,8]:
(2)
where x is the air stoichiometry, xO2 is the molar fraction of oxygen in the inlet gas-phase, p is the operated pressure, RH is the relative humidity, A is the active area, M is the molar weight of water, t is the time and psat is the saturated pressure of the water given by 5096:23 psat ¼ exp 13:669 , (3) T where T is the temperature ð CÞ. However, amount of water may be removed by the gas convection in the gas flow channel. The amount of water carried out at the outlet is calculated according to [11,12] ! iA x psat N out ¼ 1 Mt. (4) 4F xO2 p psat The water transport from anode to cathode, through the membrane, is governed by three processes: flow of water
(7)
where 0 is the initial porosity of the GDL and S is the liquid water saturation given by S¼
Np . NT
(8)
N T is the mass of the liquid water corresponding to a complete flooding: N T ¼ rl 0 hl wl M.
(9)
2.3. Oxygen transport: diffusional approach In this model, we determine the molar diffusion flux of oxygen J O according to Fick’s law: ~ J~O ¼ Deff rC,
(10)
where C is the molar concentration of the oxygen and Deff its effective diffusion coefficient. The effective diffusion coefficient of the oxygen can be termed from the diffusion coefficient in the pure gas phase D0 and the void fraction of the GDL, g using Bruggeman’s correction [9] Deff ¼ ðg Þ3=2 D0 ¼ ð0 ð1 SÞÞ3=2 D0 .
(11)
The basic transport equation for oxygen molar concentration can be stated as qC ¼ Deff DC. (12) qt The dependence of Deff and g with respect to time is derived from the production of water in the GDL (Eqs. (1)–(9) and (11)). g
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
2.4. Boundaries conditions The GDL is connected to the membrane through catalyst layer CL corresponding to an infinitely thin layer and often modelled as an interface. The oxygen consumption rate by the electrochemical reaction is given by qC i Deff . (13) ¼ qx x¼0 4F The local current density in the catalyst layer is modelled by the Butler–Volmer equation and expressed by the Tafel kinetics [10,16]: Cð0; yÞ i ¼ ð1 SÞAi0 dc expðac F Zc =RTÞ, C ref
(14)
where i0 is the reference current density, dc is the thickness of the catalyst layer, C ref is the reference oxygen concentration and ac is the cathode transfer coefficient. In Eq. (14), Zc is the activation overpotential corresponding to losses that are associated with the kinetics of the electrochemical reaction. Neglecting the overpotential at the anode side, Zc is calculated as Zc ¼ E r V cell Zohm ,
(15)
where E r is the reversible cell potential (Nernst equation), V cell is the working cell potential (fixed in the present study) and Zohm is the ohmic overpotential due to losses from the protonic and the electronic resistances in the cell. Zohm is determined by the following equation: Zohm ¼ iRcell ,
lm þ r e hl , sm
(17)
where l m and sm are, respectively, the thickness and the protonic conductivity of the membrane and re is the electrical resistivity of the cathode paper. Since the protonic conductivity is dependant on water content in the membrane [17], it will be affected by the flooding effects of the cathode. The membrane conductivity will be determined using the following expression [6]: sm ¼ s0 ð1 SÞ,
(18)
where s0 is the initial conductivity of the membrane. It follows that the boundary condition at the GDLcatalyst layer interface is written as qC ¼ KC, (19) qx x¼0 where ð1 SÞAv i0 dc expðac F Zc =RTÞ K¼ . 4 FDeff C ref
Along the GDL–gas channel interface, the boundary condition is obtained from the convective mass transport analysis of the gas channel, using heat and mass transfer analogy [16,18]. The corresponding boundary condition is then specified as qC Deff ¼ hm ðCðx ¼ hl ; yÞ C 0 Þ, (23) qx x¼hl 0pypwc
where hm is the convective mass transfer coefficient and C 0 is the molar oxygen concentration in the channel. At all the remaining faces, impermeable boundary condition is assigned: qC ¼0 qx
at x ¼ 0 and 0pypwc ,
(24)
qC ¼0 qy
at y ¼ 0 and 0pxphl ,
(25)
qC ¼0 qy
at y ¼ wl and 0pxphl .
(26)
2.5. Numerical resolution The numerical resolution of Eq. (12) is based on the finite volume method [19]. For the present computations 61 51 uniform mesh has been used. Convergence is considered to be reached when the relative changes between successive iterations, at each time step Dt in C field, is less than 105 Dt ¼ 0:05 s is chosen for all numerical resolutions. 3. Results and discussions The developed two-dimensional transient model for the cathode GDL is used to examine the effects of many important factors on the transient phenomena in PEMFC system. Calculations were conducted for a PEMFC in a constant voltage discharge mode at ambient pressure. Inlet humidity and temperature are adjusted to 60% and 23 C, respectively. Geometrical and operating parameters are listed in Table 1. 3.1. Model validation
(20)
The local current density varies with y and is calculated by iðyÞ ¼ 4FKCð0; yÞ.
The average current density is then determined from Z 1 wl iðyÞ dy. (22) im ¼ wl 0
(16)
where Rcell is the ohmic resistance which includes the electronic and the membrane resistances: Rcell ¼ Rmembrane þ Relectrode ¼
1827
(21)
Before discussing the results, the model must be validated by comparing the numerical predictions with experimental data from the literature provided by Tuber et al. [1].
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
1828 Table 1 Physical properties and parameters
2800
1
1
Universal gas constant, R ðJ mol K Þ Water molar weight, M ðkg mol1 Þ Faraday constant, F ðC mol1 Þ Cell temperature, T (K) Operating pressure, p (Pa) Inlet molar fraction of oxygen, xO2 Oxygen diffusivity in pure gas phase, D0 ðm2 s1 Þ Air stoichiometry, x Relative humidity of inlet gas, RH Initial GDL porosity, 0 GDL thickness, hl (m) Gas flow channel width, wc (m) Convective mass transfer coefficient, hm ðm2 s1 Þ Cell voltage V cell (mV) Molar oxygen concentration in the gas channel, C 0 ðmol M3 Þ Reference current density times area, Ai0 ðA m3 Þ Thickness of the catalyst layer, dc (m) Cathode transfer coefficient, ac Initial conductivity of the membrane s0 ðO1 cm1 Þ Net water transport through membrane a Electrical resistivity of the cathode re ðmO cmÞ Membrane thickness l m (cm) Reversible cell potential E r (V) Referenceoxygen concentration, C ref ðmol m3 Þ
Reference
8.314 0.018 96 487 300 1 105 0.21 2 105
[1]
2 60% 0.5 0:23 103 103 5:395 102
[1] [1] [1] [1] [1] [1]
500 0.368
[1] Calculated
5 102
[10]
5 106 2 0.2
[10] [10] [6]
0.5 80
[13] Carbon paper TORAY Assumed Assumed Assumed
Averaged current density (A/m2)
Value
Experiments Tuber et al. present model
2400
2000
1600
1200 0
100
200
300
400
Time (s) Fig. 2. Transients profiles for current density from Tuber et al. [1] and the present model.
4000
0.015 1.17 0.368
Experiments were conducted by the visualization of the water build-up in the cathode of a transparent PEMFC to explain the influence of the water flooding under a constant discharge mode. Fig. 2 shows a comparison concerning the evolution of the average current density (Eq. (22)). It can be found that calculated values agree well with the experimental results in the first period of operating conditions. The current density decreases in the transient process. This can be explained by a partial flooding of the GDL which leads to a reduction of the effective porosity. However, experimental data are not well predicted by the present model for t4150 s. This discrepancy may be explained by: (i) the reduction of the liquid saturation due to the water that leaves the GDL to the gas flow channel and (ii) the effects of the gas flow channel flooding and the reduction of the corresponding cross-section [1] by produced water, which are not taken into account in this model. 3.2. Effects of the initial GDL porosity 0 The effects of the initial GDL porosity 0 on the evolution of the averaged current density are shown in Fig. 3 for three values of 0 . It is found that the current
e0 = 0.75 ; hl = 2.30 e-4 m e0 = 0.50 ; hl = 2.30 e-4 m e0 = 0.25 ; hl = 2.30 e-4 m
Averaged current density (A/m2)
Parameter
3000
2000
1000
0 0
50
100 150 Time (s)
200
250
Fig. 3. Effects of the GDL porosity on transients profiles of the averaged current density.
density is reduced if 0 decreases. In fact, a larger porosity 0 significantly reduces the resistance of oxygen transport through the GDL and greatly raises the oxygen concentration at the reaction sites due to the increase of the effective diffusion coefficient as indicated in Eq. (11). However, for large values of 0 , more oxygen is consumed and more water is produced. Since, the effects of flooding will be important in this case and the current density decreases with time more and more sharply as 0 increases.
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
3.3. Effects of the GDL thickness
Since the electrochemical reaction is sensitive to the oxygen concentration at the reaction sites, convective mass
4000
Averaged current density (A/m2)
hl = 1.15 e-4 m
5000
e0 = 0.85 ; hl = 0.90 e-4 m e0 = 0.85 ; hl = 2.30 e-4 m
4000
3000
2000 0
40
80
120 Time (s)
160
200
Fig. 5. Effects of the GDL thickness on the transient profiles of the averaged current density: 0 ¼ 0:85.
4000
Averaged current density (A/m2)
3.4. Effects of the convective mass transfer
Averaged current density (A/m2)
6000
Fig. 4 shows the influence of the GDL thickness hl on the corresponding dynamic results for the averaged current density. A thicker GDL leads to a higher resistance for the oxygen transport. As a result, the current density is reduced at high values of hl . However, the current density decreases rapidly with time for a thinner GDL. This is explained by the fact that a thicker GDL signifies a higher value of the void space in the porous GDL. Consequently, the rate of increase in the liquid saturation with generated liquid water is smaller for large hl . With further a decrease in hl value, the void space is reduced so that the liquid saturation increases sharply. Moreover, due to the filling of pores with the liquid water, the oxygen mass-transport to the reaction sites is reduced which leads to a weak current production by the electrochemical reaction. The transient behaviour of the averaged current density is then the result of the combined effects of the resistance to the oxygen mass transport and the flooding effects. As can be seen, in thinner GDL, the mass transport resistance decreases while the flooding effects are more pronounced. Fig. 5 shows the transient profiles of the averaged current density for different GDL thickness and at high values of 0 . In this case, it can be seen that the flooding rate is so important for hl ¼ 0:90 104 m that the current density values become lower than that corresponding to hl ¼ 2:30 104 m for t450 s.
3500
1829
hm = 5 x 5.395 e-2 mol / m s
3500
hm = 5.395 e-2 mol / m s hm = 0.5 x 5.395 e-2 mol / m s
3000
2500
2000
hl = 2.30 e-4 m hl = 3.45 e-4 m
1500 0
3000
50
100 150 Time (s)
200
250
Fig. 6. Transients profiles of the averaged current density for different values of the convective heat mass transfer coefficient.
2500
2000
1500 0
50
100 150 Time (s)
200
250
Fig. 4. Effects of the GDL thickness on transients profiles of the averaged current density: 0 ¼ 0:5.
transfer at the GDL–gas channel interface can have major impact on the flooding of the GDL. Fig. 6 shows the influence of the convective mass transfer coefficient hm on the transient profiles of the current density. These profiles show that under higher values of hm the generated current is increased, however, it drops slightly faster. The flooding degree in the GDL was increased while the produced water by the reaction is increased.
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
1830
In the case of low values of hm , the mass transfer of the reactants from the gas channel is too slow to produce a significant effect on blocking the GDL. For high values of hm , the consumption rate of oxygen increases. As a consequence, the porous GDL will be rapidly filled with the liquid water. 3.5. Effects of the current collector rib For all preceding results, we have ignored the portion covered by the rib ðwc ¼ wl Þ. In this section, we will consider the presence of a zone covered by the rib, which is represented by an impervious boundary at the GDL–gas channel interface. Presented results have been carried out using a constant gas channel width wc and changing the GDL width wl . As it can be seen in Fig. 7, for high values
of wl , the averaged current density decreases which is the result of the limitation of the oxygen mass transfer at the GDL–gas channel interface. Since the water generation is proportional to the current generation, a faster current drop is depicted for low values of wl . The current collector rib has a large impact on the local current density distribution and on the oxygen concentration field as shown in Figs. 8 and 9. The molar concentration of oxygen has lower values in the portion covered by the rib (Fig. 8). The local current is not evenly distributed on the GDL-catalyst layer interface and decreases from the portion facing the gas channel toward the portion covered by the rib (Fig. 9). As opposing to the averaged current density, the maximum value of the local current density increases with wl , reflecting the reduction of the flooding effects.
3000 3000
wl = 2 x wc wl = 5 x wc
Local current density (A/m2)
Averaged current density (A/m2)
wl = wc
2000
1000
wl = 5 x wc
2000
wl = 2 x wc ; t = 40 s wl = 5 x wc ; t = 400s wl = 2 x wc ; t = 400s
1000
0 0
100
200
300
400
500
Time (s)
0 0.00
0.20
0.40
0.60
0.80
1.00
(y/wl)
Fig. 7. Transient profiles of the averaged current density for different values of the gas channel width.
wl = 2 x wc
Fig. 9. Local current distribution on the interface GDL-catalyst layer.
wl = 5 x wc
Fig. 8. Oxygen molar concentration field for different values of wl t ¼ 400 s.
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
4. Conclusions Excess liquid water in the GDL may dramatically affect the cell performance. The purpose of this study is to predict the consequences of the flooding on the current generation as well as the reactant concentration. For this purpose, we developed a two dimensional model, based on the resolution of the diffusion equation. The dynamics of the system are captured through the dependence of the effective porosity on the liquid water production from the electrochemical reaction at the catalyst layer. Results from the simulations illustrate the evolution of the generated current under the influence of different parameters. We found that the cell performance may be affected by the flooding at; high porosity, thinner GDL and important convective mass transfer of the reactant at the GDL–gas channel interface. We also note that the presence of the current collector rib can reduce the flooding effects however it enhances the resistance to the oxygen diffusion. References [1] Tu¨ber K, Po´cza D, Hebling C. Visualization of water buildup in the cathode of a transparent PEM fuel cell. J Power Sources 2003;124(2): 403–14. [2] Um S, Wang CY. Three-dimensional analysis of transport and electrochemical reactions in polymer electrolyte fuel cells. J Power Sources 2004;125(1):40–51. [3] Wang Y, Wang CY. Transient analysis of polymer electrolyte fuel cells. Electrochim Acta 2005;50(6):1307–15. [4] Wang Y, Wang CY. Modeling polymer electrolyte fuel cells with large density and velocity changes. J Electrochem Soc 2005;152(2): A445–53.
1831
[5] Chang MH, Chen F, Teng HS. Effects of two-phase transport in the cathode gas diffusion layer on the performance of a PEMFC. J Power Sources 2006;160(1):268–76. [6] Maggio G, Recupero V, Pino L. Modeling polymer electrolyte fuel cells: an innovative approach. J Power Sources 2001;101(2):275–86. [7] Pasaogullari U, Wang CY. Two-phase transport and the role of micro-porous layer in polymer electrolyte fuel cells. Electrochim Acta 2004;49(25):4359–69. [8] Roshandel R, Farhanieh B, Saievar-Iranizad E. The effects of porosity distribution variation on PEM fuel cell performance. Renew Energy 2005;30(10):1557–72. [9] Baschuk JJ, Li X. Modelling of polymer electrolyte membrane fuel cells with variable degrees of water flooding. J Power Sources 2000; 86(1–2):181–96. [10] Jeng KT, Lee SF, Tsai GF, Wang CH. Oxygen mass transfer in PEM fuel cell gas diffusion layers. J Power Sources 2004;138(1–2):41–50. [11] Janssen GJM, Overvelde MLJ. Water transport in the protonexchange-membrane fuel cell: measurements of the effective drag coefficient. J Power Sources 2001;101(1):117–25. [12] Chang M-H, Chen F, Teng H. Effects of two-phase transport in the cathode gas diffusion layer on the performance of a PEMFC. J Power Sources 2006;160:268–76. [13] You L, Liu H. A two-phase flow and transport model for the cathode of PEM fuel cells. Int J Heat Mass Transfer 2002;45:2277–87. [14] You L, Liu H. A two-phase flow and transport model for PEM fuel cells. J Power Sources 2006;155:219–30. [15] Cai Y, Hu J, Ma H, Yi B, Zhang H. Effect of water transport properties on a PEM fuel cell operating with dry hydrogen. Electrochim Acta 2006;51(28):6361–6. [16] Pasaogullari U, Wang CY. Liquid water transport in gas diffusion layer of polymer electrolyte fuel cells. J Electrochem Soc 2004;151(3): A399–406. [17] Wang ZH, Wang CY, Chen KS. Two-phase flow and transport in the air cathode of proton exchange membrane fuel cells. J Power Sources 2001;94:40–50. [18] Springer TE, Zawodzinski TA, Gottesfeld S. Polymer electrolyte fuel cell model. J Electrochem Soc 1991;138(8):2334–42. [19] Patankar V. Numerical heat transfer fluid flow. New York: Hemisphere/MacGraw-Hill; 1980.
Renewable Energy 33 (2008) 1824–1831 www.elsevier.com/locate/renene
The effects of the cathode flooding on the transient responses of a PEM fuel cell Mustapha Najjaria,b,, Faycel Khemilib, Sassi Ben Nasrallahb a
Institut Pre´paratoire aux e´tudes d’Inge´nieur, Av. Ibn Eljazzar, Monasir 5019, Tunisia Laboratoire d’Etudes des Syste`mes Thermiques et Energe´tiques, ENIM, Av. Ibn Eljazzar, Monasir 5019, Tunisia
b
Received 10 April 2007; accepted 12 October 2007 Available online 19 November 2007
Abstract This study investigates the effects of the flooding of the gas diffusion layer (GDL), as a result of liquid water accumulation, on the performance of a proton exchange membrane fuel cell (PEMFC). The transient profiles of the current generated by the cell are obtained using the numerical resolution of the transport equation for the oxygen molar concentration in the unsteady state. The dynamics of the system are captured through the reduction of the effective porosity of the GDL by the liquid water which accumulates in the void space of the GDL. The effects of the GDL porosity, GDL thickness and mass transfer at the GDL–gas channel interface on the evolution with time of the averaged current density are reported. The effects of the current collector rib on the evolution of the molar concentration of oxygen are also examined in detail. r 2007 Elsevier Ltd. All rights reserved. Keywords: Proton exchange membrane fuel cell; Gas diffusion layer; Flooding; Averaged current density
1. Introduction The proton exchange membrane fuel cell (PEMFC) is one of the most promising technologies (stationary power generation, small portable applications, powering cars, y) and will compete with many other types of energy devices. PEMFC is an electrochemical energy conversion device that converts hydrogen and oxygen into water and, in the process, it produces heat and electricity. The oxygen required for a fuel cell comes from the air pumped into the cathode. However, hydrogen has some limitations that make it impractical for use in most applications. One of the most critical aspects of a PEMFC is the water management. There is a delicate balance between the membrane hydration and avoiding cathode from flooding. Cathode flooding occurs when the water produced exceeds the water removal rate. Flooding in the cathode reduces the Corresponding author. Institut Pre´paratoire aux e´tudes d’Inge´nieur, Av. Ibn Eljazzar, Monasir 5019, Tunisia. Tel.: +216 97472458; fax: +216 73500512. E-mail address: [email protected] (M. Najjari).
0960-1481/$ - see front matter r 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.renene.2007.10.003
oxygen transport to the reaction sites and decreases the effective catalyst area. In modern fuel cells, the cathode (or the anode) consists of: a catalyst layer, gas diffusion layer (GDL) and bipolar plates. The GDL is a porous medium electrically conductor. Liquid water that builds up at a fuel cell cathode decreases the performances and inhibits robust operation. A number of studies on this problem have been reported in the literature. Experimental work, for the case of a transparent PEMFC under a constant discharge mode, is reported by Tuber et al. [1]. Images of water formed inside the cathode gas channel are presented to explain the phenomena of water flooding. Results concerning the current density as a function of time are also presented. Theoretical analysis of the effects of the cathode flooding on the dynamics of the cell performance is reported by Wang et al. [2–4] and Chang et al. [5]. These models was based on the computational fluid dynamics approach, obtained by solving transport equations governing conservation of mass, momentum, species, energy, charge and taking account on a possible liquid–vapour phase-change. In this case, calculations become difficult to explore and
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
Nomenclature A C C ref Deff D0 Er F hl hm i i0 im JO lm M N el N in Nm N out Np NT p
electrode surface area ðm2 Þ oxygen concentration ðmol m3 Þ reference oxygen concentration ðmol m3 Þ effective diffusion coefficient ðm2 s1 Þ diffusion coefficient in pure gas phase ðm2 s1 Þ reversible cell potential (V) Faraday constant GDL thickness (m) interfacial mass transfer coefficient ðm2 s1 Þ local current density ðA m2 Þ reference current density ðA m2 Þ averaged current density ðA m2 Þ molar diffusion flux of oxygen membrane thickness (m) molar weight of water ðkg mol1 Þ mass of water produced by electrochemical reaction (kg) mass of water fed to the cell by humidification (kg) mass of water fed to the cathode through the membrane (kg) mass of water removed by gas convection (kg) net mass of water produced (kg) mass of liquid water corresponding to a complete flooding (kg) operating pressure (Pa)
numerical resolution is sensitive to different parameters (source terms, boundary conditions,y). However, simplified models are presented by Maggio et al. [6], Pasaogullari et al. [7], and Roshandel et al. [8]. In these models, an effective porosity of the GDL is defined which changes with the amount of liquid water generated by the PEMFC. This approach has been used to describe water transport in proton conductive membranes [6] and to study the effects of the compression on the porosity distribution [8]. Numerous theoretical models have used a diffusional approach for oxygen transport in the porous GDL [9,10]. Jeng [10] developed a pure diffusion steady-state equation that introduced equivalent oxygen diffusivity. The results show the influence of gas channel width and the GDL porosity. However, this study has considered a singlephase-flow and does not account for emergence of liquid water which accumulates in the GDL. The objective of the present paper is to provide a comprehensive model that describes the evolution of the current density under transient conditions in the presence of a partial flooding of the porous cathode. For this purpose, a two dimensional, transient, single-phase oxygen model is developed. Changes of the molar concentration of oxygen gas-phase with respect to the time and the liquid water rate production will be derived in terms of an effective GDL porosity.
psat R RH re Rcell S T t V cell wc wl x O2
1825
saturated pressure of the water (Pa) gas constant relative humidity electrical resistivity ðO m1 Þ ohmic resistance ðOÞ liquid water saturation temperature (K) time (s) cell voltage (V) gas flow channel width (m) GDL height (m) inlet molar fraction of oxygen
Greek symbols ac a dc 0 g Zc rl sm x
the cathode transfer coefficient net water transport coefficient through the membrane the thickness of the catalyst layer (m) initial porosity of the GDL effective gas porosity cathode activation overpotential (V) liquid water density ðkg m3 Þ protonic conductivity of the membrane ðO1 m1 Þ air stoichiometry
This article is organized as follows: first we formulate the problem, next the derived model is validated with experimental results. Results concerning the effects of various parameters such as the GDL thickness, the rib’s width and the mass transfer at the GDL–gas channel interface are presented and discussed.
2. Model development The present study focuses on oxygen mass transport, under transient conditions that takes place in the cathode of a PEMFC. The problem domain under consideration (Fig. 1) is confined to the porous GDL where the major mass transport limitations occur. The oxygen phase diffuses into the GDL from the humidified air supplied at the inlet of the gas flow channel. The dynamics of the system are captured through the dependence of the effective porosity on liquid water production from the electrochemical reaction at the catalyst layer. The model employs the following basic assumptions: (i) Homogenous porous media with uniform porosity. (ii) At the beginning, the vapour pressure in the GDL exceeds its saturated value. The product water condenses in the liquid phase.
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
1826
due to the electro-osmotic drag, diffusion of water from higher to lower concentration and convective flux due to the pressure gradient. The net water transport through the membrane can be written as a function of a net water transport coefficient a [13–15]. The amount of water crossing the membrane is then given by
y wl
Current collector rib
Membrane + catalyst layer
wc
Porous GDL
lm
0
iA Mt. (5) F A positive a means that the net water transport is from the anode side to the cathode side. The net mass of water produced at time t is
Gas flow channel
hl
Nm ¼ a
x
Fig. 1. Schematic of the problem showing coordinates and dimensions of the porous GDL.
N p ¼ N in þ N el þ N m N out .
(6)
2.2. Effective gas porosity (iii) The oxygen model flow rate is related to the gradient of its molar concentration. We assume equivalent oxygen diffusivity. (iv) Isothermal model for the PEMFC cathode.
g ¼ 0 ð1 SÞ,
2.1. Water transport in the cathode Water transport in the cathode of a PEMFC can originate from three sources: the reactant gases (air), electrochemical reaction at the catalyst layer and the net water transport through the membrane. The produced water by electrochemical reaction is related to the current density i and is equal to [11] iA Mt, (1) 2F where F is the Faraday constant. The amount of water that is fed to the cell by humidification is given by [11,12] N el ¼
N in ¼
iA x RHpsat Mt, 4F xO2 pRHpsat
The presence of liquid water in the porous GDL modifies the effective porosity for the oxygen gas transport. The effective gas porosity is determined from the following expression [6,8]:
(2)
where x is the air stoichiometry, xO2 is the molar fraction of oxygen in the inlet gas-phase, p is the operated pressure, RH is the relative humidity, A is the active area, M is the molar weight of water, t is the time and psat is the saturated pressure of the water given by 5096:23 psat ¼ exp 13:669 , (3) T where T is the temperature ð CÞ. However, amount of water may be removed by the gas convection in the gas flow channel. The amount of water carried out at the outlet is calculated according to [11,12] ! iA x psat N out ¼ 1 Mt. (4) 4F xO2 p psat The water transport from anode to cathode, through the membrane, is governed by three processes: flow of water
(7)
where 0 is the initial porosity of the GDL and S is the liquid water saturation given by S¼
Np . NT
(8)
N T is the mass of the liquid water corresponding to a complete flooding: N T ¼ rl 0 hl wl M.
(9)
2.3. Oxygen transport: diffusional approach In this model, we determine the molar diffusion flux of oxygen J O according to Fick’s law: ~ J~O ¼ Deff rC,
(10)
where C is the molar concentration of the oxygen and Deff its effective diffusion coefficient. The effective diffusion coefficient of the oxygen can be termed from the diffusion coefficient in the pure gas phase D0 and the void fraction of the GDL, g using Bruggeman’s correction [9] Deff ¼ ðg Þ3=2 D0 ¼ ð0 ð1 SÞÞ3=2 D0 .
(11)
The basic transport equation for oxygen molar concentration can be stated as qC ¼ Deff DC. (12) qt The dependence of Deff and g with respect to time is derived from the production of water in the GDL (Eqs. (1)–(9) and (11)). g
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
2.4. Boundaries conditions The GDL is connected to the membrane through catalyst layer CL corresponding to an infinitely thin layer and often modelled as an interface. The oxygen consumption rate by the electrochemical reaction is given by qC i Deff . (13) ¼ qx x¼0 4F The local current density in the catalyst layer is modelled by the Butler–Volmer equation and expressed by the Tafel kinetics [10,16]: Cð0; yÞ i ¼ ð1 SÞAi0 dc expðac F Zc =RTÞ, C ref
(14)
where i0 is the reference current density, dc is the thickness of the catalyst layer, C ref is the reference oxygen concentration and ac is the cathode transfer coefficient. In Eq. (14), Zc is the activation overpotential corresponding to losses that are associated with the kinetics of the electrochemical reaction. Neglecting the overpotential at the anode side, Zc is calculated as Zc ¼ E r V cell Zohm ,
(15)
where E r is the reversible cell potential (Nernst equation), V cell is the working cell potential (fixed in the present study) and Zohm is the ohmic overpotential due to losses from the protonic and the electronic resistances in the cell. Zohm is determined by the following equation: Zohm ¼ iRcell ,
lm þ r e hl , sm
(17)
where l m and sm are, respectively, the thickness and the protonic conductivity of the membrane and re is the electrical resistivity of the cathode paper. Since the protonic conductivity is dependant on water content in the membrane [17], it will be affected by the flooding effects of the cathode. The membrane conductivity will be determined using the following expression [6]: sm ¼ s0 ð1 SÞ,
(18)
where s0 is the initial conductivity of the membrane. It follows that the boundary condition at the GDLcatalyst layer interface is written as qC ¼ KC, (19) qx x¼0 where ð1 SÞAv i0 dc expðac F Zc =RTÞ K¼ . 4 FDeff C ref
Along the GDL–gas channel interface, the boundary condition is obtained from the convective mass transport analysis of the gas channel, using heat and mass transfer analogy [16,18]. The corresponding boundary condition is then specified as qC Deff ¼ hm ðCðx ¼ hl ; yÞ C 0 Þ, (23) qx x¼hl 0pypwc
where hm is the convective mass transfer coefficient and C 0 is the molar oxygen concentration in the channel. At all the remaining faces, impermeable boundary condition is assigned: qC ¼0 qx
at x ¼ 0 and 0pypwc ,
(24)
qC ¼0 qy
at y ¼ 0 and 0pxphl ,
(25)
qC ¼0 qy
at y ¼ wl and 0pxphl .
(26)
2.5. Numerical resolution The numerical resolution of Eq. (12) is based on the finite volume method [19]. For the present computations 61 51 uniform mesh has been used. Convergence is considered to be reached when the relative changes between successive iterations, at each time step Dt in C field, is less than 105 Dt ¼ 0:05 s is chosen for all numerical resolutions. 3. Results and discussions The developed two-dimensional transient model for the cathode GDL is used to examine the effects of many important factors on the transient phenomena in PEMFC system. Calculations were conducted for a PEMFC in a constant voltage discharge mode at ambient pressure. Inlet humidity and temperature are adjusted to 60% and 23 C, respectively. Geometrical and operating parameters are listed in Table 1. 3.1. Model validation
(20)
The local current density varies with y and is calculated by iðyÞ ¼ 4FKCð0; yÞ.
The average current density is then determined from Z 1 wl iðyÞ dy. (22) im ¼ wl 0
(16)
where Rcell is the ohmic resistance which includes the electronic and the membrane resistances: Rcell ¼ Rmembrane þ Relectrode ¼
1827
(21)
Before discussing the results, the model must be validated by comparing the numerical predictions with experimental data from the literature provided by Tuber et al. [1].
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
1828 Table 1 Physical properties and parameters
2800
1
1
Universal gas constant, R ðJ mol K Þ Water molar weight, M ðkg mol1 Þ Faraday constant, F ðC mol1 Þ Cell temperature, T (K) Operating pressure, p (Pa) Inlet molar fraction of oxygen, xO2 Oxygen diffusivity in pure gas phase, D0 ðm2 s1 Þ Air stoichiometry, x Relative humidity of inlet gas, RH Initial GDL porosity, 0 GDL thickness, hl (m) Gas flow channel width, wc (m) Convective mass transfer coefficient, hm ðm2 s1 Þ Cell voltage V cell (mV) Molar oxygen concentration in the gas channel, C 0 ðmol M3 Þ Reference current density times area, Ai0 ðA m3 Þ Thickness of the catalyst layer, dc (m) Cathode transfer coefficient, ac Initial conductivity of the membrane s0 ðO1 cm1 Þ Net water transport through membrane a Electrical resistivity of the cathode re ðmO cmÞ Membrane thickness l m (cm) Reversible cell potential E r (V) Referenceoxygen concentration, C ref ðmol m3 Þ
Reference
8.314 0.018 96 487 300 1 105 0.21 2 105
[1]
2 60% 0.5 0:23 103 103 5:395 102
[1] [1] [1] [1] [1] [1]
500 0.368
[1] Calculated
5 102
[10]
5 106 2 0.2
[10] [10] [6]
0.5 80
[13] Carbon paper TORAY Assumed Assumed Assumed
Averaged current density (A/m2)
Value
Experiments Tuber et al. present model
2400
2000
1600
1200 0
100
200
300
400
Time (s) Fig. 2. Transients profiles for current density from Tuber et al. [1] and the present model.
4000
0.015 1.17 0.368
Experiments were conducted by the visualization of the water build-up in the cathode of a transparent PEMFC to explain the influence of the water flooding under a constant discharge mode. Fig. 2 shows a comparison concerning the evolution of the average current density (Eq. (22)). It can be found that calculated values agree well with the experimental results in the first period of operating conditions. The current density decreases in the transient process. This can be explained by a partial flooding of the GDL which leads to a reduction of the effective porosity. However, experimental data are not well predicted by the present model for t4150 s. This discrepancy may be explained by: (i) the reduction of the liquid saturation due to the water that leaves the GDL to the gas flow channel and (ii) the effects of the gas flow channel flooding and the reduction of the corresponding cross-section [1] by produced water, which are not taken into account in this model. 3.2. Effects of the initial GDL porosity 0 The effects of the initial GDL porosity 0 on the evolution of the averaged current density are shown in Fig. 3 for three values of 0 . It is found that the current
e0 = 0.75 ; hl = 2.30 e-4 m e0 = 0.50 ; hl = 2.30 e-4 m e0 = 0.25 ; hl = 2.30 e-4 m
Averaged current density (A/m2)
Parameter
3000
2000
1000
0 0
50
100 150 Time (s)
200
250
Fig. 3. Effects of the GDL porosity on transients profiles of the averaged current density.
density is reduced if 0 decreases. In fact, a larger porosity 0 significantly reduces the resistance of oxygen transport through the GDL and greatly raises the oxygen concentration at the reaction sites due to the increase of the effective diffusion coefficient as indicated in Eq. (11). However, for large values of 0 , more oxygen is consumed and more water is produced. Since, the effects of flooding will be important in this case and the current density decreases with time more and more sharply as 0 increases.
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
3.3. Effects of the GDL thickness
Since the electrochemical reaction is sensitive to the oxygen concentration at the reaction sites, convective mass
4000
Averaged current density (A/m2)
hl = 1.15 e-4 m
5000
e0 = 0.85 ; hl = 0.90 e-4 m e0 = 0.85 ; hl = 2.30 e-4 m
4000
3000
2000 0
40
80
120 Time (s)
160
200
Fig. 5. Effects of the GDL thickness on the transient profiles of the averaged current density: 0 ¼ 0:85.
4000
Averaged current density (A/m2)
3.4. Effects of the convective mass transfer
Averaged current density (A/m2)
6000
Fig. 4 shows the influence of the GDL thickness hl on the corresponding dynamic results for the averaged current density. A thicker GDL leads to a higher resistance for the oxygen transport. As a result, the current density is reduced at high values of hl . However, the current density decreases rapidly with time for a thinner GDL. This is explained by the fact that a thicker GDL signifies a higher value of the void space in the porous GDL. Consequently, the rate of increase in the liquid saturation with generated liquid water is smaller for large hl . With further a decrease in hl value, the void space is reduced so that the liquid saturation increases sharply. Moreover, due to the filling of pores with the liquid water, the oxygen mass-transport to the reaction sites is reduced which leads to a weak current production by the electrochemical reaction. The transient behaviour of the averaged current density is then the result of the combined effects of the resistance to the oxygen mass transport and the flooding effects. As can be seen, in thinner GDL, the mass transport resistance decreases while the flooding effects are more pronounced. Fig. 5 shows the transient profiles of the averaged current density for different GDL thickness and at high values of 0 . In this case, it can be seen that the flooding rate is so important for hl ¼ 0:90 104 m that the current density values become lower than that corresponding to hl ¼ 2:30 104 m for t450 s.
3500
1829
hm = 5 x 5.395 e-2 mol / m s
3500
hm = 5.395 e-2 mol / m s hm = 0.5 x 5.395 e-2 mol / m s
3000
2500
2000
hl = 2.30 e-4 m hl = 3.45 e-4 m
1500 0
3000
50
100 150 Time (s)
200
250
Fig. 6. Transients profiles of the averaged current density for different values of the convective heat mass transfer coefficient.
2500
2000
1500 0
50
100 150 Time (s)
200
250
Fig. 4. Effects of the GDL thickness on transients profiles of the averaged current density: 0 ¼ 0:5.
transfer at the GDL–gas channel interface can have major impact on the flooding of the GDL. Fig. 6 shows the influence of the convective mass transfer coefficient hm on the transient profiles of the current density. These profiles show that under higher values of hm the generated current is increased, however, it drops slightly faster. The flooding degree in the GDL was increased while the produced water by the reaction is increased.
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
1830
In the case of low values of hm , the mass transfer of the reactants from the gas channel is too slow to produce a significant effect on blocking the GDL. For high values of hm , the consumption rate of oxygen increases. As a consequence, the porous GDL will be rapidly filled with the liquid water. 3.5. Effects of the current collector rib For all preceding results, we have ignored the portion covered by the rib ðwc ¼ wl Þ. In this section, we will consider the presence of a zone covered by the rib, which is represented by an impervious boundary at the GDL–gas channel interface. Presented results have been carried out using a constant gas channel width wc and changing the GDL width wl . As it can be seen in Fig. 7, for high values
of wl , the averaged current density decreases which is the result of the limitation of the oxygen mass transfer at the GDL–gas channel interface. Since the water generation is proportional to the current generation, a faster current drop is depicted for low values of wl . The current collector rib has a large impact on the local current density distribution and on the oxygen concentration field as shown in Figs. 8 and 9. The molar concentration of oxygen has lower values in the portion covered by the rib (Fig. 8). The local current is not evenly distributed on the GDL-catalyst layer interface and decreases from the portion facing the gas channel toward the portion covered by the rib (Fig. 9). As opposing to the averaged current density, the maximum value of the local current density increases with wl , reflecting the reduction of the flooding effects.
3000 3000
wl = 2 x wc wl = 5 x wc
Local current density (A/m2)
Averaged current density (A/m2)
wl = wc
2000
1000
wl = 5 x wc
2000
wl = 2 x wc ; t = 40 s wl = 5 x wc ; t = 400s wl = 2 x wc ; t = 400s
1000
0 0
100
200
300
400
500
Time (s)
0 0.00
0.20
0.40
0.60
0.80
1.00
(y/wl)
Fig. 7. Transient profiles of the averaged current density for different values of the gas channel width.
wl = 2 x wc
Fig. 9. Local current distribution on the interface GDL-catalyst layer.
wl = 5 x wc
Fig. 8. Oxygen molar concentration field for different values of wl t ¼ 400 s.
ARTICLE IN PRESS M. Najjari et al. / Renewable Energy 33 (2008) 1824–1831
4. Conclusions Excess liquid water in the GDL may dramatically affect the cell performance. The purpose of this study is to predict the consequences of the flooding on the current generation as well as the reactant concentration. For this purpose, we developed a two dimensional model, based on the resolution of the diffusion equation. The dynamics of the system are captured through the dependence of the effective porosity on the liquid water production from the electrochemical reaction at the catalyst layer. Results from the simulations illustrate the evolution of the generated current under the influence of different parameters. We found that the cell performance may be affected by the flooding at; high porosity, thinner GDL and important convective mass transfer of the reactant at the GDL–gas channel interface. We also note that the presence of the current collector rib can reduce the flooding effects however it enhances the resistance to the oxygen diffusion. References [1] Tu¨ber K, Po´cza D, Hebling C. Visualization of water buildup in the cathode of a transparent PEM fuel cell. J Power Sources 2003;124(2): 403–14. [2] Um S, Wang CY. Three-dimensional analysis of transport and electrochemical reactions in polymer electrolyte fuel cells. J Power Sources 2004;125(1):40–51. [3] Wang Y, Wang CY. Transient analysis of polymer electrolyte fuel cells. Electrochim Acta 2005;50(6):1307–15. [4] Wang Y, Wang CY. Modeling polymer electrolyte fuel cells with large density and velocity changes. J Electrochem Soc 2005;152(2): A445–53.
1831
[5] Chang MH, Chen F, Teng HS. Effects of two-phase transport in the cathode gas diffusion layer on the performance of a PEMFC. J Power Sources 2006;160(1):268–76. [6] Maggio G, Recupero V, Pino L. Modeling polymer electrolyte fuel cells: an innovative approach. J Power Sources 2001;101(2):275–86. [7] Pasaogullari U, Wang CY. Two-phase transport and the role of micro-porous layer in polymer electrolyte fuel cells. Electrochim Acta 2004;49(25):4359–69. [8] Roshandel R, Farhanieh B, Saievar-Iranizad E. The effects of porosity distribution variation on PEM fuel cell performance. Renew Energy 2005;30(10):1557–72. [9] Baschuk JJ, Li X. Modelling of polymer electrolyte membrane fuel cells with variable degrees of water flooding. J Power Sources 2000; 86(1–2):181–96. [10] Jeng KT, Lee SF, Tsai GF, Wang CH. Oxygen mass transfer in PEM fuel cell gas diffusion layers. J Power Sources 2004;138(1–2):41–50. [11] Janssen GJM, Overvelde MLJ. Water transport in the protonexchange-membrane fuel cell: measurements of the effective drag coefficient. J Power Sources 2001;101(1):117–25. [12] Chang M-H, Chen F, Teng H. Effects of two-phase transport in the cathode gas diffusion layer on the performance of a PEMFC. J Power Sources 2006;160:268–76. [13] You L, Liu H. A two-phase flow and transport model for the cathode of PEM fuel cells. Int J Heat Mass Transfer 2002;45:2277–87. [14] You L, Liu H. A two-phase flow and transport model for PEM fuel cells. J Power Sources 2006;155:219–30. [15] Cai Y, Hu J, Ma H, Yi B, Zhang H. Effect of water transport properties on a PEM fuel cell operating with dry hydrogen. Electrochim Acta 2006;51(28):6361–6. [16] Pasaogullari U, Wang CY. Liquid water transport in gas diffusion layer of polymer electrolyte fuel cells. J Electrochem Soc 2004;151(3): A399–406. [17] Wang ZH, Wang CY, Chen KS. Two-phase flow and transport in the air cathode of proton exchange membrane fuel cells. J Power Sources 2001;94:40–50. [18] Springer TE, Zawodzinski TA, Gottesfeld S. Polymer electrolyte fuel cell model. J Electrochem Soc 1991;138(8):2334–42. [19] Patankar V. Numerical heat transfer fluid flow. New York: Hemisphere/MacGraw-Hill; 1980.