Abstract

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.

1 Introduction

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 [1]. Although FCEVs are currently available, challenges in infrastructure, performance, durability, and cost hinder wide-spread adoption [26].

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 [79] 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 [12] suffer from performance losses. Electrodes with lower Pt loadings exhibit large transport losses that are not yet sufficiently explained in the literature [13]. 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 [16]. 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 [1719]. 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,2426]. 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 [3033], 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 [37], 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 (ρO2(Naf), kgmNafion3) is tracked in the y and r directions. The surface site coverage of the Pt catalyst (θk,surf), when considered, is tracked in y and r (note that θk,surf 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 (ρk,g, kgmgas3) 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 [28], 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:

  1. Isothermal operation.

  2. Water is in equilibrium between the Nafion and gas phases, and the gas remains at a fixed, steady relative humidity (RH).

  3. Solid core–shell and/or flooded agglomerate structures are spherical in shape.

  4. Uniform porosity, microstructure, and Pt distribution throughout the CL.

  5. Electric potential is radially uniform within the spherical solid particles.

  6. 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.

The gas phase is treated as an ideal gas. The mass of each species k is conserved by applying the expression 
dρk,gasdt=1εg(Jk+ANafWks˙k,Naf)
(1)
where εg is the volume fraction of the gas phase and ANaf is the specific area of Nafion–gas interface per unit volume of CL. ρk,gas, Jk, Wk, and s˙k,Naf are the mass density (kgkmgas3), mass flux (kgkmCL2s1), molecular weight (kgmolk1), and net species molar production rate (molkmNafgas2s1) due to reactions at the Nafion–gas phase interface for species k, respectively. The mass flux in the gas phase captures transport via diffusion and convection as 
Jk=ρk,gasKgμPεgτfac,gDkρk,gasYk
(2)
where Kg is the gas-phase permeability (m2), μ is the viscosity (Pa·s), P is the pressure (Pa), τfac,g is the tortuosity factor (estimated using a Bruggeman correlation [42] with an exponent of −0.5), and where Dk and Yk are the diffusion coefficient (m2 s−1) and mass fraction of species k, respectively. Since the model is isothermal, the mass density of each species determines state of the gas phase, via the temperature, total phase density (ρgas=ρk,gas), and mass fractions (Yk=ρk,gas/ρgas).

2.3.2 Electrolyte Mass and Species.

Because water is assumed to be in equilibrium between the electrolyte and gas phase, the only two electrolyte species with variable concentrations are absorbed oxygen and protons. The mass of protons in the electrolyte is held constant based on the equivalent weight of the electrolyte Nafion and the assumption of charge neutrality. For the oxygen, conservation of mass is employed according to 
dρO2(Naf)dt=DO2(Naf)eff(1rddr(r2dρO2(Naf)dr))+WO2(ANafs˙O2,Naf+APts˙O2,Pt)
(3)
where most terms are as previously defined but with O2 replacing the generalized species index k. The new terms for DO2(Naf)eff, s˙O2,Pt, and APt represent the effective diffusion coefficient (m2 s−1) of oxygen absorbed in Nafion, net molar production rate of oxygen (molO2mint2s1) from the reaction occurring at the Nafion–Pt interface, and specific area of Pt per unit CL volume, respectively. Since the gradients in oxygen occur radially due to the Pt reactions within the spherical domains, diffusion of absorbed oxygen in the y direction is not considered. For any radial nodes in the electrolyte that do not contain Pt, s˙O2,Pt is zero. Likewise, the absorption reaction rate, s˙O2,Naf, is nonzero only at the outermost radius where the Nafion–gas interface exists.

2.3.3 Conservation of Charge.

Although protons are assumed here to have constant density within the electrolyte phase, charge conservation is used to calculate the electric potential difference at electrode–electrolyte interfaces. Charge neutrality is enforced, such that the electronic current density into the electrode (iext) equals the ionic current density into the membrane. Moreover, in each discretized CL volume, charge neutrality in the ionomer phase is preserved via the double layer current density idl (Amint2), which compensates for the non-zero sum of the divergence of the ionic current density plus the local Faradaic current density: 
idl=iioAPtiFar
(4)
where the Faradaic current density iFar (Amint2) is defined below in Sec. 2.4 and iio is the ionic current density (AmNafion2). The double layer current determines the double layer potential, ΔΦdl, which is the difference between the electrolyte and electrode material electric potentials within the CL, ΔΦdl = Φelyte,CL − Φed,CL. This is modeled as capacitive in nature: 
ΔΦdl=QdlCdlAdl
(5)
where Qdl is the double layer charge (Cmint2), which leads to the following for the rate of change of ΔΦdl: 
dΔΦdldt=idlCdlAdl
(6)
where Cdl and Adl are the double layer capacitance (Fmint2) and specific double layer area per unit volume of CL, respectively. The cathode potential is chosen as the reference potential (Φca,CL = 0 V), such that ΔΦdl = Φelyte,CL.
The ionic current density, iio, is Ohmic in nature: 
iio=σioεNafτfac,NafΦelyte
(7)
where σio is the ionic conductivity of the ionomer (S m−1) and εNaf and τfac,Naf are the Nafion volume fraction and tortuosity factor in the CL. In actual PEMFC CLs, some of the ionomer volume in εNaf exists not as thin-film coatings or within flooded agglomerates, but rather as discrete, single-phase “globules” with higher ionic conductivity [43]. To compensate for this unmodeled heterogeneity, the tortuosity factor of the ionomer phase within the catalyst layer is set equal to unity.

2.4 Kinetics and Reactions.

In the model, reactions take place at the Nafion–gas and Nafion–Pt interfaces. At the Nafion–gas interface, oxygen moves between the gas and electrolyte phases via absorption: 
O2(g)O2(Naf)
(8)
where (g) and (Naf) denote gas and electrolyte-absorbed phases, respectively. The kinetics for reaction (8) are assumed fast such that the outermost Nafion shell remains near equilibrium with the gas phase.
At the Pt surface, two reaction mechanisms are considered: a one-step global reaction and a two-step surface-mediated mechanism. The one-step mechanism is written as 
2H+(Naf)+12O2(Naf)+2e(Pt)H2O(Naf)
(9)
For the two-step mechanism, charge transfer proceeds through an intermediate oxygen species adsorbed on the Pt surface, followed by a global, two-electron charge transfer reaction: 
O2(Naf)+2Pt(s)2O(s)
(10)
 
O(s)+2H+(Naf)+2e(ca)H2O(Naf)+Pt(s)
(11)
The label (s) denotes a Pt surface species and e(ca) is an electron in a cathode electronic conductor phase (Pt or C). Thermochemical reactions (e.g., reactions (8) and (10)) assume concentration-dependent reversible mass-action kinetics: 
s˙k,j=νk,j(kfwd,j[Cac,k]νk,jkrev,j[Cac,k]νk,j)
(12)
where νk,j, νk,j, and νk,j are the net, forward, and reverse stoichiometric coefficients for species k in reaction j, and where νk,j=νk,jνk,j. [Cac,k] is the “activity concentration” of species k. This study assumes ideal thermodynamics (activity coefficients equal unity), such that the activity concentration equals the molar concentration, [Ck] (mol m−3 for Nafion species, mol m−2 for surface species), and kfwd,j and krev,j are the forward and reverse reaction rate coefficients, whose ratio is dictated by thermodynamic reversibility. For the charge transfer reactions in either mechanism (e.g., reactions (9) and (11)), a concentration-dependent Butler–Volmer expression calculates iFar: 
iFar,j=i,j*[Cac,k]νk,j[eβnelecFη/R¯Te(1β)nelecFη/R¯T]
(13)
where nelec is the elementary electrical charge transferred by the charge transfer reaction, β is the charge transfer symmetry factor, T is the temperature (K), and R¯ is the universal gas constant (J mol−1 K−1). The composition-independent portion of the exchange current density i,j* for charge transfer reaction j is calculated using Arrhenius parameters. The molar production rate for a species k associated with a charge transfer reaction j is 
s˙k,j=νk,jiFar,jnelecF
(14)
The total species production rates s˙k in Eqs. (1) and (3) come from summing s˙k,j over all j reactions for a given interface.

2.5 Surface Coverages.

When using the two-step mechanism, the surface coverages of species k, θk evolve with time: 
dθkdt=1Γs˙k,Pt
(15)
where Γ is the site density on the Pt surface (molmPt2), and [Ck]=Γθk for any surface species k. When using the one-step mechanism, the intermediate species that cover the Pt surface are not resolved, and Eq. (15) is neglected.

2.6 Membrane.

As a bulk material, Nafion has been studied extensively and is well documented. The experimentally determined area-specific resistance RNaf (Ω m2) can be used to describe the Ohmic losses across the membrane as 
ΔΦOhmic=Φelyte,anΦelyte,ca=iextRNaf
(16)
where RNaf is on the order of what has been experimentally measured [44].

2.7 Anode.

Although the anode is not explicitly modeled here, it must still be considered, to compare the model output to experimental results. We therefore assume that the anode operates at equilibrium. For the anode reaction: 
H2(Naf)2H+(Naf)+2e(an)
(17)
the Gibbs free energy of reaction Δgrxn is used to calculate the equilibrium potential 
ΔΦaneq=Φelyte,anΦan=ΔgrxnnelecF
(18)

2.8 Cell Voltage.

The cell voltage at each current density is the difference between the cathode and anode electric potentials: 
VCell=ΦcaΦan=ΦcaΦelyte,ca(ΦanΦelyte,ca)
(19)
Substituting in Eqs. (16) and (18) and setting Φca = 0, we get 
Vcell=ΔΦaneqΦelyte,CLiextRNaf
(20)

2.9 Boundary Conditions.

Boundary conditions assume a fixed temperature, pressure, and composition Yk at the gas channel/GDL interface. Between the CL and membrane, two boundary conditions are applied. For all gas phase species k, there is no flux in the y direction into the membrane: 
Jk|mem=0
(21)
and the ionic current equals the external current: 
iio|mem=iext
(22)
within the Nafion–Pt–Carbon agglomerates, the radial boundary conditions are dependent upon the chosen CL structure. All models assume a radial diffusion flux of zero at the Nafion/gas interface. In the core–shell model, zero diffusive flux is assumed at the Nafion interface with Carbon/Pt, which occurs at a single radius, rcarbon. For the flooded agglomerate model, the Pt–Nafion and carbon–Nafion interfaces exist at all radial locations except for the outermost shell. This moves the boundary condition of zero diffusive flux to r = 0, which is treated as a symmetry plane.

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. [29] 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:

  1. Empirical relationships for conductivity as a function of hydration λ from bulk Nafion samples [49,50].

  2. Ionic conductivity data from Nafion films at relevant thicknesses [20].

  3. Composition depth profiles of similarly thin Nafion, measured via NR [22].

Because the existence and transport properties associated with layered structures in PEMFC CLs are as yet unconfirmed, we propose here three hypothetical structure–property relationships for thin film Nafion and explore their impacts on predicted polarization behavior. These three models are illustrated in Fig. 2. In our previous publication [22], we correlated structural data from NR with thin film conductivity values measured by Paul and Karan [20] to infer “transport scaling factors” for a range of Nafion layers, representing both the anisotropic interfacial multilayers and the thicker isotropic layers. These scaling factors represent reductions in ionic conductivity for a layer, beyond what would be predicted based on the local water content. In the present model, the local ionic conductivity is obtained by multiplying these scaling factors by conductivity values calculated according to the local layer water content, λ: 
σl=σbulk(λ)fl
(23)
where σl and fl are the ionic conductivity and scaling factor for layer l, respectively, and σbulk(λ) is the reference ionic conductivity of a bulk Nafion membrane with the same λ as the layer [49].
For O2 diffusion coefficients, experimentally measured DO2 values for thick films [50] are scaled by the local water content, relative to that in bulk membranes at similar RH, and are then multiplied by the transport scaling factor: 
DO2,Naf,l=DO2,Nafλlλfl
(24)
where DO2,Naf is the experimentally measured diffusion coefficient for thick Nafion membranes, λ° is the water uptake of the reference membrane, and λl and fl are the water uptake and scaling factor of layer l, respectively.

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. [21], 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 [22]. 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 [22].

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, DO2,Naf,l 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 [53]. 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. [54].

4.2 Cantera.

The open-source software library Cantera is used to manage the model’s thermodynamic, kinetic, and gas-phase transport parameter calculations [55]. 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 [56].

5 Accessing the Model

To advance open science practices related to accessibility, transparency, and repeatability, the model code is made freely available via GitHub [57]. 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), ΔgO2,abs), 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 [53]. 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, DO2 values of roughly 8 × 10−10 m2/s were adopted from Ref. [25]. 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 DO2 values are required to capture Pt loading effects. In the flooded agglomerate model, the high DO2 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 DO2 values, for accurate limiting current predictions. However, these high DO2 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 DO2 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 εg) 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 θO(s)=1θPt(s). 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 [58].

7 Conclusions

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 [3436].

Footnote

Acknowledgment

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.

References

References
1.
Thomas
,
C. E.
,
2009
, “
Fuel Cell and Battery Electric Vehicles Compared
,”
Int. J. Hydrogen Energy
,
34
(
15
), pp.
6005
6020
. 10.1016/j.ijhydene.2009.06.003
2.
Apostolou
,
D.
, and
Xydis
,
G.
,
2019
, “
A Literature Review on Hydrogen Refuelling Stations and Infrastructure. Current Status and Future Prospects
,”
Renew. Sust. Energy Rev.
,
113
, p.
109292
. 10.1016/j.rser.2019.109292
3.
Peighambardoust
,
S. J.
,
Rowshanzamir
,
S.
, and
Amjadi
,
M.
,
2010
, “
Review of the Proton Exchange Membranes for Fuel Cell Applications
,”
Int. J. Hydrogen Energy
,
35
(
17
), pp.
9349
9384
. 10.1016/j.ijhydene.2010.05.017
4.
Alaswad
,
A.
,
Baroutaji
,
A.
,
Achour
,
H.
,
Carton
,
J.
,
Makky
,
A. A.
, and
Olabi
,
A. G.
,
2016
, “
Developments in Fuel Cell Technologies in the Transport Sector
,”
Int. J. Hydrogen Energy
,
41
(
37
), pp.
16499
16508
. 10.1016/j.ijhydene.2016.03.164
5.
Zhang
,
T.
,
Wang
,
P.
,
Chen
,
H.
, and
Pei
,
P.
,
2018
, “
A Review of Automotive Proton Exchange Membrane Fuel Cell Degradation Under Start–Stop Operating Condition
,”
Appl. Energy
,
223
, pp.
249
262
. 10.1016/j.apenergy.2018.04.049
6.
Eichman
,
J.
,
Brouwer
,
J.
, and
Samuelsen
,
S.
,
2010
, “
Exploration and Prioritization of Fuel Cell Commercialization Barriers for Use in the Development of a Fuel Cell Roadmap for California
,”
ASME J. Fuel Cell Sci. Technol.
,
7
(
5
), p.
1189
. 10.1115/1.4000689
7.
Sakai
,
K.
,
Sato
,
K.
,
Mashio
,
T.
,
Ohma
,
A.
,
Yamaguchi
,
K.
, and
Shinohara
,
K.
,
2009
, “Analysis of Reactant Gas Transport in Catalyst Layers; Effect of Pt-Loadings,”
Proton Exchange Membrane Fuel Cells 9, Vol. 25 of ECS Transactions
,
Electrochem. Soc. Inc.
, pp.
1193
1201
.
8.
Ohma
,
A.
,
Mashio
,
T.
,
Sato
,
K.
,
Iden
,
H.
,
Ono
,
Y.
,
Sakai
,
K.
,
Akizuki
,
K.
,
Takaichi
,
S.
, and
Shinohara
,
K.
,
2011
, “
Analysis of Proton Exchange Membrane Fuel Cell Catalyst Layers for Reduction of Platinum Loading At Nissan
,”
Electrochim. Acta
,
56
(
28
), pp.
10832
10841
. 10.1016/j.electacta.2011.04.058
9.
Leimin
,
X.
,
Shijun
,
L.
,
Lijun
,
Y.
, and
Zhenxing
,
L.
,
2009
, “
Investigation of a Novel Catalyst Coated Membrane Method to Prepare Low-Platinum-Loading Membrane Electrode Assemblies for PEMFCs
,”
Fuel Cells
,
9
(
2
), pp.
101
105
. 10.1002/fuce.200800114
10.
Othman
,
R.
,
Dicks
,
A. L.
, and
Zhu
,
Z.
,
2012
, “
Non Precious Metal Catalysts for the PEM Fuel Cell Cathode
,”
Int. J. Hydrogen Energy
,
37
(
1
), pp.
357
372
. 10.1016/j.ijhydene.2011.08.095
11.
Bonham
,
D.
,
Choi
,
J.-Y.
,
Kishimoto
,
T.
, and
Ye
,
S.
,
2019
, “
Integrating PGM-Free Catalysts Into Catalyst Layers and Proton Exchange Membrane Fuel Cell Devices
,”
Adv. Mater.
,
31
(
31, SI
), p.
1804846(1)
.
12.
The US Department of Energy (DOE)
, “
DOE Technical Targets for Polymer Electrolyte Membrane Fuel Cell Components
,” https://www.energy.gov/eere/fuelcells/doe-technical-targets-polymer-electrolyte-membrane-fuel-cell-components
13.
Weber
,
A. Z.
, and
Kusoglu
,
A.
,
2014
, “
Unexplained Transport Resistances for Low-Loaded Fuel-Cell Catalyst Layers
,”
J. Mater. Chem. A
,
2
(
41
), pp.
17207
17211
. 10.1039/C4TA02952F
14.
Hao
,
L.
,
Moriyama
,
K.
,
Gu
,
W.
, and
Wang
,
C.-Y.
,
2015
, “
Modeling and Experimental Validation of Pt Loading and Electrode Composition Effects in PEM Fuel Cells
,”
J. Electrochem. Soc.
,
162
(
8
), pp.
F854
F867
. 10.1149/2.0221508jes
15.
Kongkanand
,
A.
,
Gu
,
W.
, and
Mathias
,
M. F.
,
2017
, “Proton-Exchange Membrane Fuel Cells With Low-Pt Content,”
Encyclopedia of Sustainability Science and Technology
,
Meyers
,
R. A.
, ed.,
Springer
,
New York
, pp.
1
20
.
16.
Liu
,
H.
,
Epting
,
W. K.
, and
Litster
,
S.
,
2015
, “
Gas Transport Resistance in Polymer Electrolyte Thin Films on Oxygen Reduction Reaction Catalysts
,”
Langmuir
,
31
(
36
), pp.
9853
9858
. 10.1021/acs.langmuir.5b02487
17.
Li
,
H.
,
Tang
,
Y.
,
Wang
,
Z.
,
Shi
,
Z.
,
Wu
,
S.
,
Song
,
D.
,
Zhang
,
J.
,
Fatih
,
K.
,
Zhang
,
J.
,
Wang
,
H.
,
Liu
,
Z.
,
Abouatallah
,
R.
, and
Mazza
,
A.
,
2008
, “
A Review of Water Flooding Issues in the Proton Exchange Membrane Fuel Cell
,”
J. Power Sources
,
178
(
1
), pp.
103
117
. 10.1016/j.jpowsour.2007.12.068
18.
Luo
,
Z.
,
Chang
,
Z.
,
Zhang
,
Y.
,
Liu
,
Z.
, and
Li
,
J.
,
2010
, “
Electro-Osmotic Drag Coefficient and Proton Conductivity in Nafion® Membrane for PEMFC
,”
Int. J. Hydrogen Energy
,
35
(
7
), pp.
3120
3124
. 10.1016/j.ijhydene.2009.09.013
19.
Nonoyama
,
N.
,
Okazaki
,
S.
,
Weber
,
A. Z.
,
Ikogi
,
Y.
, and
Yoshida
,
T.
,
2011
, “
Analysis of Oxygen-Transport Diffusion Resistance in Proton-Exchange-Membrane Fuel Cells
,”
J. Electrochem. Soc.
,
158
(
4
), p.
B416
. 10.1149/1.3546038
20.
Paul
,
D. K.
, and
Karan
,
K.
,
2014
, “
Conductivity and Wettability Changes of Ultrathin Nafion Films Subjected to Thermal Annealing and Liquid Water Exposure
,”
J. Phys. Chem. C
,
118
(
4
), pp.
1828
1835
. 10.1021/jp410510x
21.
Paul
,
D. K.
,
McCreery
,
R.
, and
Karan
,
K.
,
2014
, “
Proton Transport Property in Supported Nafion Nanothin Films by Electrochemical Impedance Spectroscopy
,”
J. Electrochem. Soc.
,
161
(
14
), pp.
F1395
F1402
. 10.1149/2.0571414jes
22.
DeCaluwe
,
S. C.
,
Baker
,
A. M.
,
Bhargava
,
P.
,
Fischer
,
J. E.
, and
Dura
,
J. A.
,
2018
, “
Structure–Property Relationships at Nafion Thin-Film Interfaces: Thickness Effects on Hydration and Anisotropic Ion Transport
,”
Nano Energy
,
46
, pp.
91
100
. 10.1016/j.nanoen.2018.01.008
23.
DeCaluwe
,
S. C.
,
Kienzle
,
P. A.
,
Bhargava
,
P.
,
Baker
,
A. M.
, and
Dura
,
J. A.
,
2014
, “
Phase Segregation of Sulfonate Groups in Nafion Interface Lamellae, Quantified Via Neutron Reflectometry Fitting Techniques for Multi-Layered Structures
,”
Soft Matter
,
10
(
31
), pp.
5763
5776
. 10.1039/C4SM00850B
24.
Harvey
,
D.
,
Pharoah
,
J. G.
, and
Karan
,
K.
,
2008
, “
A Comparison of Different Approaches to Modelling the PEMFC Catalyst Layer
,”
J. Power Sources
,
179
(
1
), pp.
209
219
. 10.1016/j.jpowsour.2007.12.077
25.
Sun
,
W.
,
Peppley
,
B. A.
, and
Karan
,
K.
,
2005
, “
An Improved Two-Dimensional Agglomerate Cathode Model to Study the Influence of Catalyst Layer Structural Parameters
,”
Electrochim. Acta
,
50
(
16–17
), pp.
3359
3374
. 10.1016/j.electacta.2004.12.009
26.
Yoon
,
W.
, and
Weber
,
A. Z.
,
2011
, “
Modeling Low-Platinum-Loading Effects in Fuel-Cell Catalyst Layers
,”
J. Electrochem. Soc.
,
158
(
8
), pp.
B1007
B1018
. 10.1149/1.3597644
27.
Weber
,
A. Z.
, and
Newman
,
J.
,
2004
, “
Modeling Transport in Polymer-Electrolyte Fuel Cells
,”
Chem. Rev.
,
104
(
10
), pp.
4679
4726
. 10.1021/cr020729l
28.
Weber
,
A. Z.
,
Borup
,
R. L.
,
Darling
,
R. M.
,
Das
,
P. K.
,
Dursch
,
T. J.
,
Gu
,
W.
,
Harvey
,
D.
,
Kusoglu
,
A.
,
Litster
,
S.
,
Mench
,
M. M.
,
Mukundan
,
R.
,
Owejan
,
J. P.
,
Pharoah
,
J. G.
,
Secanell
,
M.
, and
Zenyuk
,
I. V.
,
2014
, “
A Critical Review of Modeling Transport Phenomena in Polymer-Electrolyte Fuel Cells
,”
J. Electrochem. Soc.
,
161
(
12
), pp.
F1254
F1299
. 10.1149/2.0751412jes
29.
Zenyuk
,
I. V.
,
Das
,
P. K.
, and
Weber
,
A. Z.
,
2016
, “
Understanding Impacts of Catalyst-Layer Thickness on Fuel-Cell Performance Via Mathematical Modeling
,”
J. Electrochem. Soc.
,
163
(
7
), pp.
F691
F703
. 10.1149/2.1161607jes
30.
Zhou
,
J.
,
Putz
,
A.
, and
Secanell
,
M.
,
2017
, “
A Mixed Wettability Pore Size Distribution Based Mathematical Model for Analyzing Two-Phase Flow in Porous Electrodes I. Mathematical Model
,”
J. Electrochem. Soc.
,
164
(
6
), pp.
F530
F539
. 10.1149/2.0381706jes
31.
Das
,
P. K.
,
Li
,
X.
, and
Liu
,
Z.-S.
,
2010
, “
Analysis of Liquid Water Transport in Cathode Catalyst Layer of PEM Fuel Cells
,”
Int. J. Hydrogen Energy
,
35
(
6
), pp.
2403
2416
. 10.1016/j.ijhydene.2009.12.160
32.
Zhang
,
G.
,
Fan
,
L.
,
Sun
,
J.
, and
Jiao
,
K.
,
2017
, “
A 3d Model of PEMFC Considering Detailed Multiphase Flow and Anisotropic Transport Properties
,”
Int. J. Heat Mass Transfer
,
115
, pp.
714
724
. 10.1016/j.ijheatmasstransfer.2017.07.102
33.
Dujc
,
J.
,
Forner-Cuenca
,
A.
,
Marmet
,
P.
,
Cochet
,
M.
,
Vetter
,
R.
,
Schumacher
,
J. O.
, and
Boillat
,
P.
,
2018
, “
Modeling the Effects of Using Gas Diffusion Layers With Patterned Wettability for Advanced Water Management in Proton Exchange Membrane Fuel Cells
,”
ASME J. Electrochem. Energy Convers. Storage
,
15
(
2
), p.
021001
. 10.1115/1.4038626
34.
Xing
,
L.
,
Shi
,
W.
,
Das
,
P. K.
, and
Scott
,
K.
,
2017
, “
Inhomogeneous Distribution of Platinum and Ionomer in the Porous Cathode to Maximize the Performance of a PEM Fuel Cell
,”
AIChE J.
,
63
(
11
), pp.
4895
4910
. 10.1002/aic.15826
35.
Xing
,
L.
,
Wang
,
Y.
,
Das
,
P. K.
,
Scott
,
K.
, and
Shi
,
W.
,
2018
, “
Homogenization of Current Density of PEM Fuel Cells by In-Plane Graded Distributions of Platinum Loading and GDL Porosity
,”
Chem. Eng. Sci.
,
192
, pp.
699
713
. 10.1016/j.ces.2018.08.029
36.
Xing
,
L.
,
Xu
,
Y.
,
Das
,
P. K.
,
Mao
,
B.
,
Xu
,
Q.
,
Su
,
H.
,
Wu
,
X.
, and
Shi
,
W.
,
2019
, “
Numerical Matching of Anisotropic Transport Processes in Porous Electrodes of Proton Exchange Membrane Fuel Cells
,”
Chem. Eng. Sci.
,
195
, pp.
127
140
. 10.1016/j.ces.2018.11.034
37.
He
,
C.
,
Desai
,
S. H.
,
Brown
,
G. G.
, and
Bollepalli
,
S. B.
,
2005
, “
PEM Fuel Cell Catalysts: Cost, Performance, and Durability
,”
Electrochem. Soc. Interface
,
14
(
3
), pp.
41
44
.
38.
Cong
,
Z.
, and
Xinghu
,
L.
,
2010
, “
A New One Dimensional Steady State Model for PEM Fuel Cell
,”
25TH World Battery, Hybrid And Fuel Cell Electric Vehicle Symposium and Exhibition Proceedings
,
Beijing, China
, Vols.
1–2
, pp.
956
960
.
39.
Vasilyev
,
A.
,
Andrews
,
J.
,
Jackson
,
L. M.
,
Dunnett
,
S. J.
, and
Davies
,
B.
,
2017
, “
Component-Based Modelling of PEM Fuel Cells With Bond Graphs
,”
Int. J. Hydrogen Energy
,
42
(
49
), pp.
29406
29421
. 10.1016/j.ijhydene.2017.09.004
40.
Pasaogullari
,
U.
, and
Wang
,
C. Y.
,
2004
, “
Liquid Water Transport in Gas Diffusion Layer of Polymer Electrolyte Fuel Cells
,”
J. Electrochem. Soc.
,
151
(
3
), p.
A399
. 10.1149/1.1646148
41.
Santamaria
,
A. D.
,
Das
,
P. K.
,
MacDonald
,
J. C.
, and
Weber
,
A. Z.
,
2014
, “
Liquid–Water Interactions With Gas-Diffusion-Layer Surfaces
,”
J. Electrochem. Soc.
,
161
(
12
), pp.
F1184
F1193
. 10.1149/2.0321412jes
42.
Bruggeman
,
G.
,
1935
, “
Calculation of Various Physics Constants in Heterogeneous Substances I Dielectricity Constants and Conductivity of Mixed Bodies From Isotropic Substances
,”
Annalen der Physik
,
416
(
7
), pp.
636
664
. 10.1002/andp.19354160705
43.
Paul
,
D. K.
,
Karan
,
K.
,
Docoslis
,
A.
,
Giorgi
,
J. B.
, and
Pearce
,
J.
,
2013
, “
Characteristics of Self-Assembled Ultrathin Nafion Films
,”
Macromol.
,
46
(
9
), pp.
3461
3475
. 10.1021/ma4002319
44.
Slade
,
S.
,
Campbell
,
S.
,
Ralph
,
T.
, and
Walsh
,
F.
,
2002
, “
Ionic Conductivity of an Extruded Nafion 1100 EW Series of Membranes
,”
J. Electrochem. Soc.
,
149
(
12
), pp.
A1556
A1564
. 10.1149/1.1517281
45.
Kusoglu
,
A.
,
Dursch
,
T. J.
, and
Weber
,
A. Z.
,
2016
, “
Nanostructure/Swelling Relationships of Bulk and Thin-Film PFSA Ionomers
,”
Adv. Funct. Mater.
,
26
(
27
), pp.
4961
4975
. 10.1002/adfm.201600861
46.
Eastman
,
S. A.
,
Kim
,
S.
,
Page
,
K. A.
,
Rowe
,
B. W.
,
Kang
,
S.
,
DeCaluwe
,
S. C.
,
Dura
,
J. A.
,
Soles
,
C. L.
, and
Yager
,
K. G.
,
2013
, “
Correction to Effect of Confinement on Structure, Water Solubility, and Water Transport in Nafion Thin Films
,”
Macromolecules
,
46
(
2
), pp.
571
571
. 10.1021/ma302140d
47.
Kusoglu
,
A.
,
Kwong
,
A.
,
Clark
,
K. T.
,
Gunterman
,
H. P.
, and
Weber
,
A. Z.
,
2012
, “
Water Uptake of Fuel-Cell Catalyst Layers
,”
J. Electrochem. Soc.
,
159
(
9
), pp.
F530
F535
. 10.1149/2.031209jes
48.
Dura
,
J. A.
,
Murthi
,
V. S.
,
Hartman
,
M.
,
Satija
,
S. K.
, and
Majkrzak
,
C. F.
,
2009
, “
Multilamellar Interface Structures in Nafion
,”
Macromolecules
,
42
(
13
), pp.
4769
4774
. 10.1021/ma802823j
49.
Yadav
,
R.
, and
Fedkiw
,
P. S.
,
2012
, “
Analysis of EIS Technique and Nafion 117 Conductivity as a Function of Temperature and Relative Humidity
,”
J. Electrochem. Soc.
,
159
(
3
), pp.
B340
B346
. 10.1149/2.104203jes
50.
Sethuraman
,
V. A.
,
Khan
,
S.
,
Jur
,
J. S.
,
Haug
,
A. T.
, and
Weidner
,
J. W.
,
2009
, “
Measuring Oxygen, Carbon Monoxide and Hydrogen Sulfide Diffusion Coefficient and Solubility in Nafion Membranes
,”
Electrochim. Acta
,
54
(
27
), pp.
6850
6860
. 10.1016/j.electacta.2009.06.068
51.
Shrivastava
,
U. N.
,
Fritzsche
,
H.
, and
Karan
,
K.
,
2018
, “
Interfacial and Bulk Water in Ultrathin Films of Nafion, 3m PFSA, and 3m PFIA Ionomers on a Polycrystalline Platinum Surface
,”
Macromolecules
,
51
(
23
), pp.
9839
9849
. 10.1021/acs.macromol.8b01240
52.
Murthi
,
V. S.
,
Dura
,
J.
,
Satija
,
S.
, and
Majkrzak
,
C.
,
2008
, “
Water Uptake and Interfacial Structural Changes of Thin Film Nafion® Membranes Measured by Neutron Reflectivity for PEM Fuel Cells
,”
ECS Trans.
,
16
(
2
), pp.
1471
1485
.
53.
Virtanen
,
P.
,
Gommers
,
R.
,
Oliphant
,
E.
,
Haberland
,
M.
,
Reddy
,
T.
,
Cournapeau
,
D.
,
Burovski
,
E.
,
Peterson
,
P.
,
Weckesser
,
W.
,
Bright
,
J.
,
van der Walt
,
S.
,
Brett
,
M.
,
Wilson
,
J.
,
Millman
,
K.
,
Mayorov
,
N.
,
Nelson
,
A.
,
Jones
,
E.
,
Kern
,
R.
,
Larson
,
E.
,
Carey
,
C.
,
Polat
,
I.
,
Feng
,
Y.
,
Moore
,
E.
,
VanderPlas
,
J.
,
Lexalde
,
D.
,
Perktold
,
J.
,
Cimrman
,
R.
,
Henriksen
,
I.
,
Quintero
,
E.
,
Harris
,
C.
,
Archibald
,
A.
,
Ribeiro
,
A.
,
Pedregosa
,
F.
, and
van Mulbregt
,
P.
,
2020
, “
SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python
,”
Nature Methods
.
54.
Owejan
,
J. P.
,
Owejan
,
J. E.
, and
Gu
,
W.
,
2013
, “
Impact of Platinum Loading and Catalyst Layer Structure on PEMFC Performance
,”
J. Electrochem. Soc.
,
160
(
8
), pp.
F824
F833
. 10.1149/2.072308jes
55.
Goodwin
,
D. G.
,
Speth
,
R. L.
,
Moffat
,
H. K.
, and
Weber
,
B. W.
,
2018
, “
Cantera: An Object-Oriented Software Toolkit for Chemical Kinetics, Thermodynamics, and Transport Processes
,”
Version 2.4.0
, https://www.cantera.org
56.
DeCaluwe
,
S. C.
,
Weddle
,
P. J.
,
Zhu
,
H.
,
Colclasure
,
A. M.
,
Bessler
,
W. G.
,
Jackson
,
G. S.
, and
Kee
,
R. J.
,
2018
, “
On the Fundamental and Practical Aspects of Modeling Complex Electrochemical Kinetics and Transport
,”
J. Electrochem. Soc.
,
165
(
13
), pp.
E637
E658
. 10.1149/2.0241813jes
57.
Randall
,
C. R.
, and
DeCaluwe
,
S. C.
,
2019
, “
Pseudo-2D PEMFC Catalyst Layer Model
,”
Version 1.1
, https://github.com/coresresearch/p2d_pemfc
58.
Dickinson
,
E. J. F.
, and
Hinds
,
G.
,
2019
, “
The Butler–Volmer Equation for Polymer Electrolyte Membrane Fuel Cell (PEMFC) Electrode Kinetics: A Critical Discussion
,”
J. Electrochem. Soc.
,
166
(
4
), pp.
F221
F231
. 10.1149/2.0361904jes