Effects of anisotropic permeability and electrical conductivity of gas diffusion layers on the performance of proton exchange membrane fuel cells

Applied Energy 95 (2012) 50–63

Contents lists available at SciVerse ScienceDirect

Applied Energy journal homepage: www.elsevier.com/locate/apenergy

Effects of anisotropic permeability and electrical conductivity of gas diffusion layers on the performance of proton exchange membrane fuel cells M.S. Ismail ⇑, K.J. Hughes, D.B. Ingham, L. Ma, M. Pourkashanian Centre for Computational Fluid Dynamics, Energy Technology and Innovation Initiative (ETII), University of Leeds, Leeds LS2 9JT, UK

a r t i c l e

i n f o

Article history: Received 25 July 2011 Received in revised form 17 January 2012 Accepted 3 February 2012 Available online 14 March 2012 Keywords: PEM fuel cells Gas diffusion layers Gas permeability Electrical conductivity CFD model Anisotropic transport properties

a b s t r a c t A 3-dimensional model for an in-house proton exchange membrane (PEM) fuel cell with serpentine channels has been developed in order to investigate the sensitivity of the fuel cell performance to the anisotropic gas permeability and electrical conductivity of gas diffusion layers (GDLs). For a realistic range of transport properties being investigated, the fuel cell performance was found to be very sensitive to the electrical conductivity but almost insensitive to the gas permeability of the GDL. For the given operating conditions, the current density was found to be a maximum in the vicinity of the edge between the flow channel and the rib of the current collector. Since the most common GDL materials present a rather significant anisotropy in the in-plane directions, the effects of such anisotropy has been evaluated. Given that the through-plane conductivity is maintained constant for all the cases investigated, for a realistic range of the in-plane electrical conductivity, the fuel cell performance was found to be almost insensitive to this parameter. Therefore such anisotropy can be practically ignored. Finally, for single phase operating conditions, the U-bend in the serpentine channel has no effect on the overall performance of the fuel cell. Hence, only a straight channel of the fuel cell may be modelled and used as a quick performance indicator. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Over the past few years, there have been numerous computational fluid dynamics (CFD) models for proton exchange membrane (PEM) fuel cells. The benefits of these models, among others, are that they act as a fast and cost-effective design and optimisation tool, and assist in a better understanding of the physics associated with the transport of reactant gases and liquid water within and at the interfaces between the several components of the fuel cell. Notably, in most CFD models, a considerable number of the parameters are frequently assumed or theoretically estimated. For example, the transport properties of the gas diffusion layers (GDLs) such as the thermal conductivity and the gas permeability are assumed to be isotropic [1–10]. This assumption is questionable as the structure of the non-woven GDLs (often known as carbon papers), the most commonly used GDL [11], dictates a significant anisotropy in its transport properties. They are mainly made from carbon fibres that are normally oriented in the in-plane directions, see Fig. 1. Therefore, it is expected that the transport properties in the in-plane directions are different from those in the through-plane direction. ⇑ Corresponding author. Tel.: +44 113 343 2512; fax: +44 113 246 7310. E-mail address: [email protected] (M.S. Ismail). 0306-2619/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.apenergy.2012.02.003

It must be stressed that most of the developers of the CFD models have adopted such a simplification because of the lack of experimental data concerning the in-plane transport properties of the GDLs. As reported in our previous papers [12–15], the gas permeability and electrical conductivity in the in-plane directions are at least one order of magnitude higher than those in the through-plane direction. This significant anisotropy in the transport properties of GDLs has also been confirmed by other investigators [16–18]. There have been a number of numerical investigations in the literature that have studied the effect of the GDL anisotropy on the global performance of PEM fuel cells and on the distribution of the key variables within and at the interfaces between the components of the fuel cell. Bapat and Thynell [19,20] developed a 2-dimensional model to examine the effect of the anisotropic electrical conductivity and thermal conductivity on the distribution of the current density and the temperature within the cathode electrodes. They [19] showed that, for a given through-plane electrical conductivity, as the in-plane electrical conductivity increases, the local current density increases over the channel and decreases over the rib. On the other hand, for a given in-plane electrical conductivity, the current density over both the rib and the channel increase with increasing through-plane conductivity. They [20] also found that, for a given moderately high temperature, a high current density is obtained for low through-plane and high in-plane thermal

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63

51

Nomenclature a ac b Cp D E Er F h i io ~ J K M n p pv psat ~ q R

specific area, m2 water activity Tafel slope, V specific heat capacity, J kg1 K1 diffusion coefficient, m2 s1 cell potential, V equilibrium cell potential, V Faraday constant, C mol1 specific enthalpy, J kg1 current density, A m2 exchange current density, A m2 diffusion flux, kg m2 s1 permeability, m2 molecular weight, kg mol1 number of transferred electrons per mole pressure, Pa partial pressure of water vapour, Pa saturation pressure of water, Pa heat flux, W m2 universal gas constant, J mol1 K1

conductivity. Meng and Wang [21] developed a 3-dimensional model for a typical PEM fuel cell. They showed that, at high cell voltage, the minimum current density always occurs under the middle of the gas channel and increases towards the rib of the collector. At low cell voltage, except for the inlet section of the cell, the maximum current density occurs under the gas channel and decreases towards the rib [21]. Pharoah et al. [22] developed a 2-dimensional model for a PEM fuel cell cathode and showed that the local current density and the overall fuel cell performance are underestimated if the anisotropy in the oxygen diffusivity is ignored. Also, they showed that ignoring the anisotropy of the GDL electrical conductivity results in shifting the maximum local current density from the region over the channel to the regions over the rib of the collector. However, they [22] noticed that the polarisation curve of the ‘anisotropic electrical conductivity’ case is very similar to that of the ‘isotropic electrical conductivity’ case; this signifies the importance of equally investigating the global and local performance of the fuel cell. Zhou and Liu [23] developed a 3-dimensional PEM fuel cell model to examine the effect of the anisotropic electrical conductivity on the distribution of the local current and the overall performance of the fuel cell. Given that realistic values for the through- and

Fig. 1. A scanning image of a typical GDL.

S T u Y

source term in the conservation equations temperature, K velocity, m s1 mass fraction, –

Greek symbols e porosity of the porous medium b inertial coefficient, m1 q density, kg m3 r electrical conductivity, S m1 rm ionic conductivity, S m1 /s solid phase potential, V /m membrane phase potential, V j thermal conductivity, W m1 K1 l dynamic viscosity, Pa s a transfer coefficient gact activation overpotential, V s shear stress, Pa n tortuosity

in-plane electrical conductivity are used, they showed that the highest local current density is at the boundary between the channel and the rib. Further, Zhou and Liu [23] argued that the average current density is larger over the channel where there is relatively higher amount of oxygen. Zhang et al. [24] developed a 2-dimensional PEM fuel cell model and showed that if the in-plane transport properties are lower than the through-plane ones, then the output power is lower than that of the base case where the transport properties are assumed to be isotropic. Also, they found that the performance of the fuel cell slightly improves with increasing in-plane gas permeability and in-plane thermal conductivity, and significantly improves with increasing in-plane mass diffusivity. However, they stated that the increase in the in-plane electrical conductivity has almost no effect on the performance of the fuel cell but it has an effect on the distribution of the local current density. Apart from the work of Zhou and Liu, there have been no investigations that have conducted a parametric study to investigate the effect of the anisotropy in the transport properties of GDLs on the fuel cell performance using experimentally characterised and realistic anisotropic transport properties. Also, it has been observed that these investigations were conducted for only a single, straight short channel. For serpentine flow configuration, such a geometry may not be appropriate since the convection flow from one turn to another is not captured [25]. Furthermore, there have been no investigations that have thoroughly studied the effect of the gas permeability of GDLs on the fuel cell performance and on the local distribution of the key variables. In this investigation, the effects of the anisotropy in the permeability and the electrical conductivity of the GDL on the global and local performance of the fuel cell has been investigated. In order to achieve this, the previously characterised through- and in-plane permeability and electrical conductivity, as reported in Refs. [12–15], have been employed in a CFD model for an in-house PEM fuel cell. Furthermore, the effect of the U-bend of the serpentine channel on the performance of the fuel cell has been probed. The layout of this paper is as follows: the in-house PEM fuel cell and the corresponding CFD model are described in Sections 2 and 3, respectively. Section 4 presents and discusses the results of the run simulations. The main conclusions of this study are presented in Section 5.

52

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63

2. Experimental procedure As described in Section 2.1, a commonly-used GDL, namely SGL 10BC, has been used in fabricating a membrane electrode assembly (MEA). This MEA has been used in an assembled single PEM fuel cell as described in Section 2.2. 2.1. MEA fabrication For the cathode electrode, the amount of the catalyst, 60% Pt/C (Alfa Aesar, UK), was calculated based on 0.4 mg cm2 platinum loading and 11.56 cm2 active area. The required amount of the catalyst was first ultrasonically mixed with a few drops of deionised water for a few seconds in order to break the catalyst powder into small pieces and prevent it from being deactivated when blending it with the dispersion agent, ethanol. A calculated amount of 50 wt.% ethanol (calculated based on 15 ml ethanol per 50 mg Pt/C) was added to the wet catalyst and the resulting mixture was ultrasonically blended for about 3 min. The resulting blend was then mixed with an amount of 5 wt.% Nafion solution (Sigma Aldrich, UK) that was calculated based on 0.1 mg Nafion cm2. This final mixture was ultrasonically blended for 15 min to form what is normally known as an ‘electrode ink’. The electrode ink was then manually sprayed onto the surface of the GDL, mounted on a PTFE plate heated-to-75 °C until the desired platinum loading, 0.4 mg cm2, was realised. The same procedures were repeated for the anode side which has slightly different requirements for the catalyst layer, namely 0.20 and 0.25 mg cm2 for the platinum and Nafion loadings, respectively. It should be noted that the catalyst used for the anode side was different to that used for cathode side, namely 20% Pt/C (Alfa Aesar, UK). Both the cathode and anode electrodes, i.e. the catalysed GDLs, were positioned on the faces of the 115 Nafion membrane and hot-pressed (using a hydraulic press) at 130 °C and 50 kg cm2 for 3 min to form a unit of MEA. 2.2. Electrochemical tests A single PEM fuel cell was assembled by sandwiching the MEA between two composite graphite bipolar plates (Bac2, UK) with an 11-turn and single-pass serpentine flow channel of depth 2 mm and width 2 mm, and a rib of width 0.8 mm. The fuel cell, after being placed between two copper end plates, was connected to an in-house fuel cell test station and fed with humidified hydrogen and air at the anode and cathode sides respectively. The flow rates of both the hydrogen and air were kept constant equivalent to a stoichiometry of 2 which was calculated based on 500 mA cm2 current density. Hydrogen and air gases were humidified by flowing them through bubbler-type humidifiers. The temperatures of the inlet gases and the fuel cell were maintained at 30 °C. The pressures on the cathode and anode sides of the fuel cell were both 2.5 bar gauge. The polarisation curve of the fuel cell was obtained using an electronic load, LD300 DC (TTI, UK). The fuel cell was run twice, with a separation period of a day, and the respective polarisation curves were found to be repeatable.

3. Model formulation The computational domain was limited to one turn of a serpentine channel as shown in Fig. 2. This has been done since the domain of the entire fuel cell is computationally very expensive and, as a consequence, takes a significantly longer time to achieve convergence. For the sake of simplicity, the model assumes the following:

 The fuel cell operates under steady-state conditions.  The flow in the channels is laminar since the Reynolds number is small.  Water exists in vapour phase only, i.e. a single phase water transport is assumed.  The compression applied on the fuel cell components is uniform. 3.1. Governing equations The 3-dimensional model being developed in this investigation accounts for (i) heat transport within the fuel cell, (ii) mass, momentum, and species in the flow channels, gas diffusion layers, and the catalyst layers, (iii) electrochemical reactions in the catalyst layer, (iv) electronic conduction in the electrically conducting solid parts of the fuel cell, and (v) ionic conduction in the membrane. The transport phenomena in a fuel cell are governed by a number of conservation equations. After consideration of the assumptions mentioned above and the substitution of the source terms, these equations are as follows [26–30]: 3.1.1. Conservation of mass, momentum, and species equations

Mass equation : r  ðeq~ uÞ ¼ 0 Momentum equation : r  ðeq~ u~ uÞ 2 ~ e lu 3 ~~ ¼ erp þ r  ðlre~ þ e qbuu uÞ þ K uY k Þ ¼ r  ðqDeff Species equation : r  ðeq~ k rY k Þ þ Sk

ð1Þ

ð2Þ ð3Þ

where q is the fluid density, p is the pressure, l is the dynamic viscosity of the fluid, ~ u is the volume-averaged velocity vector, e, K, b are the porosity, the permeability, and the inertial coefficient of the porous medium, Yk and Deff k are the mass fraction and the effective mass diffusivity of species k, and Sk is the source term of species k resulting from the electrochemical reactions. It should be noted that the porosity of the flow channels is set to unity. Deff k can be calculated using a number of correlations, such as those described by Tomadakis and Sotirchos, Nam and Kaviany, and Bruggemann [31,32]. The latter is typically used by the CFD-ACE+ software to calculate the effective mass diffusivity: n Deff k ¼ e Dk

ð4Þ

where n is the tortuosity of the porous medium and Dk is the ordinary diffusion coefficient of species k. 3.1.2. Electrochemical reactions Hydrogen and oxygen are consumed in the anode and the cathode catalyst layers, respectively. The corresponding sink terms in Eq. (3) are given by:

SH2 ¼ 

ia aa M H2 2F

ð5Þ

SO2 ¼ 

ic ac M O2 4F

ð6Þ

On the other hand, water is produced at the cathode and therefore the corresponding source term in Eq. (3) is given by:

SH2 O ¼

ic ac M H2 O 2F

ð7Þ

where ia and ic are the current densities for the reactions at the anode and cathode electrodes, respectively, aa and ac are the specific areas of the anode and cathode electrodes, respectively, M is the

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63

53

Fig. 2. (a) A schematic for the operating PEM fuel cell, and (b) the computational domain.

molecular weight, and F is the Faraday’s constant (96,485 C/mol). The specific area is calculated as follows [33]:

a¼

Catalyst loading  Surface area Thickness of catalyst layer

ð8Þ

The surface areas for the cathode and anode catalysts were provided by the supplier (Alfa Aesar, UK) as 440 and 1185 cm2 mg1 respectively. As mentioned in Section 2, the catalyst loadings for the anode and cathode electrodes were 0.2 and 0.4 mg cm2 respectively. The current densities are obtained from the Butler– Volmer equation:

ia aa ¼

ref ia aa

ref

ic ac ¼ ic ac

Y H2

!0:5 

Y ref H2 Y O2 Y ref O2



exp ! exp

aa;a F RT



ac;a F RT





gact;a  exp  



gact;c  exp 

aa;c F RT

ac;c F RT

gact;a

gact;c

 ð9Þ

 ð10Þ

where iref is the exchange current density at a reference temperature and pressure of 25 °C and 1 atm, respectively, and aa,a and aa,c are the anodic and the cathodic transfer coefficients for the reaction at the anode, and ac,a and ac,c are the anodic and the cathodic transfer coefficients for the reaction at the cathode. The values of the exchange current density and the transfer coefficients are

ref listed in Table 1. Y ref H2 and Y O2 are the mass fractions for hydrogen and oxygen at a reference temperature and pressure of 25 °C and 1 atm, respectively. R is the universal gas constant (8.314 J mol1 K1), and T is the temperature. The activation overpotential, gact, is given by [23,26]:

gact ¼ /s  /m  Eo

ð11Þ

where /s is the solid phase potential and it is given by Eq. (13), /m is the membrane phase potential and it is given by Eq. (14), and Eo is the reference potential of the electrode. Eo is zero at the anode side and equal to the equilibrium cell potential, Er, at the cathode side [34]:

  Er ¼ 1:482  0:000845T þ 0:0000431T ln pH2 p0:5 O2

ð12Þ

3.1.3. Conservation of charge equations The electronic and ionic conductions are given by:

r  ðrs r/s Þ ¼ S/;s

ð13Þ

r  ðrm r/m Þ ¼ S/;m

ð14Þ

where rs, /s, and S/,s are the electrical conductivity, electrical potential, and source term of the solid phase, respectively, and rm, /m, and S/,m are the electrical conductivity, electrical potential,

54

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63 Table 1 Geometrical, operational, and physical properties for the base case. Property

Value

Channel length Channel height Channel width Land area width GDL thickness Catalyst layer thickness Membrane thickness Operating temperature Gauge pressure at anode Gauge pressure at cathode Relative humidity of inlet gases Oxygen/nitrogen molar ratio in air Catalyst layer porosity GDL porosity GDL through-plane permeability GDL in-plane permeability GDL through-plane inertial coefficient GDL in-plane inertial coefficient GDL through-plane electrical conductivity GDL in-plane electrical conductivity Catalyst layer permeability Electrical conductivity of catalyst layers Membrane permeability Thermal conductivity of GDLs Thermal conductivity of catalyst layers Thermal conductivity of membrane Faraday constant Universal gas constant Active area Density of current collectors Thermal conductivity of current collector Specific heat capacity of current collectors Electrical conductivity of current collectors Tortuosity for GDLs Tortuosity for catalyst layers Tortuosity for membrane Anode inlet gas velocity Cathode inlet gas velocity Anode inlet mass fraction of hydrogen Anode inlet mass fraction of water Cathode inlet mass fraction of oxygen Cathode inlet mass fraction of nitrogen Cathode inlet mass fraction of water Hydraulic diameter Density of gases, kg m3 Viscosity of gases, Pa s Diffusion coefficient, m2 s1 Specific heat capacity of gases Thermal conductivity of gases

2.8  102 m 2.0  103 m 2.0  103 m 8.0  104 m 3.0  104 m 3.0  105 m 1.5  104 m 303 K 2.5 bar 2.5 bar 100% 0.21/0.79 0.4 [34] 0.7a 4.97  1013 m2 [14] 1.87  1012 m2 [12] 4.22  107 m1 [14] b 4.05  106 m1 [12]b 48 S/mc 4000 S m1 [13] 1  1013 m2 [38] 300 S m1 [39] 1.8  1018 m2 [37] 200 W m1 K1 [40] 200 W m1 K1 [40] 200 W m1 K1 [40] 96,485 C mol1 8.314 J mol1 K1 11.56  104 m2 1860 kg m3 [41] 200 W m1 K1 [40] 865 J kg1 K1 [41] 3200 S m1 [41] 1.5 [40] 1.5 [40] 5 [40] 0.42 m s1 1.06 m s1 0.37 0.63 0.22 0.72 0.06 2  103 m Computed using Ideal Gas Law [40]d Computed using Mix Kinetic Theory [40]d Computed using Schmidt number = 0.7 [40]d Computed using Mix JANNAF Method [40]d Computed using Prandtl number = 0.707 [40]d 1.93  105 A m3e

Reference exchange current density  specific area at cathode, ac iref o;c Reference exchange current density  specific area at anode, aa io;a

0.512e 9.3  108 A m3 [40]

Transfer coefficients for anode reaction

0.5 [40]

Transfer coefficients for cathode reaction ref

a

Calculated using Eq. (18) in [15]. The areal weight of GDL, which is required to be substituted in the above equation, was measured as 160 g m2. b The effects of the inertial coefficient have been investigated and found to have negligible effects on the fuel cell performance as well as the distribution of the variables investigated in this study. c Calculated following the experimental procedures described in Section 2.3.3 in [13]. Since there are no inputs for contact resistances in the CFD-ACE+, the through-plane conductance was assumed to be only due to the ‘bulk’ conductances of the current collector and the GDL. Based on this, the through-plane conductivity of the GDL material was calculated as 48 S m1. d Details of the equations are available in Appendix A. e Obtained by nonlinearly regressing the experimental data shown in Fig. 4 to Eq. (1) in [42]. This regression was performed under an assumption that the activation losses associated with the anode reaction, hydrogen oxidation reaction (HOR), are negligible and this is due the considerably fast kinetics of the HOR compared to that of the cathode reaction, oxygen reduction reaction (ORR). For a realistic range of kinetic parameters of the HOR, the fuel cell performance was found to be insensitive to the exchange current density and the transfer coefficients of the HOR [43].

and source term of the ionic (membrane) phase, respectively. The source terms are non-zero in the catalyst layers and equal to either iaaa (at the anode side) or icac (at the cathode side). They are zero elsewhere, i.e. GDLs, current collectors, and the membrane. The

ionic conductivity for Nafion membranes rm was measured and fitted by Springer et al. [35] as follows:





rm ¼ ð0:005139k  0:00326Þ exp 1268

1 1  303 T



ð15Þ

55

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63

where k is the membrane water conductivity and is given by:

8 2 3 > < 0:043 þ 17:18ac  39:85ac þ 36:09ac k ¼ 14 þ 1:4ðac  1Þ for 1 < ac 6 3 > : 16:8 for ac > 3

for 0 < ac 6 1

ð16Þ The water activity, ac, is given by:

ac ¼

pv psat

ð17Þ

where pv is the partial pressure of water vapour, and psat is the saturation pressure of water which can be obtained from the following formula [35]:

Y

log psat ¼ 2:1794 þ 0:02953ðT  273:15Þ  9:1837  105 ðT  273:15Þ2 þ 1:4454  107 ðT  273:15Þ3

ð18Þ

Z

X

3.1.4. Conservation of energy The conservation equation for energy may be written as follows:

r  ðeq~ q þ e~ sr  ~ u  ðiaÞgact þ uhÞ ¼ r  ~

j~i ~ij

r

ð19Þ

where h is the specific enthalpy of the gas mixture and ~ s is the shear stress tensor. The heat flux, ~ q, is given by:

~ q ¼ jeff rT þ

n X ~ J k hk

ð20Þ

k¼1

where n is the total number of gas phase species in the system, and ~ Jk and hk are the diffusion flux and enthalpy of species k, respectively. The effective thermal conductivity, jeff, can be theoretically calculated using a number of correlations such as the ones developed by Sadeghi et al. [36] and Dagan [27]. The latter was used in this study:

jeff ¼ 2jS þ



e 2jS þ jF

þ

1 1e 3jS

ð21Þ

where jS and jF are the thermal conductivities of the solid matrix and fluid regions (pores) respectively. 3.2. Boundary conditions As mentioned in Section 2, the stoichiometric flow ratio for hydrogen and air was fixed at 2 based on a current density of 500 mA/cm2. Therefore, the inlet velocities for the humidified gases were calculated to be 0.42 and 1.06 m/s, respectively. It should be noted that both the fuel and air entered the fuel cell at 100% relative humidity and the inlet mass fractions for all the gaseous species are listed in Table 1. Outlet pressure boundary conditions were used at the outlets and symmetry boundary conditions were prescribed at the left and right hand sides of the domain shown in Fig. 2b, namely at x = 0 mm and x = 5.60 mm. On all walls, the no-slip boundary condition is applied for the momentum equation. The above set of conservation equations were solved by adopting a finite-volume scheme on a structural grid within the framework of a commercial CFD software (CFD-ACE+). The number of control volumes for the base case was about 100,000, see Fig. 3, and this was changed to about 200,000 in order to check for mesh independency. The average current density at 0.55 V was calculated for both grids and a negligible relative error was obtained, less than 0.1%. Therefore, the computational mesh with only

Fig. 3. The discretised computational domain for the base case (about 100,000 control volumes).

100,000 control volumes was used in order to save computation time without a significant loss in accuracy. In the computations, the iteration for the solutions progresses until the following convergence criteria is satisfied: 6 l jðwli;j;k  wl1 i;j;k Þ=wi;j;k j 6 10

ð22Þ

where w represents the transported variables described in each of the conservation Eqs. (2), (3), (13), (14), and (19), and i, j, and k are the spatial indices for the x-, y-, and z-directions, respectively. 4. Results and discussion The operating in-house fuel cell described in Section 2 has been computationally modelled. The model was run at different cell potentials in order to generate the polarisation curve. For the given operating conditions, the modelling outputs are in very good agreement with the experimental data, see Fig. 4. 4.1. Sensitivity to anisotropic permeability and electrical conductivity As the main objective of this work is to investigate the sensitivity of the fuel cell performance to the anisotropy in the gas permeability and the electrical conductivity of the GDL, the computation cases shown in Table 2 were simulated. It should be noted that the through- and in-plane permeability and conductivity values listed in Case 1 in the table were experimentally estimated as described in [12–15]. They were fed into the model used to generate the

56

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63

Fig. 4. Experimental and modelling polarisation curves for 10BC-MEA.

Fig. 5. Polarisation curves for the computation cases listed in Table 2. It should be noted that the polarisation curves for Cases 1, 2, and 3 overlap each other.

numerical polarisation curve shown in Fig. 4. The effect of the gas permeability on the fuel cell performance was investigated by running Cases 2 and 3. The permeability of the GDLs was assumed to be isotopic and having either the value of the characterised through-plane (Case 2) or in-plane (Case 3) permeability. On the other hand, the effect of the electrical conductivity of the GDLs was investigated by running Cases 4 and 5. Again, the electrical conductivity was assumed to be isotropic and having either the value of the characterised through-plane (Case 4) or in-plane (Case 5) conductivity. It should be noted that the transport properties are not independent of each other and are related to the microstructure and the properties of the material. Therefore it is hypothetical that the gas permeability of the GDL changes while its electrical conductivity remains constant (as is the situation in Cases 2 and 3) and vice versa (as is the situation in Cases 4 and 5). Such an assumption has been made in order to only investigate the effects of the variable under investigation, e.g. the gas permeability of the GDL in Cases 2 and 3. Fig. 5 shows the corresponding polarisation curves for the Cases 1–5. For a realistic range of permeability and electrical conductivity values, it can be seen that the fuel cell performance does not change with the gas permeability of the GDL. On the other hand, it shows a significant sensitivity to the electrical conductivity of the GDL. To illustrate, the average current density at 0.55 V is underestimated (compared to that of the base case) by about 23% if the electrical conductivity is assumed to be isotropic and having the value of the through-plane conductivity (Case 4), i.e. 48 S m1. In contrast, the same average current density is overestimated by about 30% (compared to that of the base case) if the electrical conductivity is assumed to be isotropic and having the value of the inplane conductivity (Case 5), i.e. 4000 S m1. The distribution of the current density (at 0.55 V) within the cathode GDL, halfway the

length of the channel, was generated for all the cases listed in Table 2, see Fig. 6. For a cell potential of 0.55 V, it can be seen from Fig. 6 that the distributions of the local current density in all the simulated cases have a similar trend; the local current density is a minimum over the central region of the channel and increases towards the rib of the collector. Also, the maximum current densities are generally located in the vicinity of the edge between the rib and the flow channel. Except for Case 4, where the electrical resistance of the GDL is considerably higher, this indicates that the distribution of the current density is not dominantly controlled by the ohmic losses; also it is affected by the concentration losses. This is manifested by the maximum current densities being away from the central region of the ribs; they are located where the combined supply of charge and reactant to the catalyst layer is optimum. Fig. 6a shows that the permeability has virtually no effect on the distribution of the current density. For the given set of realistic and relatively low permeability values, it is evident that convection, which scales with permeability, has a minimal contribution to the transport in the through-plane direction (i.e. from the flow channel to the catalyst layer). In contrast, the electrical conductivity plays an influential role in determining the amount and the distribution of the current density, see Fig. 6b. In general, the average local current increases with increasing electrical conductivity of the GDL. The increased electrical conductivity reduces the ohmic losses, and consequently, the solid potential of the GDL is decreased. This, according to Eq. (11), increases the activation overpotential which in turn increases the amount of current being consumed in the catalyst layer; this is evident from Butler–Volmer Equation, Eq. (10). For a given through-plane conductivity of 48 S m1, if the inplane conductivity is decreased from 4000 S m1 (Case 1) to

Table 2 List of the computation cases investigated. Case number

1 2 3 4 5

(Base Case) (K= = K\) (K\ = K=) (r= = r\) (r\ = r=)

Permeability (m2)

Electrical conductivity (S m1)

Through-plane

In-plane

Through-plane

In-plane

4.97  1013 4.97  1013 1.87  1012 4.97  1013 4.97  1013

1.87  1012 4.97  1013 1.87  1012 1.87  1012 1.87  1012

48 48 48 48 4000

4000a 4000 4000 48 4000

a It is the arithmetic mean for the measured in-plane conductivities for SGL 10BC, namely 3000 and 5000 S m1 [13]. The validity of this simplification will be elaborated in Section 4.2.

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63

a

a

b

b

Fig. 6. Distribution of the local current density within the cathode GDL for (a) Cases 1, 2 and 3, and (b) Cases 1, 4 and 5. Straight lines in (b) represent the average local current densities.

48 S m1 (Case 4), the current over the flow channel is decreased. This is due to the increased in-plane resistance over the flow channel. In the same context, the maximum current density moves to be over the rib since the corresponding paths for transmitting electrons to the catalyst layer are shorter than that over the channel. For Case 4, where the performance of the fuel cell is dominantly controlled by the ohmic losses, these paths offer the least electrical resistance. On the other hand, for a given in-plane conductivity of 4000 S m1, if the through-plane conductivity is increased from 48 S m1 (Case 1) to 4000 S m1 (Case 5), the current over the rib is substantially increased and this is due to the significant decrease in the through-plane electrical resistance over the rib. Also, due to the enhanced through-plane conductance, the maximum current density moves more towards the rib of the collector where the path to transmit electrons is shorter. Interestingly, the distribution of the current density is more uniform in Case 1. The reason behind this is that (i) the throughplane electrical conductivity is not very high to bring the current density up over the rib, and (ii) the in-plane conductivity is not very low to bring the current density down over the channel. It is worth noting that, for Cases 1 and 5, the current decreases towards the centre of the rib where the concentration of oxygen is a minimum, see Fig. 7b. In Case 4, this phenomena is less profound since the increased ohmic resistance dictates a slow oxygen reduction reaction in the catalyst layer and subsequently lower consumption of oxygen. This translates into presence of sufficient

57

Fig. 7. (a) Distribution of activation overpotential (absolute value) within the cathode catalyst layer for Cases 1, 4 and 5 and (b) distribution of oxygen concentration within the cathode GDL for Cases 1, 4, and 5.

Fig. 8. Distribution of temperature within cathode catalyst layer for Cases 1, 4, and 5.

oxygen in the adjacent GDL in such a way that the effects of concentrations losses become minimal. Fig. 7 presents the distributions of activation overpotential and oxygen concentration within the cathode catalyst layer and the GDL, respectively, for Cases 1, 4 and 5. It can be seen from Fig. 7a that the activation overpotential is more uniform for Cases 1 and

58

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63

Fig. 9. Distribution of diffusive flux within cathode GDL for Cases 1, 4 and 5.

5 compared to that for Case 4. This is basically due to the reduced ohmic losses over the channel, as a result of the enhanced electronic conduction in the in-plane direction. Notably, the activation

overpotential is higher over the ribs than over the flow channels for all the simulated cases. The electrons take longer paths over the flow channels to reach the catalyst layer than over the ribs. This leads to higher electrical resistance, higher solid potential, and subsequently lower activation overpotential over the flow channel. Fig. 7b shows that the concentration of oxygen in the GDL is the lowest for Case 5 which is in accordance with the highest current density demonstrated by the same case, see Fig. 6b. This implies that, relative to Cases 1 and 4, the electrochemical reaction is fast and consequently the consumption of oxygen in the catalyst layer is high. Ultimately, this translates into the presence of less amount of oxygen in the adjacent GDL. As expected, the concentration of oxygen over the second channel is lower than that over the first channel. This can be explained by the fact that, according to the oxygen reduction reaction, the oxygen and water are consumed and produced respectively as the gas mixture moves towards the outlet. In contrast, the activation overpotential is higher over the second channel than over the first channel. This can be explained as follows. According to the Tafel Equation, gact = b ln (i/io), which is a valid approximation for Butler–Volmer Equation at the cathode side [21], the activation overpotential increases with decreasing exchange current density, io (b is the Tafel slope in the above equation). The exchange current

a

b

Inlet

Outlet Fig. 10. (a) Velocity vector plot within the cathode GDL, and (b) pressure distribution in the cathode side.

59

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63

Fig. 11. Distribution of convective flux within cathode GDL for Cases 1, 2 and 3.

density itself increases exponentially with temperature [44]. Fig. 8 shows the distribution of temperature for Cases 1, 4, and 5. Although it is not clear, the temperature profile over the second channel is slightly lower than that over the first channel; this is most likely due to the decreased rate of the exothermic reaction and subsequently the decreased amount of heat generated. Based on the above, the effects of such a slight decrease in temperature are amplified to yield a decreased exchange current density and increased activation overpotential over the second channel. As expected, the temperature profile is a maximum for Case 5, where the rate of the exothermic reaction (and subsequently the rate of heat generated) is a maximum; this is inferred from the maximum average current density presented by this case. Fig. 9 shows the profile of diffusive flux within the cathode GDL for Cases 1, 4 and 5. The diffusive flux demonstrates a considerable change with the electrical conductivity of the GDL. The diffusive flux is a maximum for Case 5. This is again in line with the maximum average current density (i.e. the maximum rate of reaction) shown by this case, see Fig. 6b. Simply, the higher is the rate of reaction, the higher is the rate of oxygen consumption in the catalyst layer, the higher is the concentration gradient, and the higher is the diffusive flux. Since the gas mixture is supplied through the flow channel, the diffusive flux of oxygen is more profound over the channel than over the rib. It can be also seen that the diffusive flux increases towards the edge between the rib and the flow channel. As mentioned earlier in this section, the rate of reaction is a maximum in the vicinity of the above edge where the supply of both oxygen and charge is optimum. Such a high rate of reaction dictates a high supply of oxygen via the main mode of transport in the through-plane direction, diffusion. On the other hand, as expected, the diffusive flux was found to be insensitive to changes in the permeability of the GDL (not shown). As for the convective flux, Fig. 10a shows that this type of transport is important in the regions over the rib since the pressure drop present between the inflow channel and the outflow channel (see

Fig. 10b) forces the gas mixture to flow in the plane of the GDL. As expected, Fig. 11 shows that the convective flux decreases with decreasing permeability – it is a minimum when the permeability of the GDL is 4.97  1013 m2. This is in agreement with the finding of Pharaoh’s computational study [25] that shows that almost no convective flow passes through the plane of the GDL if the permeability of the latter is of the order of, or less than, 1  1013 m2. It is worth noting that the convective flux was found to be insensitive to changes in the electrical conductivity of the GDL (not shown). It must be stressed that the parametric study was deliberately conducted at an intermediate cell potential (e.g. 0.55 V) as it is normally the cell potential at which most of the real-life PEM fuel cells are operated. However, the above distributions were also generated at a higher cell potential, i.e. 0.78 V, and they were found to maintain the same relations. At low cell potentials, the concentrations losses are expected to prevail and consequently the maximum current density is likely to be more towards the central region of the flow channel, where the amount of the reactant is a maximum. Having stated this, the rationales used in this study are valid for all the low- and high-cell potential cases. 4.2. Anisotropic in-plane electrical conductivity It has been experimentally demonstrated in Ref. [13] that there is anisotropy in the in-plane electrical conductivity. Namely, the electrical conductivity in the in-plane direction in which most of the fibres are oriented (0° direction) is larger than that in which the fibres are less oriented (90° direction) by a factor of about 2. To exemplify, for SGL 10BC GDL (which has been employed in the MEA being investigated in the current study), the in-plane electrical conductivity was measured in the 0° and 90° directions to be about 5000 and 3000 S m1, respectively [13]. For the sake of simplicity, the in-plane electrical conductivity in Case 1 shown in Table 2 has been assumed to be isotropic and has a value of 4000 S m1, which is the arithmetic mean of the measured in-plane electrical conductivies for SGL 10BC. In order to investigate the sensitivity of the fuel cell performance to the anisotropy in the in-plane electrical conductivity, the cases listed in Table 3 were run at a cell potential of 0.55 V and the corresponding current densities were obtained as shown in the fourth column of the table. It is clear that the resulting current density changes very slightly from case to case. In general, for a realistic range of values for the in-plane electrical conductivity, the fuel cell performance becomes slightly better with an increase in the in-plane electrical conductivity. For the above computational cases, Fig. 12 shows the distribution of the current density within the cathode GDL and the activation overpotential within the cathode catalyst layer respectively. The local current density changes very slightly from case to case, although they maintain the same trends. Such a slight change is attributed, as explained in the last section, to the changes of activation overpotential with the electrical conductivity of GDL; the higher is the electrical conductivity, the higher is the activation overpotential. Practically, it could be assumed that, for the given characterised in-plane electrical conductivities, the in-plane conductivity is isotropic.

Table 3 List of the computation cases for the in-plane electrical conductivity of the tested GDL.

a

Case

Conductivity in 0° in-plane direction (S m1)

Conductivity in 90° in-plane direction (S m1)

Resulting current density at 0.55 V (mA cm2)

A (Base Case) B ar¼ ¼ ðr¼;0 þ r¼;90 Þ=2 C (r¼;90 ¼ r¼;0 ) D (r¼;0 ¼ r¼;90 )

5000 4000 5000 3000

3000 4000 5000 3000

212.77 214.15 215.69 212.65

Case 1 in Table 2.

60

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63

4.3. Effect of U-bend Fig. 10a shows that there are recirculations in the region of the cathode GDL laying over the U-bend of the channel and this is

due to the deceleration of the flow in the U-bend. The question here is whether the performance of the fuel cell would be affected by such phenomena. As mentioned in Section 4.1, the main mode of transport for reactant gases from the flow channel to the

a

b

Fig. 12. (a) Distribution of local current density within cathode GDL and (b) distribution of activation overpotential (absolute value) within cathode catalyst layer.

Fig. 13. Vector plot for the oxygen diffusive flux inside the cathode GDL. The encircled area represents the region of the GDL which is located over the bend.

Fig. 14. Straight channel-based computational domain.

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63





 Fig. 15. Polarisation curves of a serpentine channel – and straight channel – based models for the in-house PEM fuel cell.

 catalyst layers is diffusion. It can be clearly seen from Fig. 13 that the diffusive flux in the region over the bend are of a similar intensity to that of the fluxes available elsewhere in the cathode GDL – this signifies that the diffusive flux has not been virtually influenced by the recirculations present in the bend of the channel. Therefore, one may deduce that such recirculations would not affect the performance of the fuel cell (but they may increase the pressure drop along the channel as shown, for example, by Ref. [45] and this ultimately leads to the use of high pumping power). If this is the situation, a straight channel-based model could satisfactorily represent the operating fuel cell and it will enable the user to quickly obtain the global performance of the fuel cell. Fig. 14 presents the computational domain of the straight channel. Upon checking for mesh independence and being simulated at different cell potentials, Fig. 15 shows that the polarisation curve of the straight channel-based model is more or less identical to that of the serpentine channel-based model. Accordingly, such a considerably less computationally expensive model could be used if the intention is only to obtain the global performance of the fuel cell. However, if the convection from one turn to another is to be captured in the model, then the serpentine channel-based model must be used. Further, it is most likely to be inappropriate to use the straight channel-based model if the fuel cell is operated under multiphase operating conditions. Under these operating conditions, the recirculation and separation regions may cause a performance loss with a local accumulation of liquid water in the bend [46].

61

GDL over the rib of the current collector (due to the presence of pressure drop between the inflow and outflow channels), and (ii), as expected, increases with an increase in the permeability of the GDL. In contrast, the fuel cell performance was found to change significantly with the electrical conductivity of the GDL. The greater is the electrical conductivity, the better is the performance of the fuel cell. This is due to the fact that the activation overpotential, with which the current density scales, increases with an increase in the electrical conductivity. For the reported cell potential, although it tends to be more towards the regions over the rib of the current collector, the maximum local current density within the GDL is neither over the middle of the flow channel nor the middle of the rib. This indicates that the fuel cell performance is not entirely dominated by the ohmic losses in the intermediate current density region. For the given experimental data, the fuel cell performance is almost insensitive to the anisotropy in the in-plane electrical conductivity and therefore, in practice, such anisotropy can be ignored. For single phase operating conditions, a straight channel section of the operating fuel cell can be modelled and used as a quick tool to generate the global performance curves. If the cross-convection is to be captured and/or the fuel cell is operated under multi-phase operating conditions, the serpentine channel– based model should be used.

Recent studies [47–49] have shown that the transport properties may locally change by orders of magnitude as a result of the non-uniform compression that a GDL typically experiences inside the fuel cell. Therefore an interesting future research investigation would be the effects of non-uniform compression on the performance of the fuel cell. Acknowledgments The authors gratefully acknowledge financial support from a Dorothy Hodgkin Postgraduate Award, the UK Engineering and Physical Sciences Research Council (EPSRC), and Shell UK. The authors are grateful to ESI-CFD for providing licences of CFDACE+. Also, the authors would like to thank SGL Technologies GmbH, Germany, for providing the GDL materials. Appendix A In the CFD-ACE solver, the following equations have been used to calculate the physical properties of the gases [40]:

5. Conclusions

A.1. Density

In this paper, a parametric study has been conducted in order to investigate the effects of the anisotropic experimentally-determined permeability and electrical conductivity of the GDL on the fuel cell performance. For the reported range of the permeability and electrical conductivity, the main conclusions can be summarised as follows:

When compressible effects are not negligible, the Ideal Gas Law should be used. The Ideal Gas Law is given by:

 The performance of the fuel cell was found to be insensitive to the anisotropy in the permeability of the GDL. The sufficiently low values for the given realistic permeability impose a minimal contribution of convection towards the mass transport in the through-plane direction, i.e. from the flow channel to the catalyst layer (where the reaction takes place). However, convection was found to (i) be significant in the regions of the

q¼

ðp þ pref ÞM RT

ðA:1Þ

where pref is the reference pressure, p is the calculated static pressure, M is the molecular weight of the gaseous mixture, R is the universal gas constant, and T is the temperature. A.2. Viscosity The Mix Kinetic Theory option in CFD-ACE+ uses the kinetic theory of gases to calculate the viscosity of the gas or mixture of gases. For a pure monatomic gas, the viscosity is defined as follows:

62

li

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63

pffiffiffiffiffiffiffiffiffi Mi T ¼ 2:6693  105 2

ðA:2Þ

ci Xl

where li is the dynamic viscosity of species i, Mi is the molecular weight of species i, T is the temperature in Kelvin, ci is the characteristic diameter of the molecule in Angstroms, and Xl is the collision integral. The latter is given by

Xl ¼

1:16145 0:52487 2:16178 þ þ e0:7732T€ e2:43787T€ T€ 0:14874

ðA:3Þ

where T€ is the dimensionless temperature and is given by:

kB T T€ ¼ d

ðA:4Þ

where d is the characteristic energy, and kB is the Boltzmann constant. To calculate the mixture viscosity using kinetic theory, the following equation is used.

lmix ¼

N X i¼1

Yl Pi i Y j /ij

ðA:5Þ

where Yi and Yj are the mass fractions of species i and j, respectively, and /ij is a dimensionless quantity, which is given by:

2  1=2 1 Mi 41 þ /ij ¼ pffiffiffi 1 þ Mj 8

li lj

!1=2   32 1=4 Mj 5 Mi

ðA:6Þ

A.3. Specific heat capacity The mix JANNAF method is used to calculate the specific heat capacity, CP, of the gases:

CP ¼ a1 þ a2 T þ a3 T 2 þ a4 T 3 þ a5 T 4 þ a6 T 5 R

ðA:7Þ

The coefficients a1  a6 are obtained by curve fitting the experimental data. A.4. Thermal conductivity If the Prandtl number is given (as is in our case, Pr = 0.707), the CFD-ACE solver uses this value to calculate the thermal conductivity as follows:

j¼

CP l Pr

ðA:8Þ

where j is the thermal conductivity, CP is the specific heat capacity, and l is the dynamic viscosity. A.5. Diffusion coefficient If the Schmidt number is given (as is in our case, Sc = 0.7), the mass diffusivity of a gas is given by:

D¼

l qSc

ðA:9Þ

References [1] Wang Y, Wang CY. Simulation of flow and transport phenomena in a polymer electrolyte fuel cell under low-humidity operation. J Power Sources 2005;147:148–61. [2] Zhou TH, Liu HT. A 3D model for PEM fuel cells operated on reformate. J Power Sources 2004;138:101–10. [3] Mazumder S, Cole JV. Rigorous 3-d mathematical modeling of PEM fuel cells – II. Model predictions with liquid water transport. J Electrochem Soc 2003;150:A1510–7. [4] Baschuk JJ, Li X. A comprehensive, consistent and systematic mathematical model of PEM fuel cells. Appl Energy 2009;86:181–93.

[5] Nam J, Chippar P, Kim W, Ju H. Numerical analysis of gas crossover effects in polymer electrolyte fuel cells (PEFCs). Appl Energy 2010;87:3699–709. [6] Perng SW, Wu HW. Non-isothermal transport phenomenon and cell performance of a cathodic PEM fuel cell with a baffle plate in a tapered channel. Appl Energy 2011;88:52–67. [7] Perng SW, Wu HW. Effect of the prominent catalyst layer surface on reactant gas transport and cell performance at the cathodic side of a PEMFC. Appl Energy 2010;87:1386–99. [8] Perng SW, Wu HW, Jue TC, Cheng KC. Numerical predictions of a PEM fuel cell performance enhancement by a rectangular cylinder installed transversely in the flow channel. Appl Energy 2009;86:1541–54. [9] Wu H, Berg P, Li X. Steady and unsteady 3D non-isothermal modeling of PEM fuel cells with the effect of non-equilibrium phase transfer. Appl Energy 2010;87:2778–84. [10] Wu HW, Ku HW. The optimal parameters estimation for rectangular cylinders installed transversely in the flow channel of PEMFC from a threedimensional PEMFC model and the Taguchi method. Appl Energy 2011; 88:4879–90. [11] Lim C, Wang CY. Effects of hydrophobic polymer content in GDL on power performance of a PEM fuel cell. Electrochim Acta 2004;49:4149–56. [12] Ismail MS, Damjanovic T, Ingham DB, Ma L, Pourkashanian M. Effect of polytetrafluoroethylene-treatment and microporous layer-coating on the inplane permeability of gas diffusion layers used in proton exchange membrane fuel cells. J Power Sources 2010;195:6619–28. [13] Ismail MS, Damjanovic T, Ingham DB, Pourkashanian M, Westwood A. Effect of polytetrafluoroethylene-treatment and microporous layer-coating on the electrical conductivity of gas diffusion layers used in proton exchange membrane fuel cells. J Power Sources 2010;195:2700–8. [14] Ismail MS, Borman D, Damjanovic T, Ingham DB, Pourkashanian M. On the through-plane permeability of microporous layer-coated gas diffusion layers used in proton exchange membrane fuel cells. Int J Hydrogen Energy 2011;36:10392–402. [15] Ismail MS, Damjanovic T, Hughes K, Ingham DB, Ma L, Pourkashanian M, et al. Through-plane permeability for untreated and PTFE-treated gas diffusion layers in proton exchange membrane fuel cells. J Fuel Cell Sci Technol 2010;7:051016–7. [16] Becker J, Fluckiger R, Reum M, Buchi FN, Marone F, Stampanoni M. Determination of material properties of gas diffusion layers: experiments and simulations using phase contrast tomographic microscopy. J Electrochem Soc 2009;156:B1175–81. [17] Zamel N, Li X, Shen J, Becker J, Wiegmann A. Estimating effective thermal conductivity in carbon paper diffusion media. Chem Eng Sci 2010;65:3994–4006. [18] Gostick JT, Fowler MW, Pritzker MD, Ioannidis MA, Behra LM. In-plane and through-plane gas permeability of carbon fiber electrode backing layers. J Power Sources 2006;162:228–38. [19] Bapat CJ, Thynell ST. Effect of anisotropic electrical resistivity of gas diffusion layers (GDLs) on current density and temperature distribution in a Polymer Electrolyte Membrane (PEM) fuel cell. J Power Sources 2008;185: 428–32. [20] Bapat CJ, Thynell ST. Effect of anisotropic thermal conductivity of the GDL and current collector rib width on two-phase transport in a PEM fuel cell. J Power Sources 2008;179:240–51. [21] Meng H, Wang CY. Electron transport in PEFCs. J Electrochem Soc 2004;151:A358–67. [22] Pharoah JG, Karan K, Sun W. On effective transport coefficients in PEM fuel cell electrodes: anisotropy of the porous transport layers. J Power Sources 2006;161:214–24. [23] Zhou T, Liu H. Effects of the electrical resistances of the GDL in a PEM fuel cell. J Power Sources 2006;161:444–53. [24] Zhang X, Song DT, Wang QP, Huang C, Liu ZS. Influence of anisotropic transport properties of the GDL on the performance of PEMFCs. In: Fuller T, Shinohara K, Ramani V, Shirvanian P, Uchida H, Cleghorn S, et al., editors. Proton exchange membrane fuel cells. ECS Trans, vol. 8; 2008. p. 913–23. [25] Pharoah JG. On the permeability of gas diffusion media used in PEM fuel cells. J Power Sources 2005;144:77–82. [26] Barbir F. PEM fuel cells: theory and practice. 1st ed. Oxford: Elsevier Academic Press; 2005 [Chapter 7]. [27] Dagan G. Flow and transport in porous formations. Berlin: Springer-Verlag; 1989. [28] Wang CY. Fundamental models for fuel cell engineering. Chem Rev 2004;104:4727–65. [29] Carcadea E, Ene H, Ingham DB, Lazar R, Ma L, Pourkashanian M, et al. A computational fluid dynamics analysis of a PEM fuel cell system for power generation. J Numer Methods Heat Fluid Flow 2007;17:302–12. [30] Ma L, Ingham DB, Pourkashanian M, Carcadea E. Review of the computational fluid dynamics modeling of fuel cells. J Fuel Cell Sci Technol 2005;2:246–57. [31] Zamel N, Li X, Shen G. Correlation for the effective gas diffusion coefficient in carbon paper diffusion media. J Energy Fuels 2009;23:6070–8. [32] Tomadakis MM, Sotirchos SV. Transport properties of random arrays of freely overlapping cylinders with various orientation distributions. J Chem Phys 1993;98:616–26. [33] Barbir F. PEM fuel cells: theory and practice. 1st ed. Oxford: Elsevier Academic Press; 2005 [Chapter 4]. [34] Barbir F. PEM fuel cells: theory and practice. 1st ed. Oxford: Elsevier Academic Press; 2005 [Chapter 2].

M.S. Ismail et al. / Applied Energy 95 (2012) 50–63 [35] Springer TE, Zawodzinski TA, Gottesfeld S. Polymer electrolyte fuel-cell model. J Electrochem Soc 1991;138:2334–42. [36] Sadeghi E, Bahrami M, Djilali N. Analytic determination of the effective thermal conductivity of PEM fuel cell gas diffusion layers. J Power Sources 2008;179:200–8. [37] Bernardi DM, Verbrugge MW. A mathematical-model of the solid-polymerelectrolyte fuel-cell. J Electrochem Soc 1992;139:2477–91. [38] Kim GS, Sui PC, Shah AA, Djilali N. Reduced-dimensional models for straightchannel proton exchange membrane fuel cells. J Power Sources 2010;195:3240–9. [39] Cheng CH, Lin HH, Lai GJ. Numerical prediction of the effect of catalyst layer Nafion loading on the performance of PEM fuel cells. J Power Sources 2007;164:730–41. [40] CFD-ACE+ User’s Manual, ESI-CFD, Huntsville, AL; 2009. [41] Bac2 Company. . [42] Srinivasan S, Ticianelli EA, Derouin CR, Redondo A. Advances in solid polymer electrolyte fuel-cell technology with low platinum loading electrodes. J Power Sources 1988;22:359–75.

63

[43] Larmine J, Dicks A. Fuel cell systems explained. 2nd ed. Chichester: John Wiley and Sons; 2005 [Chapter 3]. [44] Barbir F. PEM fuel cells: theory and practice. 1st ed. Oxford: Elsevier Academic Press; 2005 [Chapter 3]. [45] Maharudrayya S, Jayanti S, Deshpande AP. Pressure losses in laminar flow through serpentine channels in fuel cell stacks. J Power Sources 2004;138:1–13. [46] Yoon SY, Ross JW, Mench MM, Sharp KV. Gas-phase particle image velocimetry (PIV) for application to the design of fuel cell reactant flow channels. J Power Sources 2006;160:1017–25. [47] Tamayol A, McGregor F, Bahrami M. Single phase through-plane permeability of carbon paper gas diffusion layers. J Power Sources 2012;204:94–9. [48] Tamayol A, Bahrami M. In-plane gas permeability of proton exchange membrane fuel cell gas diffusion layers. J Power Sources 2011;196:3559–64. [49] Jiao K, Park J, Li X. Experimental investigations on liquid water removal from the gas diffusion layer by reactant flow in a PEM fuel cell. Appl Energy 2010;87:2770–7.