This work presents a pseudo-two-dimensional proton exchange membrane fuel cell (PEMFC) model incorporating Nafion ionomer structure–property relationships in the cathode catalyst layer (CL) to capture and explain losses at low Pt loading. Structural data from neutron reflectometry and thin film Nafion conductivity measurements predict variations in the oxygen diffusion coefficient and ionic conductivity with changing CL ionomer thickness and Pt loading. By including these structure–property relationships, predicted polarization curves agree closely with previously published experimental data from cells with Pt loadings between 0.025 and 0.2 mg/cm2. Results demonstrate that structure–property relationships based on physically measurable ionomer and CL properties provide a feasible interpretation of PEMFC CL phenomena for a range of Pt loadings and help explain previously unaccounted-for losses at low Pt. Results also show that simulations must account for surface species coverage variations in order to properly capture the kinetic losses. Finally, results suggest that an increase in ionomer thickness surrounding the C/Pt surfaces may lead to improved cell performance due to improved ionic conductivity.
Proton exchange membrane fuel cells (PEMFCs) provide clean and efficient conversion of chemical energy into electrical energy. Fuel cell electric vehicles (FCEVs) are proposed as a competitor to battery electric vehicles (B-EVs) for decarbonized transportation. FCEVs have an advantage in categories such as range, weight, and refueling times when compared to B-EVs . Although FCEVs are currently available, challenges in infrastructure, performance, durability, and cost hinder wide-spread adoption [2–6].
The PEMFC performance, cost, and durability are all closely tied to the Pt catalyst used in the fuel cell electrodes. Efforts to lower the amount of Pt loading [7–9] or to replace the catalyst altogether [10,11] have shown unacceptable cell performance or low durability. The U.S. Department of Energy (DOE) has set a goal to lower the PEMFC Pt loading to 0.125 mg/cm2 or 0.15 mg/kW. However, current state-of-the-art PEMFCs with loadings of roughly 0.13 mg/cm2  suffer from performance losses. Electrodes with lower Pt loadings exhibit large transport losses that are not yet sufficiently explained in the literature . Some researchers propose resistances at the ionomer/gas interface and ionomer/Pt interface as a source of these resistances [14,15], although other work suggests limited impact for the former . Although incorporating these resistances provides reasonable agreement to data for existing models, their cause is not yet fully understood.
Multiple PEMFC limiting phenomena occur in the cathode catalyst layer (CL), where nano-thin sheets of ionomer (<10 nm thick) coat carbon-supported Pt nanoparticles. Voltage losses in the CL are attributed to water transport, flooding, electro-osmotic drag, resistances associated with the Pt surface, and transport through the nano-thin ionomer [17–19]. Because of the difficulty in experimentally measuring transport properties in films at this scale, modelers commonly use parameters taken from bulk ionomer materials. However, recent literature demonstrates significant differences between nano-thin and bulk-scale ionomers such as Nafion. Limited direct measurements of the ionic conductivity have been reported [20,21], and the thin film oxygen diffusion coefficients are largely unreported due to experimental limitations. That said, recent studies provide significant insight into the structure of thin film Nafion. Results show complex, layered interfacial structures, reduced water uptake, and reduced plasticity for films thinner than 60 nm due to confinement and substrate interactions [22,23]. However, there is limited experimental evidence directly connecting these structural properties and associated PEMFC CL fabrication/processing conditions to the thin film ionomer transport properties.
Detailed physical models provide insight into these losses and can help guide future PEMFC designs. Limiting processes in the CL occur over a wide range of length scales (ranging from nanometers to millimeter). Therefore, models typically implement effective cathode CL microstructures to capture relevant nm-scale phenomena within a continuum-scale model framework, using one of the three approaches: discrete-catalyst, thin-interface, and heterogeneous [14,24–26]. These types of models are mature in the PEMFC community, and a wide range of phenomena have been incorporated into them since the 1970s. The history of physically based PEMFC model development is well covered by a pair of reviews by Weber et al. [27,28]. Recent models have focused on phenomena such as the impact of microstructure [25,29,30], water management [30–33], Pt distribution [34,35], and species transport [32,36]. Broadly speaking, heterogeneous models offer the most detail and include effective microstructures such as the “core–shell” and “flooded agglomerate” models. As illustrated in Fig. 1, the core–shell microstructure takes as the computational domain a Pt-covered carbon core wrapped in an ionomer shell. The flooded agglomerate model considers a collection of Pt-covered carbon particles packed inside an ionomer shell that has been flooded with ionomer. Each approach has strengths and limitations. The core–shell model assumes that reactants and products only need to travel through a thin film of ionomer before reaching the gas phase. In reality, some of the Pt surface is“buried” within agglomerates, making reactant access more challenging. The flooded agglomerate model, on the other hand, underestimates the degree of Pt which is readily accessible to gas-phase reactants. The flooded agglomerate model requires artificially large agglomerates (hundreds of nanometers, encompassing tens of carbon support particles) in order to be compatible with continuum-scale finite differencing.
Actual microstructures, as observed in transmission electron microscope (TEM) images , lie somewhere between the two approaches. Agglomerates of roughly 3–5 carbon particles, flooded with ionomer, are typically observed. These agglomerates provide easy access to most of the Pt catalyst, but are typically highly non-spherical, with characteristic lengths ranging from roughly 50 to 300 nm. Such shapes are not readily represented in continuum-scale simulations, at present. Detailed computational fluid dynamics simulations can accurately represent high-fidelity microstructures, but scaling these up to the PEMFC stack level is a challenge. The effective geometries in Fig. 1, while not fully descriptive of the true CL microstructure, therefore provide reasonable computational performance at device-relevant scales.
In this work, we present pseudo-two-dimensional (2D) PEMFC simulations to explore limiting cathode CL phenomena at low Pt loading. Previously published structural data and ionic conductivity measurements are correlated with one another to develop structure–property relationships for the CL ionomer. These relationships capture transport losses related to variations in Pt loading, ionomer thickness, and CL microstructure, for both the core–shell and flooded agglomerate model frameworks. Results show that the core–shell geometry provides good agreement with experimental results for Pt loadings in the range of interest by the DOE. The flooded agglomerate models, on the other hand, require unrealistic ionomer transport properties and do not provide good agreement with data at low Pt loadings. Moreover, results demonstrate that microstructure models informed by TEM images and parameters based on physically meaningful structure–property relationships can capture complex transport and electrochemical kinetic phenomena in the PEMFC CL. With additional refinement, such modeling tools therefore provide a viable pathway to identify future PEMFC designs and operating conditions for low-cost, high-performance PEMFCs for EVs and other consumer applications.
2 Model Formulation
The model solves physically based conservation equations using a finite-volume approach to discretize the cathode GDL and CL. Within the GDL, one-dimensional porous flux expressions are used to compute evolution of the gas-phase species. For the CL, a pseudo-2D Newman-type model is used to track the PEMFC state variables as a function of the electrode depth (y) and the spherical solid particle radius (r). The mass density of oxygen absorbed in the Nafion (, ) is tracked in the y and r directions. The surface site coverage of the Pt catalyst (), when considered, is tracked in y and r (note that is not tracked in r for the core–shell model, since there is no radial distribution of Pt in the core–shell model; see below for details). The Nafion electric potential (ΦNaf, V) and the gas-phase species mass densities (, ) are tracked only in the y dimension.
The model domain is illustrated in Fig. 1. Conservation equations are derived as time-based derivatives, solved numerically on a discretized series of finite volumes representing the model domain. The model equations capture the following physical processes: (i) gas-phase convective and diffusive transport, (ii) transport of ions and dissolved neutral species in the ionomer electrolyte, and (iii) reversible chemical and electrochemical reactions at ionomer-gas and ionomer-Pt interfaces. Similar model implementations are common in the literature , but what sets this model apart is the derivation and implementation of the ionomer structure–property relationships, as presented in Sec. 3.1. Evaluating the unsteady governing equations for these processes at a constant, user-defined current density and integrating over a long time span (100 s) until steady state provides the model output. Although the proton exchange membrane (PEM) is shown in Fig. 1, it is modeled here as a simple resistor with constant resistance. The anode is additionally not shown or explicitly modeled here, but its reversible electric potential is calculated and used as a reference so that the half-cell model can be directly compared to full-cell PEMFC experiments.
2.1 Model Assumptions.
The following assumptions are employed in the model:
Water is in equilibrium between the Nafion and gas phases, and the gas remains at a fixed, steady relative humidity (RH).
Solid core–shell and/or flooded agglomerate structures are spherical in shape.
Uniform porosity, microstructure, and Pt distribution throughout the CL.
Electric potential is radially uniform within the spherical solid particles.
The carbon and Pt phases are isopotential through the CL depth.
Aside from neglecting fluctuations in water content via the second assumption, these assumptions are common throughout PEMFC modeling [32,38,39]. Although water transport and flooding are critical phenomena which limit cell performance with low Pt loading [31,40,41], they are neglected here so that the impact of structure–property relationships can be evaluated in the absence of confounding impacts from uncertain water dynamics. Incorporating more complex water dynamics and electrolyte phenomena will be the topic of future modeling work from our group.
2.2 Catalyst Layer Microstructure Modeling.
Two different structures are considered as representative domains within the CL: core–shell and flooded agglomerate. Each domain has been commonly used throughout PEMFC modeling literature [14,26]. Since the actual microstructure of a PEMFC CL is heterogeneous and non-uniform over very small length scales, these structures attempt to capture the average behavior within each discretized volume. The zoomed in section of the CL in Fig. 1 illustrates each domain. Note that the domains are not drawn to scale, but scale bars in the figure show the approximate relative sizes. A variety of values have been used with each of these geometries in previous literature, even beyond the scales provided here.
2.3 Conservation Equations
2.3.1 Gas-Phase Mass and Species.
2.3.2 Electrolyte Mass and Species.
2.3.3 Conservation of Charge.
2.4 Kinetics and Reactions.
2.5 Surface Coverages.
2.8 Cell Voltage.
2.9 Boundary Conditions.
3 Parameter Estimation
Based on TEM images and literature, a minimal number of parameters, all physically based, are used in the model. Briefly, the CL microstructure is defined by user inputs for the carbon particle radius, Pt loading, Pt radius, porosity, CL thickness, and Nafion shell thickness. In addition, the flooded agglomerate model requires additional inputs for the agglomerate radius and the volume fraction of carbon within the agglomerate. All other microstructure parameters are derived from these inputs. GDL and CL gas permeabilities Kg are taken from Ref.  and scaled linearly with the porosity. Input parameters that describe the core–shell and flooded agglomerate microstructures are summarized in Table 1. The parameters used in the model that are derived from this base set are included in Tables 2 and 3 for the core–shell and flooded agglomerate domains, respectively.
3.1 Structure–Property Relationships for Thin Film Catalyst Layer Ionomers.
Experiments show that transport properties for both protons and oxygen within the ionomer change as a function of the temperature, RH, thickness, and substrate interactions [21,45]. Specifically, variations in thickness and substrate type result in confinement effects and nano-structural rearrangements that impact water absorption and overall material structure [46,47]. Moreover, structural data from neutron reflectometry (NR) observe multi-layered Nafion structures at hydrophilic substrate interfaces [22,23,48]. These layers alternate between water-rich and water-poor compositions, with an oscillation that damps toward a uniform, thicker outer layer away from the substrate, beyond ca. 10 nm.
This model is unique in its use of structure–property relationships within the thin film CL ionomer, i.e., Nafion. These structure–property relationships are derived by combining three separate datasets:
The three structure–property methods used here make varying use of these scaling factors, as described as follows.
3.1.1 Method 1: Layered.
The “layered” model hypothesizes that lamellae persist at every Nafion/solid interface. Effective diffusion coefficients for O2 absorbed in Nafion are calculated for each layer, as described in Eq. (24). Because the film is much thinner than the carbon particle radius, layer diffusion coefficients are summed using a linear series resistance to calculate an effective diffusion coefficient for the film. The ionic conductivities are taken directly from measurements on similar thickness films on silicon, as measured by McCreery et al. , interpolating when necessary.
3.1.2 Method 2: Uniform.
Although layered interfacial structures have been observed on Si and Pt substrates using NR [22,23,51,52], similar definitive studies do not currently exist for carbon substrates. While it is possible that lamellae do not exist in PEMFC CLs, it is still known that confinement effects impact water absorption and transport through the Nafion electrolyte . To account for this possible no-lamellae case, a “uniform film” method is implemented, which calculates uniform, isotropic transport properties for the thin Nafion films, using Eqs. (23) and (24), where λ is taken from NR measurements of similarly thick layers .
3.1.3 Method 3: Mixed.
In all likelihood, if lamellae exist in the PEMFC CL, they are present to a limited extent. Layered structures have been observed at Nafion/Pt interfaces [51,52], consistent with the more hydrophilic Pt surface. The “mixed” model, then, proposes that lamellae are present at the Pt interface but that regions of exposed carbon are covered by “uniform” films. Because O2 only diffuses in regions near the Pt, for this model is calculated as in the “layered” approach. σio, meanwhile, is calculated as an average of the “uniform” and “layered” conductivities, with the two values weighted in parallel by the surface area fractions of carbon and Pt, respectively. The conductivity is therefore a function of the Pt loading. Higher loading correlates with higher conductivity, due to the higher water content of the interfacial layers.
3.2 Transport of Dissolved Nafion Species.
Figure 2 also demonstrates one additional degree-of-freedom for dissolved oxygen transport in the Nafion. For the core–shell model, a “capture angle” is considered to represent lateral diffusion of the O2 within the Nafion film. Because the Pt does not cover the entire carbon surface, only a fraction of the Nafion coating will be active for O2 transport. If transport parallel to the film surface is negligible, the fraction of the Nafion area available for O2 transport at any given radius will be equal to the fraction of the carbon surface area covered by Pt. However, as “lateral” O2 transport becomes more prevalent, a greater fraction of the Nafion area becomes available. In reality, the transport field lines will be curved. As a first approximation, an effective capture angle is used here to represent the phenomenon linearly, as shown in Fig. 2. Low values represent the mostly radial diffusion, whereas larger results in greater Nafion area available for O2 transport. Because O2 transport in the flooded agglomerate model is isotropic, with respect to the Pt surface orientation, this parameter is only implemented in the core–shell model.
4 Simulation Procedure
4.1 Numerical Integration.
Written as a system of ordinary differential equations in python,2 the model is solved numerically using the solve_ivp integrator available in the SciPy package . A sufficiently large simulation time is used to ensure that steady-state conditions are reached. Running the model over a range of external currents between 0 and 2 A/cm2, polarization curves are compared to experimental results with Pt loadings between 0.025 and 0.2 mg/cm2 from Ref. .
The open-source software library Cantera is used to manage the model’s thermodynamic, kinetic, and gas-phase transport parameter calculations . Specifically, Cantera is used here to calculate mixture-averaged diffusion coefficients and viscosity for ideal gases, net production rates from the previously mentioned kinetic expressions, gas phase pressure and molecular weights. After setting the state of the Cantera phase objects using state variables, these properties can be calculated via simple and generalizable calls to Cantera library functions, simplifying the code structure significantly and increasing the repeatability of the model code. Cantera also makes it simple to change between the one-step and two-step charge transfer mechanisms, with no changes to the model code itself .
5 Accessing the Model
To advance open science practices related to accessibility, transparency, and repeatability, the model code is made freely available via GitHub . We encourage any interested parties to download, use, adapt, or extend this model. Contributions to the model via GitHub are also welcome.
6 Results and Discussion
6.1 Pt-Loading Effects.
A goal of this model is to provide a physically meaningful explanation for fuel cell losses observed at low Pt loading. To accomplish this task, a set of five parameters including the kinetic rate coefficient for the surface adsorption step (Eq. (10)) kfwd,ads, the concentration-independent portion of the exchange current density io*, membrane resistance RNaf, Nafion O2 thermodynamics (characterized here by the free energy of the absorption reaction in Eq. (8), ), and capture angle are varied to find a set that, when held constant across all Pt loadings, provide a fit with relatively low error. The fitting routine utilizes the L-BFGS-B (i.e., the limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm with bounds) method of the ‘minimize’ package available within SciPy . This algorithm minimizes the chi squared calculated between the model output and data, including variance, while allowing bounds for each parameter. The algorithm is in the family of quasi-Newton methods and uses an estimated hessian matrix to predict directions of improvement for each variable at each step. The optimized values from this fitting procedure are summarized in Tables 4 and 5. Figure 3 shows the best fits for both the core–shell and flooded agglomerate models. The core–shell model is shown using the “mixed” method for Nafion transport properties, as described in Sec. 3.1.3. Details for the properties used in the flooded agglomerate results are discussed below.
As observed in Fig. 3(a), the core–shell model provides a good fit across the range of Pt loadings, with slight discrepancies on the order of 10–20 mV for the lowest Pt loading. These discrepancies, in part, are likely due to the lack of modeling water transport. With the inclusion of liquid water, localized electrode flooding can limit transport at very low Pt loadings and high currents. Although the present study clearly identifies the importance of accurate structure–property relationships for the CL ionomer, incorporating water dynamics will likely improve the fit and change some of the fitted parameter values reported herein. This is the topic of forthcoming work from our group. Figure 3(b) demonstrates that the flooded agglomerate model is not capable of fitting the data. The fit captures trends for the highest two Pt loading values, but additional losses at low Pt are not captured. Moreover, obtaining the reported fit required abandoning the physically informed structure–property relationships developed herein. Using any of the three methods in Sec. 3.1 provides unrealistic limiting currents, related to oxygen transport limitations through the agglomerate. For the fit in Fig. 3, values of roughly 8 × 10−10 m2/s were adopted from Ref. . This value is two orders of magnitude larger than those generated from the structure–property relationship herein.
Comparing the results for the core–shell and flooded agglomerate models, it is apparent that physically realistic values are required to capture Pt loading effects. In the flooded agglomerate model, the high values compensate for the large agglomerate diameter. These large diameters, in turn, are required for valid application of continuum methods to the agglomerate interior. In other words, developing homogenized effective transport properties for strictly radial diffusion through a 150 nm agglomerate is not appropriate, when said agglomerate consists of only 2–5 carbon particles with diameter 50 nm. Such treatment requires larger agglomerates, which in turn must be compensated by artificially large values, for accurate limiting current predictions. However, these high values reduce the sensitivity of the model results to other important structure–property relationships, namely, the conductivity variations with Pt loading predicted by the “mixed” method. A higher affinity for liquid water near Pt surfaces (relative to carbon) correlates to higher ionic conductivities when Pt is more abundant. When is artificially high, the model is no longer sensitive to these variations in σio. Lastly, the core–shell geometry accommodates use of the effective capture angle, providing an additional relationship between the Pt loading and transport losses.
However, at the lowest Pt loading of 0.025 mg/cm2, the core–shell results differ slightly from experimental data, demonstrating a general limitation of these effective microstructure approaches. In a real PEMFC CL, not all of the Pt is as easily accessible to oxygen as the core–shell geometry allows. At very low loading, Pt wedged between carbon cores likely has lower utilization, an effect that the core–shell model does not capture. A reasonable adaptation to hybridize the two microstructures presented here may provide a means to capture this effect. Use of a flooded agglomerate microstructure that neglects Nafion volume fraction and tortuosity effects and implements a transport capture angle for the outermost radial Pt-containing node could result in more realistic limiting currents, with the added benefit of capturing additional concentration losses at very low Pt loadings caused by limited access to oxygen. Alternately, higher fidelity simulations of a single, spatially resolved agglomerate can be used to derive effective transport parameters, which can then serve as inputs into reduced-order models with effective microstructure representations. The development of alternative modeling strategies is the subject of future work from our group.
Owing to its superior fit, the remainder of the results section will present results for only the core–shell model. This microstructure approach will provide a framework for analyzing other model phenomena, including structure–property approaches, Nafion thickness effects, and the impact of surface-mediated kinetic mechanisms.
6.2 Comparing Structure–Property Methods.
Of the three derived structure–property relationships, the “mixed” method provides the greatest sensitivity to Pt loading. Once operating conditions and shell thickness are chosen, the layered and uniform method transport properties are invariant with the Pt loading. The mixed method properties, on the other hand, are an average of those from the other two methods, weighted by the Pt loading. To examine the effects of these relationships, Fig. 4 shows predicted polarization curves for 0.05 mg/cm2 loading, using the layered and uniform structure–property methods. The effective ionic conductivities are 0.9 and 1.25 S/m for the uniform and layered methods, respectively, a difference that has little impact on the overall polarization behavior. Rather, the primary variable impacted by the structure–property method is the limiting current. The diffusion coefficient for oxygen in Nafion decreases by roughly 50% for the layered case, relative to the uniform model, from 6.3 × 10−12 to 3.2 × 10−12 m2/s. This correlates with a limiting current of 2.8 A/cm2 in the layered model compared to 3.6 A/cm2 in the uniform model. The Fig. 4 inset demonstrates that this trend is consistent across Nafion shell thicknesses.
6.3 Shell Thickness Effects.
Figure 5 shows the polarization and power curves of the core–shell model at a Pt loading of 0.05 mg/cm2 with varying ionomer shell thicknesses. In this figure, as the Nafion shell thickness changes, the accompanying change in volume fraction is accommodated by changes in the gas-phase volume fraction. For the 25 nm case, the porosity is 0.1, but as the shell is reduced to 10 nm and 5 nm, the modeled CL porosity changes to approximately 0.5 and 0.6, respectively. This implementation ensures that the Pt surface area and carbon radius do not change so that any changes in polarization behavior are due purely to transport effects. It is clear from the figure that the 25 nm shell provides the highest losses, which are attributed to poor gas transport (due to low ) and large O2 radial diffusive distances through the CL ionomer.
Another interesting result is that the 10 nm shell has a higher maximum power density than the 5 nm shell. This 13% increase in power density is despite lower porosity and the larger diffusive distances for the 10 nm case (which lead to a lower limiting current in the 10 nm case). The implementation of structure–property relationships in this model predicts higher ionic conductivity and oxygen diffusion coefficients for thicker shells due to greater water absorption. The model demonstrates that, within the range of realistic property co-variations, improved transport properties can outweigh correlated deleterious microstructural impacts. This result stands out from other models that assume transport properties which are invariant with Nafion thickness and Pt loading [25,26]. These other models generally use the shell thickness as a fitting parameter to adjust transport resistance. By making the shell thicker, they can capture additional resistance; however, in this model, thicker shells can result in less resistance due to improved transport properties.
To better understand predicted polarization trends as a function of CL Nafion thickness, Figs. 6 and 7 show profiles for the ionomer electric potential, Faradaic current distribution, and O2(Naf) density as a function of CL depth. Results are presented for current density of 1.5 A/cm2 (i.e., where the performance of the 10 nm shell is noticeably better than that of the 5 nm case). Orientation along the x-axis is with the PEM interface at the left and the GDL interface at the right. To explain differences in ionic transport, Fig. 6 shows the electrolyte potential and Faradaic current distribution (fraction of total Faradaic current) throughout the CL depth. This plot demonstrates the impact of higher conductivities due to higher water absorption. The thicker 10 nm shell has higher ionic conductivity, such that utilization of protons is reasonably uniform throughout the CL depth, even near limiting currents. In the thinner 5 nm shell, we observe larger ionomer electric potential gradients due to lower σio, which in turn leads to a greater fraction of the Faradaic current near the electrolyte membrane. The larger Ohmic penalty is therefore associated with higher activation losses, in the 5 nm thick ionomer film case. The combination of losses provides an explanation as to why the thinner ionomer shell predicts lower performance, at this current density.
Looking at the oxygen depth profile in Fig. 7, we observe that the uneven Faradaic current of the 5 nm shell results in a similarly uneven distribution of O2(Naf) concentration gradients throughout the CL depth. The oxygen density of the outermost shell is plotted for both film thicknesses to show that porosity effects do not significantly impact the gas-phase transport. Differences between the outermost and innermost shells represent the radial concentration gradients. For both the 5 and 10 nm film cases, we see that the Faradaic currents are transport-limited near the electrolyte membrane (O2(Naf) densities are near zero at the Pt interface). However, even though the diffusion coefficient is higher for the 10 nm case, the O2(Naf) gradients are more severe farther from the membrane, relative to the 5 nm case. This is due in part to the larger diffusive distance that oxygen travels to reach the reactive Pt sites in the 10 nm case but is also consistent with higher catalyst utilization farther from the membrane, in the 10 nm thickness case. Despite these larger radial diffusive gradients, which lead to a lower limiting current in the 10 nm case, this is more than compensated by the increase in ionic conductivity, leading to improved power density for a range of reasonable operating current densities.
6.4 Reaction Mechanism and Surface Chemistry.
All results to this point are based on the two-step reaction mechanism (Eqs. (10) and (11)). As shown in Fig. 3, this approach captures cell polarization phenomena in the low-current “activation region,” for both the core–shell and flooded agglomerate models over the range of Pt loadings. To explore the impact of detailed surface kinetics, a separate fit to the data was run with the one-step global mechanism (Eq. (9)). Fitting results show that Pt loading effects are not correctly predicted with this mechanism. Figure 8 shows a close-up of the kinetic region from polarization curves created with the core–shell geometry using the two reaction mechanisms. While the two-step mechanism, repeated in Fig. 8(a) for comparison’s sake, matches the data well within this region, Fig. 8(b) shows the one-step mechanism, at its best, results in a grouping of curves centered near the middle of the data. This grouping is caused by poorly predicted activation overpotentials resulting from simplified kinetics.
Equation (13) provides an explanation for this effect. In the Butler–Volmer formulation, the exchange current density is a function of the local activity concentrations [Cac,k] (here, equal to the molar concentrations) of the products and reactants. In the one-step mechanism, the only such species whose activity concentration varies with current density is O2(Naf). The proton and H2O activity concentrations are constant, and so the exchange current density scales with [O2(Naf)]0.5. However, for the two-step mechanism, the intermediate species O(s) and Pt(s) participate directly in the reaction, with . Therefore, as the reactant O(s) is depleted at higher current densities, the product species Pt(s) concentration increases. The exchange current density for this mechanism scales with the ratio [O(s)]/[Pt(s)]. For fast, near-equilibrium kinetics, this ratio will be proportional to [Cac,k], and the two mechanisms are roughly equivalent. However, for slow adsorption kinetics (Eq. (10)) associated with decreased Pt loading, the ratio [O(s)]/[Pt(s)] will drop faster than [O2(Naf)]0.5 with increasing current. With decreasing Pt loading, the kinetic limitation becomes more severe, leading to the Pt loading sensitivity observed in Figs. 3 and 8(a), for the two-step mechanism. To state the obvious: while the global reaction is mathematically convenient and can accurately capture polarization in non-kinetic-limited performance, it is incumbent on the modeler to explore the implications of any kinetic mechanism assumptions. In a recent article, Dickinson and Hinds provide an insightful overview of various approaches to Butler–Volmer kinetics in PEMFCs and highlights the strengths and limitations of a number of common approaches .
To reduce the cost of PEMFCs while maintaining suitable performance, a quantitative and mechanistic understanding of the high observed overpotentials at low Pt loading is required. Two physically based models developed here identify the importance of both transport and kinetic phenomena in understanding and addressing these performance limitations. Physically informed structure–property relationships are hypothesized to describe the transport of ions and dissolved O2 through nano-thin Nafion sheets in the cathode CL. Functionally, these relationships only provide reasonable results for one of the two considered microstructures (the so-called “core–shell” model). With the use of (i) the core–shell model approach, (ii) transport properties based on experimentally observed Nafion structures which are sensitive to Pt loading, and (iii) a two-step, surface-mediated reaction mechanism, cell polarization behavior is matched over a range of low Pt loadings. This result provides a more detailed explanation for these losses compared to generalized resistances used in other models. Capturing Nafion transport properties based on the ionomer structure also provides insight into future design considerations. For example, a slight increase to the CL ionomer thickness predicts increased maximum power densities, due to improved ionic conductivity associated with greater water uptake in the CL Nafion. This work establishes that structure–property relationships based on ionomer thickness have non-negligible impact on the predicted performance and should be considered in future modeling efforts. Future development of new modeling frameworks for effective CL microstructures and incorporation of additional physical phenomena can lead to quantitatively predictive modeling tools to guide the design of advanced PEMFCs [34–36].
This research was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Materials Science and Engineering under Early Career Award #DE-SC0018019, Dr. Pappan Thiyagarajan, Program Manager.