Abstract
A one-dimensional model of a high temperature polymer electrolyte membrane fuel cell using polybenzimidazole (PBI) membranes is described. The model considers mass transport through a thin film electrolyte covering the catalyst particles as well as through the porous media. The incorporation of a thin film model describing reactant gas mass transport through electrolyte covering the electrocatalyst is shown to be an essential requirement for accurate simulation. The catalyst interface is represented using a macrohomogeneous model. The influence of carbon monoxide, carbon dioxide, and methane, which would be present in a reformate gas, is considered in terms of the effect on the anode polarisation/kinetics behaviour. The model simulates the influence of operating conditions, cell parameters, and fuel gas compositions on the cell voltage current density characteristics. The model gives good predictions of the effect of oxygen and air pressures on cell behaviour and correctly simulates the mass transport behaviour of the cell. The model with reformate gas shows that additional voltage losses associated with CO poisoning can lead to loss in voltage of tens of mV and thus reduction in power.
1. Introduction
There are numerous mathematical models of polymer electrolyte membrane fuel cells (PEMFCs) that have been reported in the literature. These include empirical (curve fitting) and zero-dimensional models which essentially coupled thermodynamic, kinetic, and resistance effects, to one/two-dimensional phenomenological models. One of the early phenomenological models of a PEMFC with a Nafion membrane was developed by Bernardi and Verbrugge [1]. Since then, significant developments have been made. Models were used to study concentration and current distributions in PEMFCs [2], to map liquid saturation and temperature distributions [3] for comparison with experimental data, to solve equations describing multicomponent flow in diffusion layers and flow channels [4], or to module mass transport in porous electrodes [5].
The mass transport behaviour in PEMFC electrodes was considered to be either in a single phase [6], where no liquid water is formed, or involving two-phases [7] involving an electrolyte/water film. It was found [8, 9] that under the assumption of no liquid water formation, the model consistently overpredicted measured polarization behaviour.
On the other hand, results showed that the inclusion of liquid water transport greatly enhanced the predictive capability of the model and was necessary to match experimental data at high current density [10]. This could be achieved in the form of one-dimensional [11], two-dimensional [12] and three-dimensional models [13, 14].
In comparison, there has been little modeling of PEMFCs based on polybenzimidazole (PBI) membranes. The first proposed model [15] was a parametric model, in which a very low value of cathode porosity of 9.2% and a high value for the transfer coefficient (α) of 2 were used to fit the model to experimental data. However such a model could not explain the limiting current observed under air operation [15].
A second one-dimensional model assumed the Tafel approximation to describe the electrode kinetics, with a transfer coefficient equal to 0.5 whilst the exchange current density was fitted to the polarisation curves using a least square error method leading to a reaction order for oxygen equal to 0.7 [16]. The simulated polarisation curve showed a better fit for air than for oxygen, where the model underestimated the performance at the higher current densities. The influence of humidity and consequently product water generation on membrane conductivity was given as a reason for the observed behaviour (the model assumed constant membrane conductivity).
A third model considered a three-dimensional structure, with transfer coefficient (α) equal to 2 and reaction order () equal to 1. However, once again the model failed to explain the observed difference between air and oxygen operation: with oxygen, simulation underestimated the experimental data, whilst with air operation, the simulated data overpredicted the experimental data at high current densities, where contrary to experimental results no limiting current was observed [17].
Scott et al. also proposed a one-dimensional model for PBI-based fuel cells. They described electrode kinetics by the Butler-Volmer equation and mass transport by the multicomponent Stefan Maxwell equations coupled with Darcy’s law. The model had a good fit with the experimental data but failed to show limiting current under air operation [18].
Similarly, one-and two-dimensional degradation models were constructed to simulate the steady state polarisation curves recorded at different times during aging test.
The models again failed to show any apparent mass transport limitations under air operation (limiting current) and used α equal to 1 and equal to 1 [19, 20].
The failure of the reported models to predict the mass transport limitations under air operation and therefore over estimating cell performance, particularly at high current densities, was caused by the assumption that mass transport solely occurred through the porous media. This is similar to the single phase mass transport observation discussed earlier in Nafion based PEMFC models. In reality an electrolyte (PBI/Acid) thin film surrounding the catalyst sites (particles) is present and mass transport through this phase should be considered. In this film reactants have to dissolve in the electrolyte media and diffuse through it to reach the catalytic sites, in a similar way to a two phase mass transport approach where diffusion through liquid water is considered. Diffusion through the acid electrolyte is much slower in comparison to that through porous media and can explain the observed PBI-PEMFC mass transport behaviour. The effect of the electrolyte thin film has been realised and modelled in phosphoric acid fuel cells [21].
Lobato et al. presented a model of a PBI-based high temperature PEMFC using neural networks [22]. The model can be used to elucidate the processes within the cells, by allowing optimization of materials, cells, stacks, and systems and support control systems.
This absence of a thin film model approach also explains the relatively thick catalyst layers and very low porosity used in previous reported models [15, 16] in attempts to compensate for the thin film effects and to try to match the experimental data. Similarly, whilst most models used a value for the reaction order equal to one, values varied from 0.5 to 2. As shown elsewhere [23], values of transfer coefficient changed with doping level. For example, at 150°C, the lowest Tafel slope observed at the minimum doping level of 4.5 PRU (4.5 mole per repeat PBI unit) was 92 mV dec−1() which increased to 104 mV dec−1() at a doping level of 10 PRU [24, 25] and to 90–135 mV dec−1() for phosphoric acid fuel cells [26–32].
Additionally, empirical models have been developed to study temperature effects in PBI-based PEMFCs. A change in the transfer coefficient (increase) with temperature was obtained using a nonlinear, least squares fitting method [33]. This result agrees with the finding of others [23] and similar documented behaviour in PAFC [26, 28, 29, 34].
Modelling of PBI PEMFCs will increase understanding of their behaviour, enable prediction of their performance, and assist with their operational control. Thus, a one-dimensional (1D) model of the high temperature PBI-based PEMFC was developed which considered multicomponent mass transport of a reformate-based fuel, through porous media and a thin electrolyte film covering the electrocatalyst in the electrode layers. The model is a simplification of a recent two-dimensional model [35] which was solved using commercial software. The simplified model was designed to be used as a platform for a high temperature PEMFC control system. In this paper we also demonstrate that by ignoring the influence of mass transport through a thin film electrolyte over the catalyst, the PBI PEMFC model can fail to give accurate simulations of fuel cell voltage performance.
2. Experimental
Gas diffusion electrodes suitable for high temperature operation (nonwoven carbon cloth) incorporated with wet proofed microporous layer (obtained from Freudenberg (FFCCT, Germany) were used as substrates to deposit the catalyst layer for both anode and cathode.
A spraying machine was built to deposit reproducible catalyst layers of uniform structure/porosity. A Computer numerical control (CNC) milling machine (Sherline 2010, USA) was modified and used to provide the desired spraying pattern, whilst a fixed stainless steel spraying (0.5 mm) nozzle (Schlick 970S8, Germany) and associated metering valve were used to control the spray mixture (with nitrogen) and catalyst ink flow rate.
The catalyst ink was prepared by sonicating the catalyst and PTFE dispersion (60% wt, Aldrich) in water-ethanol mixture, 50% Pt/C (Etek, USA), 0.4 mgPtcm−2 in the cathode and 20% Pt/C (Etek, USA), 0.2 mgPtcm−2 in the anode. The final PTFE content in the catalyst layer was 40 and 20% wt for cathode and anode, respectively. The required amount of phosphoric acid (PA) was then added to the surface by means of a micropipette, and the electrodes were left for a week to cure to achieve a uniform acid distribution. A loading of 2 and 4 mgH3PO4cm−2 was used for the cathode and anode, respectively.
The MEA was finally obtained by hot pressing the electrodes, onto a 5.6 mole% phosphoric acid doped PBI membrane at 150°C for 10 minutes with load of 40 kg cm−2.
The cell body had a 3 cm × 3 cm gold-plated parallel flow fields. Mica-filled PTFE inserts were used to surround the flow fields and provide location for the O-ring seal and a dynamic hydrogen electrode (DHE). The solid state DHE consisted of two platinum wires on each side of the membrane located outside the O-ring: a distance of 10 mm away from the MEA edge to avoid side current effects (the membrane used was ~50 μm thick). A small current of 1.0 mA cm−2 (~10 μA) was applied by means of 9.0 V battery connected in series with an appropriate resistance.
3. Mathematical Model of the Fuel Cell
The mathematical model of the fuel cell was one-dimensional, and the gas flow channels were not considered. The fuel cell consisted of two-diffusion layers, anode and cathode catalyst layers, and the membrane. The assumptions adopted in the model are as follows.(i)Steady state and isothermal operation.(ii)Mass transport was solely due to diffusion where convection effects are negligible.(iii)Ideal gas behaviour.(iv)Membrane was impermeable to hydrogen and oxygen.(v)Negligible contact resistances between components.(vi)No membrane swelling after installation in the test cell.(vii)Only gas phases present (no water condensation, temperatures were 120°C).(viii)Catalyst layer was thin and treated as uniform potential and concentration interface (i.e., 0D).(ix)Isotropic macrohomogeneous porous regions.
Isothermal operation was a reasonable assumption as the test cell temperature was controlled using electrical heating and an on-off controller. Ideal gas behaviour was appropriate, as the cell was not operated at high pressure. The cell was operated at relatively high temperatures, above the normal boiling point of water, and thus it was reasonable to assume that there was no liquid saturation and that single phase behaviour applied. The flow channel was not included in the model as the boundary conditions between the flow channel and the diffusion layer were taken as the feed gas compositions. This was justified by the high gas flow stoichiometry used (>2.2 at the highest current).
The macrohomogeneous model for the catalyst layer assumes that the porous electrode is an “average” of the solid electrode and the electrolyte (Figure 1). Thus, the effective conductance of the porous electrode is the weighted volume average of the respective conductance. Diffusion coefficients and other properties are similarly averaged.
The objective of the model is to determine the effect of a range of operating variables and parameters on the cell voltage of the PBI-based fuel cell. The overall cell voltage is given from a combination of the thermodynamic cell potential and voltage losses associated with Ohmic resistances in the electrodes and membrane, and kinetic losses at the anode and cathode which are influenced by mass transport restriction where is the reversible cell potential, refers to electrode polarisation losses, and is the Ohmic resistance loss.
3.1. Thermodynamic Equilibrium Potential
The overall electrochemical reaction in a PEM fuel cell running on H2 as fuel and O2 as oxidant at temperature above 100°C can be written as Considering that thin electrolyte film covers the catalyst surface, the partial pressure of the reactants/products should be replaced with their activity in the electrolyte. The activity of hydrogen and oxygen can be replaced by their concentration in the thin film considering their activity coefficient is close to 1 (the concentrations are very low).
The thermodynamic equilibrium potential can be calculated using the Nernst equation: where oxygen and hydrogen concentration (solubility) is obtained using Henry’s law: where (atm) is the equilibrium partial pressure of species above the electrolyte film and (atm cm3 mole−1) is Henry’s constant for the given species-electrolyte (H3PO4) pair at a given temperature . The term RT inside the logarithm expression in (3) was introduced to convert the concentration (mole cm−3) into pressure units (atm) to calculate the standard cell reversal potential.
, and are the reference oxygen, hydrogen partial pressure, and water activity, respectively. Their values are equal to unity.
was calculated using [36] where is in kJ is in J , and is in K.
The water activity in is given by [37] where is the water vapour pressure in equilibrium with the acid electrolyte, is the saturation vapour pressure of pure water at the same temperature.
PBI-based fuel cells can operate under dry conditions or very low humidification; in this study the gas reactants were passed through a humidifier at room temperature (18°C) prior to entering the cell ( at 150°C), the initial water vapour pressure (O.C.P conditions) was considered to be equal to the water saturation pressure at 18°C (291 K), that is, atm.
Saturated water vapour pressures were obtained from steam tables [36]. The following polynomial function presents this data in the temperature range of 273–500 K: where is in atm and is in Kelvin.
3.2. Gas Transport in Porous Media
3.2.1. Diffusion in the Porous Cathode
There were three species in the cathode gas stream, that is, oxygen, nitrogen and water.
For the purpose of this model only one-dimensional diffusion, normal to the face of the electrode, is considered. Diffusion of multicomponent gas streams through the porous carbon electrode can be described using the Stefan-Maxwell equation
is the mole fraction of species , is the mole flux of species , and is the effective binary diffusion coefficient for the pair - in the porous medium.
can be calculated using the corrected Slattery-Bird correlation [38] to account for the porosity/tortuosity effects with the Bruggeman correlation [18] where and are the gas critical temperature and pressure, respectively. is the gas molecular weight, is the porosity and is the tortuosity. and are constants, is 0.0002745 for diatomic gases and 0.000364 for water vapour and is 1.832 for diatomic gases and 2.334 for water vapour.
The species’ flux can be given as follows:
nitrogen (inert species):
oxygen:
water: where is the current density per geometric electrode area and F is Faraday’s constant.
Substituting the species flux in (8), we obtain
At the catalyst layer where is thickness of the gas diffusion electrode we can write [36]Oxygen molar fraction is given by
3.2.2. Diffusion in the Porous Anode
The reformate anode gas feed can consist of a mixture of CH4, CO2, , and H2. For the five component gas anode species, we can write balances as where ,, , and , Solution of the above equations leads to
3.3. Transport through Thin Film Electrolytes
The macrohomogeneous model for the catalyst layer assumed that the catalyst layer is an “average” of the solid electrode and the electrolyte. Thus, the effective conductance of the catalyst layer was the weighted volume average of the respective conductance: diffusion coefficients and film thickness are similarly averaged, and so forth.
Oxygen transport from the porous media to the catalyst active surface area occurred through a thin acid film covering the catalyst agglomerates. The film provides proton conductive paths from the catalyst active sites to the membrane. The average film thickness was estimated using the following equation:
is the total mass added of H3PO4 per unit area(loading), is the density, and is the surface area of carbon/platinum per unit area covered by the electrolyte.
Values of oxygen diffusion in phosphoric acid (98% wt) were reported to be cm2 s−1, an order of magnitude higher than that of doped PBI (doping level of 6PRU) cm2 s−1 at 150°C. Similar values of dissolved oxygen concentration was obtained for doped PBI (doping level of 6 PRU) mole cm−3 and for 95% wt phosphoric acid 0.5 mole cm−3 at 150°C and atmospheric pressure [24].
Liu et al. studied oxygen reduction at the platinum/phosphoric acid doped PBI interface, and found that oxygen diffusion was increased by increasing the volume fraction of amorphous (free) H3PO4 (i.e., doping level). They also suggested that the crystalline PBI regions were not involved in proton or oxygen transport [24, 25].
To determine the oxygen/hydrogen concentration at the catalyst surface, we can derive from Fick’s law for diffusion:
is the molar flux (geometric area), is the reactant concentration on the catalyst surface, and is the equilibrium reactant concentration in the acid film at the studied temperature. as mentioned earlier was the real platinum surface area (ESA) per unit area that was covered with an electrolyte film (also known as the roughness factor ).
Due to insufficient data on hydrogen solubility in phosphoric acid at high temperature, hydrogen solubility was considered similar to that of oxygen [39] at the same conditions (pressure, temperature and phosphoric acid concentration).
Similarly, for the hydrogen diffusion coefficient in phosphoric acid electrolyte we can write [40]:
Experimental data (not shown) showed that with air at atmospheric pressure, the product decreased slowly with temperature, whilst for the same electrode under air at 1 bar (relative) the product exhibited a maximum at 150°C. This effect was caused by variations in diffusion coefficient and Henry’s constant at a given temperature with changes in phosphoric acid concentration and consequently water partial pressure above the electrolyte film; in other words the concentration and viscosity of phosphoric acid at a given temperature depended on the humidity content.
Klinedinst et al. studied oxygen solubility and diffusivity in hot phosphoric acid [41]. They suggested that both gas diffusivities and solubilities exhibit exponential reciprocal temperature dependencies, that is,
and are pre-exponential factors (given in the Appendix), is the diffusion activation energy, and is the enthalpy of solution. Both and change with the concentration of phosphoric acid.
3.3.1. Diffusion, Temperature, and Phosphoric Acid Concentration
A second-order polynomial was fitted using a least square error technique to fit the data obtained by Klinedinst et al. [41], for the activation energy of oxygen in phosphoric acid at different acid weight concentrations (). The correlation equation is given by
The high values of the diffusion activation energy were assigned to the extensive hydrogen bond network and high viscosity of concentrated H3PO4 [42]. This also explained the increase in activation energy with increasing phosphoric acid concentration.
The calculated diffusion values were in the range of ( cm s−1) depending on the temperature and phosphoric acid concentration; values were in good agreement with values obtained in the literature [32, 41, 43, 44].
3.3.2. Solubility, Temperature, and Phosphoric Acid Concentration
A third-order polynomial was built to fit oxygen enthalpy of solution data obtained from [41]
The enthalpy of solution (, as for % wt) decreased with increase in phosphoric acid concentration until it reached negative values at 96% wt. This means that a smaller decrease in the solubility of oxygen in phosphoric acid with temperature will occur as the acid concentration increases from 85% wt to 95% wt; beyond this concentration a slow increase in solubility with temperature will occur.
and can be obtained from (26) and (25) for a given acid weight concentration. Substituting their values into (20) and (21) will finally lead to , the oxygen concentration on the catalyst surface required in the Butler-Volmer kinetic equation.
3.3.3. Phosphoric Acid Concentration, Temperature, and Water Vapour Pressure
MacDonald and Boyack [45] studied the density, conductivity and equilibrium water vapour pressure of concentrated phosphoric acid from room temperature to 170°C. Data from their work was used to determine phosphoric acid concentrations (wt %) at a given temperature, and water vapour pressure.
To obtain a good fit without using polynomials of very high order, they suggested expressing concentrations as mole percent , instead of (wt %) using the formula
At any given concentration, a linear relation was obtained between log (/mm Hg) the equilibrium water vapour pressure and (°C−1):
Two functions were built to correlate and with , % mole; values of and with the relevant functions are given in the appendix.
During simulation of the model in this work, and were continuously updated based on water vapour pressure (or relative humidity) and temperature. was obtained from using
From (28) and (29), it can be seen that (or humidity) is highly dependent on phosphoric acid concentration at the typical low humidity conditions for PBI-based PEMFCs gas feed. This means that phosphoric acid concentrations will vary greatly with the water produced by the fuel cell (logarithmic relation).
3.4. Kinetics
The Butler-Volmer equation was used to describe the kinetics at the anode and cathode where the subscripts and are for anode and cathode, respectively. α is the transfer coefficient, is the exchange current density at the studied conditions per Pt unit area given by (32) below.
Values of α can be obtained from the Tafel slope . It was shown earlier [23] that the value of α changes with doping level for example, at 150°C the lowest Tafel slope observed at the minimum doping level of 4.5 PRU was 92 mV dec−1() increasing to 104 mV dec−1() at a doping level of 10 PRU [24, 25].
Oxygen reduction in PBI-free electrodes environment (or for PBI electrodes with very high doping level) is the same as that in phosphoric acid, while for PBI-based electrodes, appropriate transfer coefficient should be chosen, depending on the doping level [23].
Various Tafel slopes have been reported for phosphoric acid fuel cells in the range of 90 and 135 mV dec−1() at 150°C [26–32]. An average value of alpha equal to 0.75 (112 mV dec−1) was considered in this work. The closest available experimental data for a PBI-free environment at 150°C is from half cells results with a high doping (16 PRU) PBI electrode at temperature of 140°C [23]; the reported alpha value was 0.75, which agrees with the chosen reported alpha’s value for phosphoric acid.
Whilst the discussed values for alpha were for a temperature of 150°C, alpha values vary with operating cell temperature (Table 1). Appleby [26] reported variations of α from 0.56 at 25°C to 0.66 at 136°C in 85% wt phosphoric acid. Similarly, O’Grady et al. [34] reported values from 0.53 at 25°C to 0.68 at 70°C in 85% wt phosphoric acid. Haung et al. [27] reported α value (85% wt H3PO4) of 0.47, 0.61, and 0.67 at temperatures of 25, 100, and 150°C, respectively. Kunz and Gruver [46] reported at 160°C (96% wt H3PO4) using Pt/C as catalyst.
The variation of dependency with temperature was expressed as [48] where and are constants. Values for , the rate change of alpha with temperature, were 0.0014 [28], 0.0015 [47], 0.0034 [34], and 0.0043 [23] (PBI 16 PRU).
The value of alpha and its variation with temperature depend on the catalyst treatment and the impurity content in the acid [29, 47].
The exchange current density was obtained using (A cmPt−2 ESA) is the exchange current density measured at a reference temperature and reference dissolved oxygen concentration (solubility) . is the reactant concentration on the catalyst surface calculated from (20) (or (21) for hydrogen). (mPt2 g−1) is the catalyst-specific accessible electrochemical surface area (covered by electrolyte) in the electrode measured using cyclic voltammetry (31.55 and 35.45 m2 g−1 for cathode and anode electrodes utilizing 50% Pt/C and 20% Pt/C, resp.), which corresponded to 50%−30% of the given value by the manufacture (ETEK, USA) of 128 m2 g−1 for 20% Pt/C and 86 m2 g−1 for 50% Pt/C when placed in the electrode structure [49]. can also be estimated by multiplying the ionomer’s volume fraction in the catalyst layer by the catalyst ESA.
is the catalyst loading, which corresponds to the weight of platinum per unit geometric area (mg cm−2). The product is the roughness factor (dimensionless), which is the electrochemical surface area divided by the electrode geometric area, referred to earlier as . The units of the product or will be A cm−2 (current density per electrode geometric area).
is the pressure coefficient or the reaction order with respect to oxygen in phosphoric acid/PBI, a value of 1 (first order) was reported in [24, 32].
, the activation energy of oxygen reduction in hot phosphoric acid, was found to be independent of phosphoric acid concentration [47 equals to 72.4 kJ mole−1 [32]]. Values for vary in the literature (Table 2), depending on the temperature, phosphoric acid concentration, and the dissolved oxygen concentration. Most values were in the range of 10−8 A cmPt−2, and value from [31, 47] was used.
Similarly, a value of = 0.144 A cm−2 at 160°C ( mole L−1), kJ mole−1, , and for hydrogen oxidation in phosphoric acid was obtained from [39].
3.5. Reformate Fuel
One of the major advantages of high temperature operation in PEMFCs is the tolerance and therefore the systems ability to operate with reformate gas feed without considerable loss in performance. Reformate gas composition varies depending on the hydrocarbon used in the reformation process. A typical diesel reformate composition contains 3% CH4, 19% CO2, 2.2% , and 20% , and the rest is hydrogen.
The model was developed to account for the poisoning effect of and possibly methane on the catalyst surface. poisoning effects on in phosphoric acids system have been studied by several groups [40, 41, 50, 51]. There are two main effects of on anode performance; the first is a dilution effect already accounted for in the Butler-Volmer equation, and the other is a kinetic effect arising from a reduction in active surface area because of adsorption on the surface.
The rate determining step for hydrogen oxidation is the dissociation reaction (Volmer) giving rise to the second-order dependency of hydrogen oxidation rate on surface coverage and Vogel et al. [39] suggested the following equation: is the exchange current density for hydrogen oxidation after poisoning, is the exchange current density for hydrogen oxidation without presence, and is the surface coverage by given byis reported to vary linearly with ln (/H2) [39, 40, 50], although various values for the constants of their relation dependency have been reported. Kohlmayr and Stonehart [51] suggested that coverage is independent of temperature, in the range from 100 to 150°C in phosphoric acid media. Considering the heat of adsorption, it is expected that coverage should exhibit an exponential dependency on temperature. Dhar et al. [50] suggested the following equation for the variation of coverage variation with temperature in H3PO4 In this work, (35) was used to express coverage on a platinum surface and its variation with content, temperature, and hydrogen content.
In the case of carbon dioxide only its dilution effect was considered in the model, ignoring any effect of reduction of carbon dioxide to carbon monoxide by hydrogen at the studied temperature [52].
Under the circumstances of insufficient data on the effect of methane on the platinum surface in phosphoric acid, two cases were considered.
Case 1. Sustersic et al. [53] showed that methane adsorption on platinum surface started from a potential 0.2 V versus NHE and reached a maximum adsorption peak at 0.25 V (versus NHE). In HT-PEMFC the anode potential does not exceed 50 mV (even at high current densities of 1.5 A cm−2) which is far below the minimum potential suggested for methane adsorption. Niedrach [54] found that methane had low adsorption on platinum (0.28 at STP) compared to that of other saturated hydrocarbons (moderate) except ethane, while unsaturated hydrocarbons exhibited very high adsorption: 0.79 (at STP) for propylene and 0.91 (at STP) for cyclopropane. This suggests that consideration should be taken for nonsaturated hydrocarbons in the reformate gas, and their concentration should be determined.
Case 2. Taylor and Brummer [55] studied methane adsorption in 12 M H3PO4 at 130°C. Similarly, they found a maximum adsorption peak at 0.25 V versus NHE; however, on the contrary, they observed methane adsorption at low potentials of 0.1 V and lower (versus NHE) some 50% of that at 0.25 V. This suggests that methane adsorption might occur even at the low operating anode potentials.
Hsieh and Chen [56] studied methane oxidation on Pt in 1 M . The rest potential was 100 mV versus NHE, the exchange current density was 10−7 A cm−2 at 80°C and 1 atm methane, and the activation energy was 125 kJ mole−1. The activation energy was 3.5 larger than that of hydrogen oxidation, and the exchange current density was six orders of magnitude lower than that of hydrogen oxidation. The rate of methane oxidation compared to hydrogen is negligible, so it can be treated as an inert species with only dilution/concentration effect (Case 1) or considering the worse case scenario as adsorption similar to that for (Case 2).
3.6. Gas Crossover Polarisation
In reality the observed cell OCP is lower than the estimated value from the thermodynamics due to the effect of crossover and other phenomena (carbon corrosion, platinum oxidation, Pt-OH, etc.). Considering crossover current density (at ) leading to an overvoltage , given by where α is the transfer coefficient and is the exchange current density.
The above equation can be expressed in the following form:
The observed open circuit will be equal to + ; typically for a PBI membrane with thickness 80–40 μm the observed OCP with oxygen at 150°C was in the range of 0.88–1.0 V. OCP values of 0.95 V [33] and 0.9 V [15, 17] have been previously used for PBI PEMFC models.
It should be stressed here that the observed OCP cannot replace as is solely measured at . The effects of crossover on polarisation curves can be ignored, as it is only important at very small currents and becomes negligible when polarising the electrode, where the total current becomes equal to () as can be seen in (38) in other words the overvoltage losses due to hydrogen crossover become negligible
3.7. Conductivity and Losses
To express the proton conductivity of acid doped PBI at different temperatures a modified Arrhenius equation has been suggested for polymer/acid (salt) complexes [57–61] Expressions for , the pre-exponential factor, and , the activation energy, which are humidity and doping level dependant are given in the appendix and the corresponding constants to calculate and are given in [58, Table 4].
Polynomial functions were fitted for and based on the relative humidity at a given doping level (membrane doping level of 5.6 PRU in this case), given in the Appendix.
The conductivity obtained from (39) along with the membrane thickness of (40 μm) were used to determine the losses through the membrane. On the other hand the losses (protonic) through the catalyst layer were obtained from the catalyst layer thickness (15 μm) and the effective phosphoric acid conductivity (product of phosphoric acid conductivity at the studied condition and its volume fraction in the electrode). The conductivity of phosphoric acid varies with temperature and relative humidity (RH) (or phosphoric acid concentration (wt%)).
The following functions for conductivity were given [45] for %:
For %:
Matlab V.7.3 and Simulink V.6.5 (The MathWorks, U.K) equipped with Ordinary Differential Equation solver (ODE 45) were used to solve the model governing equations. The model results were compared to experimental data for a 50% Pt/C (0.4 mgPt cm−2) cathode and 20% Pt/C (0.2 mgPt cm−2) anode. The experimental data were collected using a cathodic sweep from OCP to 0 V under pseudosteady-state conditions. High resolution experimental data were sampled every 1.0 mV step in order to allow for good comparison with the model results.
4. Results and Discussion
4.1. Influence of the Thin Film
We carried out simulations of the PBI-based fuel cell using a model in which thin film characteristics were not included. The model was based on a two-dimensional architecture [35] in which flow channel effects can be accounted for but with high gas stoichiometries reduces to a 1D model. The model in general predicts the effects of parameters such as temperature and gas pressure as shown in Figure 2, that is, performance improved with higher temperature and with higher oxygen partial pressure. However, as shown in Figure 3 under conditions where experimentally a limiting current density was observed, the model failed to predict this behaviour and overpredicted performance. Consequently for modelling PBI PEMFC the influence of mass transport of reactants through an electrolyte film is required as shown in the following sections.
| (a) Temperature effect. (a) ♦ : 170°C; (b) : 150°C; (c) ● : 120°C |
| (b) Pressure effect. (a) ♦ : 2 bar; (b) atmospheric pressure |
4.2. Oxygen Partial Pressure Effects
Figure 4 compares simulated and experimental data for PBI-based HT-PEMFC at different oxygen partial pressures: from air (atm) to air (1 bar) and pure oxygen. Large increases in voltage were seen in both modelled and experimental data due to enhancement in kinetics and mass transport when increasing oxygen concentration (partial pressure). This effect was due to the low oxygen permeability in hot concentrated phosphoric acid. The experimental observed limiting current increased from (~1.4 to 2.1 A cm−2) when doubling oxygen partial pressure (1 atm to 1.99 atm). The model data and the experimental data were in good agreement. However, small differences were observed in the slope of the polarisation curves, at high current densities, where the model overestimated the cell performance with air and air (1 bar pressure) whilst not under oxygen operation.
The slope of the polarisation curves is affected by a combination of kinetic, mass transport, and effects. The observed difference of the slope at high current densities in this case was probably caused by unknown kinetic or mass transport effects, and not effects, as a good match in the slope was observed at high current densities under oxygen operation. In order to improve the fit between the experimental data and the model prediction, it would seem that a correction at high current densities should be made to the kinetic effect and mass transport (higher ).
To introduce increased mass transport losses, the reaction order () with respect to oxygen in phosphoric acid was increased from the commonly reported value of 1.0 [32] to 1.375; values above 1.0 have also been reported for PBI/H3PO4 interface [25]. The reported values of reaction order were measured at constant overvoltage (and not constant voltage ), that is, they were measured from the exchange current density ratios (), where was measured at equilibrium potential which in turn is affected by a concentration change (Nernst equation) as will be discussed shortly.
Figure 5 shows the effect of the higher reaction order (1.375) on the cell voltage polarisation curve and maintaining the standard value of alpha (0.75). An improvement in the simulation and experimental data agreement at various oxygen pressures is seen, where the model predictions followed the experimental data very closely even at high current densities and different oxygen concentration
A higher value of gamma than 1.0 can be caused by factors such as the following.(i)Thermodynamic Effect of Oxygen Concentration Losses:
Whilst the oxygen concentration influence was accounted for in the Butler-Volmer (32), the cell reversible potential was considered constant over the cell polarisation where the equilibrium concentrations at rest potential (zero current) was considered. An additional change in the reversible potential will occur due to lower oxygen concentrations (see (3)). This shift is equal to () ln () as given by the Nernst equation.
Thus, rather than writing the BV equation for high current densities as:
The overpotential is given by:
This can be rewritten as: Where is equal to or in this case , the value used for gamma in the latest simulation.
(ii) Losses of Surface Area (ESA) due to Oxygen Starvation:
Losses in the accessible electrochemical surface area might occur due to oxygen starvation, where at very low oxygen concentrations it will not be able to reach parts of the catalyst surface even though it is available for reaction (covered with electrolyte). This effect would lead to a further decline in performance at high current densities and low oxygen concentrations (in comparison to gamma = 1). Unfortunately, such an effect is not accounted for in this model and is one of the limitations of the macrohomogeneous model, where it is assumed that the average accessible catalyst layer is constant.
However if such effects are present then, the ESA can be written as function of ln () which would effectively cause an increase in apparent reaction order.
4.3. Porous Media Effects
The effect of mass transport losses due to diffusion through porous media can be seen by comparing cell performances using air (O2–N2) and heleox mixture (O2–He) as oxidant at the cathode. While both mixtures contain the same oxygen concentration (), the inert gas used in the mixture is different. The binary diffusion coefficients of oxygen-helium and water-helium are four and eight times higher than those of oxygen-nitrogen and water-nitrogen, respectively.
The enhancement in cell performance observed when switching from air to heleox is caused by mass transport enhancement through the porous structure. Generally a larger difference between the limiting currents using air or heleox is due to a greater contribution of the porous media mass transport losses or the smaller the value of (). Similarly, to the earlier observation (oxygen partial pressure effect), using and , the model provided good predictions of cell performance; however it overestimated the performance at high current densities.
Figure 6 shows that when increasing the gamma value to 1.375, whilst maintaining the standard value of alpha (0.75), an improvement in the model fit with the experimental data is achieved. The model predictions followed the experimental data very closely even at high current densities and with different oxidant: air or heleox.
The model also successfully estimated the limiting current values with air and heleox. This suggests that the chosen value for porosity of 22% ( [18]) is appropriate for the studied electrode.
4.4. Thin Film and Electrolyte Distribution
As discussed earlier incorporation of the thin film in the model is essential to enhance the model ability to predict cell performance at high current densities and low oxygen concentrations. The electrolyte content is a crucial parameter for optimising the three-phase boundaries and therefore the electrode performance. Increasing the electrolyte content in the electrode will increase the accessible ESA and increase the film thickness whilst reducing the porosity. The electrolyte film thickness can be estimated using (19). Values of 50 and 160 nm were obtained for the studied cathode and anode, respectively.
Figure 7 shows the effect of cathode film thickness on cell performance under air operation. Increasing the film thickness by factors of 2 and 4 led to a severe impact on cell performance even at relatively small current densities and as a consequence the limiting current fell from 1.32 to 0.78 and 0.43 A cm−2, respectively. This shows the crucial impact of electrolyte content on cathode and cell overall performance. Similarly, decreasing the film thickness enhanced the overall performance due to enhancement in mass transport and correspondingly an increase in the observed limiting current. This behaviour agreed with reported experimental observations [62], where electrodes with the highest observed limiting current showed the best overall performance over the entire potential range.
While the above analysis assumed fixed ESA, conductivity, and porosity, in practice there is a minimum film thickness below which the ESA will fall rapidly (conductivity also will fall, while porosity will increase).
Considering a case where the film thickness is maintained constant and the ESA is halved, this corresponds to decrease in the acid content (volume fraction) in the cathode layer by half. Such a change led to a dramatic fall in cell performance in both the kinetic (activation) and the mass transport regions of the polarisation curve where the limiting current decreased from ~1.4 to 0.78 A cm−2. On the other hand if the electrolyte content is increased (for example doubled) above the optimum content, where no further significant increase in ESA is achieved, then the overall performance will decrease due to mass transport effects caused by a doubling in the film thickness.
This behaviour agrees with the experimental finding of the importance of acid loading (content) in the cathode layer [62]: excess acid (high doping case) led to severe mass transport limitations, and heat treatment (low doping case) led to dramatic impact on both the kinetic, and mass transport losses.
The model showed (data not reported) that a variation in anode film thickness had a less crucial effect on performance than that of the cathode film, even though the anode initial film thickness was ~3 times higher than that of the cathode: no limiting current was observed with oxygen operation (no mass transport limitation from the cathode) up to 2 A cm−2 (no experimental data is available beyond 2 A cm−2 due to instrument limitations). The behaviour was caused by the high hydrogen concentration in the anode feed (pure hydrogen) in comparison to air operation in the cathode. The model predictions show that for the anode catalyst layer, the acid content was far less crucial than that of the cathode; the anode catalyst layer was able to accommodate more acid before considerable mass transport limitations became visible. Also a low PBI content in the anode catalyst layer (5% wt) did not result in mass transport limitation up to the studied current density of 2 A cm−2. However, a larger PBI content did lead to severe mass transport limitations and an observed limiting current (due to anode) at low current densities.
4.5. Comparison with 1/2D Model
Figure 8 shows the effect of oxidant on the predictions of a 2D model which includes a thin film model [35]. This model as would be expected gives accurate simulations of the experimental data. The differences between predictions of this model and the simpler model described in this paper (Figures 5 and 6) are very small and fall in the range of error of estimation of the model parameters (from experimental data) and the experimental cell voltage/current density measurements. Thus under conditions when excess oxidant and air are used, it would seem acceptable to adopt a simpler model as described in this paper. When however incorporating the effect of gas supplied near stoichiometric requirements a more rigorous 2D (or 3D) model should be used.
(a)
(b)
4.6. Diesel Reformate
The use of diesel reformate instead of pure hydrogen is expected to have its major influence on the behaviour of the anode polarization, and thus this factor is considered in the simulations.
A comparison between the model and experimental and CO2 effects on the anode performance at 150°C is shown in Figure 9. The voltage losses at a current density of 1.5 A cm−2 were for both experimental and modelled data 8, 12, and 22 mV when switching from pure hydrogen to 20%, 33%vol CO2, and 2.5%vol , respectively. Similarly, Li et al. showed that small voltage losses of less than 10 mV with current densities up to 1.3 A·cm−2 were observed for PBI/H3PO4 with levels of 1 and 3 % at temperatures of 150 and 200°C, respectively [52]. This results in losses of power density of 12, 18, and 33 mW cm−2, respectively. Very good agreement was observed (Figure 9) between the experimental and the model data over the entire studied compositions. A small deviation was observed at high current densities where the experimental data approached a plateau (curve) whilst the model data continued in the expected, close to linear variation. This behaviour was also observed with pure hydrogen. Due to fast hydrogen oxidation kinetics and high hydrogen concentrations used, losses were the main contribution for anode losses at high current densities. Whilst no water flux through the membrane was considered in this model, the flux of gaseous water product from the cathode to the anode will occur at high current densities, driven by the large gradient in water content which in turn will cause the observed enhanced proton conductivity of the anode (anode feed—in this experimental case—had a very low humidity at 150°C with a humidifier temperature of ~16°C). Water permeability of to 200 cm3 (STP) m m−2 s−1 Pa−1 for phosphoric acid doped PBI membranes was reported [37] at elevated temperatures of 125–150°C.
Figure 10 shows the simulated anode performance at 150°C with various reformate compositions considering Case 1 (concentration effect only for methane).
The following observations can be drawn.(i)Methane has a worse dilution effect than CO2, at the same concentration of 22%. An increase in anode potential in comparison to pure hydrogen of 11 mV was observed at 1.5 A cm−2 when using methane (a power density loss of 16.5 mW cm−2) compared to that of 8 mV when using CO2 (loss of 12 mW cm−2).(ii)The addition of the separate dilution effect of inert species (CO2, , etc.) with the poisoning and dilution effects of poisonous species () does not add up to the simultaneous effects of a mixture of both species. This is explained by the coverage as proportional to ln ([]/[H2]) [40, 50], which means that any fall in hydrogen concentration, due to dilution by inert species, will lead to higher coverage at the same concentration and therefore enhance the poisoning effect.(iii)A 30 mV potential difference between pure hydrogen and the diesel reformate composition at a current density of 1.5 A cm−2 was obtained considering only the dilution effect for methane (Case 1), reflecting a fall in power density of 45 mW cm−2.(iv)The potential difference between pure hydrogen and diesel reformate increases to 50 mV at the same current density, when adding 20% water to the gas stream (no information was available for water concentration in the reformate mixture). However, 20% water content should be reasonable for steam reforming case, which means a loss of 75 mW cm−2 in power density at 1.5 A cm−2. Water dilution is significant as the hydrogen concentration falls from 75% to 55% leading to intensified poisoning effects, whilst initially increasing water concentration in the stream is thought to be beneficial in terms of proton conductivity enhancement. The maximum allowed water content in the gas steam should also be considered to avoid acid losses (washout) and corresponding losses of system overall conductivity.
Figure 11 shows the modelling results for anode performance at 150°C under various reformate compositions, considering Cases 1 and 2 simultaneously. The voltage loss compared to pure hydrogen increased from 30 to 35 mV when switching from Case 1 to Case 2 (methane treated similar to as poisoning species) with 5% methane in the gas feed at 1.5 A cm−2. The addition of 20% water to the stream led to a voltage loss (compared to pure hydrogen) of 50 mV and 55 mV for Cases 1 and 2, respectively. This means a power density loss of 75 to 82.5 mW cm−2. Considering a system power density peak of 250–300 mW cm−2 when operating with air at atmospheric pressure, this would mean a loss of ca. 30% of the generated power density.
5. Conclusions
A model of a high temperature PEM fuel cell using PBI membranes has been developed using thermodynamics, transport, and kinetic equations. The model considers mass transport through a thin film electrolyte as well as through the porous media. The incorporation of the thin film layer was crucial for model accuracy, particularly at high current densities and low oxygen concentrations (air).
Oxygen permeability through the thin electrolyte film varies for a given temperature depending on the equilibrium vapour pressure of the product water above the thin film and, correspondingly, the operating current density, due to variations in oxygen permeability with phosphoric acid concentration (electrolyte film).
The model showed very good agreement with experimental data under various operating conditions and oxygen concentrations. The model emphasises the importance and sensitivity of the electrolyte content on electrode performance, particularly on the cathode. The model is a suitable tool for optimisation of the electrode performance and helps in understanding the reasons behind performance limitations. For example, the model showed that acid doped PBI is not suitable as electrolyte additive for the cathode catalyst layer and that pristine phosphoric acid is preferred. Similarly, the anode can tolerate small quantities (loading) of PBI although excess PBI will lead to severe mass transport limitations.
Appendix
A summary of and constants for (28) is in given Table 3:
The pre-exponential factors and in (23) and (24) are calculated usingwhere () denotes (). Expressions for the calculation of constants and in (39):
List of Symbols
Latin Letters
| or : | Oxygen concentration in the thin film at platinum-electrolyte boundary/mole cm3 |
| : | Oxygen concentration in the thin film at gas pore-electrolyte boundary/mole cm3 |
| : | The reference standard H2 concentration (expressed as 1/RT)/mol cm−3 atm−1 |
| : | The reference standard O2 concentration (expressed as 1/RT)/mol cm−3 atm−1 |
| : | Oxygen diffusion through ionomer/cm2 s−1 |
| : | The effective binary diffusion/cm2 s−1 |
| or : | Activation energy/Kj mole−1 |
| : | Standard cell reversal potential/V |
| ESA: | Electrochemical surface area/m2 g−1 |
| F: | Faraday constant/coulomb mole−1 |
| : | Standard Gibbs free energy of a reaction/kJ |
| GDL: | Gas diffusion electrode |
| : | Henry’s constant for oxygen solubility in the electrolyte m3 atm mole−1 |
| : | The exchange current density/A cmPt−2 (ESA) |
| : | Exchange current density for hydrogen oxidation after poisoning |
| : | The exchange current density at reference temperature/A cmPt−2 |
| : | Crossover current (rate)/A cm−2 (geometric) |
| : | Limiting current density/A cm−2 (geometric) |
| : | Current density/A cm2 (geometric) |
| : | The catalyst loading, weight of platinum per unit area/g cm−2 |
| : | The molecular weight of the gas/g mole−1 |
| : | The mass of acid per unit area/mg cm2 (geometric) |
| : | The mass of PBI per unit area (PBI loading)/mg cm2 (geometric) |
| : | The molar flux of species /mole s−1 cm−2 (geometric) |
| : | Number of electrons involved in the reaction |
| : | Oxygen partial pressure/atm |
| : | The gas critical temperature and pressure/atm |
| : | The saturation vapour pressure of pure water at given temperature/atm |
| PRU: | Number of doped acid molecules per repeat polymer unit |
| : | Gas constant/J K−1 mol−1 |
| RH: | Relative humidity/% |
| : | The specific surface area of carbon per unit area/dimensionless |
| : | The surface area/m2 g−1 |
| : | Temperature/K |
| : | Gas critical temperature/K |
| : | The reference temperature where is measured/K |
| : | Cell potential/V |
| : | Potential versus SHE/V |
| : | Phosphoric acid weight concentration/wt % |
| : | The molar fraction of species |
| : | Phosphoric acid molar concentration/mole %. |
Greek Symbols
| α: | Transfer coefficient |
| and : | Transfer coefficient for anodic and cathodic reactions |
| : | The water activity |
| : | The reference standard water activity (unity) |
| : | The pressure coefficient or the reaction order |
| : | Average ionomer film thickness covering a catalyst particle/m |
| : | The porosity/% |
| : | Overpotential/V |
| : | Anode overpotential/V |
| : | Cathode overpotential/V |
| : | Cross over potential loss/V |
| : | The density/g cm3 |
| : | Conductivity/S cm−1 |
| : | The tortuosity |
| : | Surface coverage by . |
Acknowledgments
EPSRC through the SUPERGEN Fuel cell Consortium and the EU through the FP6 project FURIM supported this paper.