Communication—Modeling Polymer-Electrolyte Fuel-Cell Agglomerates with Double-Trap Kinetics

A new semi-analytical agglomerate model is presented for polymer-electrolyte fuel-cell cathodes. The model uses double-trap kinetics for the oxygen-reduction reaction, which can capture the observed potential-dependent coverage and Tafel-slope changes. An iterative semi-analytical approach is used to obtain reaction rate constants from the double-trap kinetics, oxygen concentration at the agglomerate surface, and overall agglomerate reaction rate. The analytical method can predict reaction rates within 2% of the numerically simulated values for a wide range of oxygen concentrations, overpotentials, and agglomerate sizes, while saving simulation time compared to a fully numerical approach.

A common assumption in the mathematical modeling of polymerelectrolyte fuel-cell (PEFC) cathodes is the existence of small agglomerates, first proposed by Giner and Hunter 1 and since then extensively used in literature. [2][3][4][5][6][7][8][9] Agglomerate-based models provide a way to incorporate some microstructural details within the macroscale modeling framework necessary for full-cell simulations. The analytical agglomerate models in literature typically assume a first-order Tafel kinetics; however, the literature on oxygen-reduction-reaction (ORR) order is differing. [10][11][12] Furthermore, the coverage of intermediate species is known to be potential and concentration dependent, 13,14 resulting in a nonuniform Tafel slope. 15 Recently, coverage factors and non-first order kinetics have been incorporated with Tafel kinetics to improve its accuracy; 7 however, the model still assumes a single rate-determining step (RDS) and a single Tafel slope.
Double-trap (DT) kinetics 13 takes into account multiple reaction pathways and the effect of different species coverages. The model is able to reproduce the doubling of the Tafel slope 14 observed experimentally. While Moore, et al. 9 used a numerical agglomerate model with different order Tafel kinetics and DT kinetics, only Zenyuk, et al. 16 have used a simple analytical approach to include the DT kinetics within an agglomerate model. This article is aimed at developing an analytical agglomerate model based on DT kinetics, thereby enabling the benefit of both models. The analytical results are compared to a numerical model to access its accuracy.

Theory
Double-trap kinetics.-The cathode ORR is modeled as a dualpathway reaction 13,14 Where the steps are RA: reductive adsorption, RD: reductive desorption, DA: dissociative adsorption, and RT: reductive transition. The overall kinetic current is given as the combination of current from the two pathways.
where j * is a reference scaling factor used to set the scale for activation energy with its value set to 1000 mA/cm 2 ; G * R D , and G * −R D are the potential dependent energies; and θ O H and θ Pt are potential and concentration dependent fractional coverage factors of OH ad and free Pt surface respectively. The detailed methodology to obtain these variables is discussed elsewhere in literature. 14 Agglomerate numerical model.-A spherical, ionomer-filled agglomerate with a surrounding ionomer film is assumed, with an illustration given in Figure 1. The mathematical model of the agglomerate is similar to the model by Yoon and Weber. 7 The agglomerate is assumed to have uniform ionic and electric potential, and the oxygen transport is assumed to be Fickian. The reaction rate in core is calculated using DT kinetics, and no reaction is assumed in the surrounding film. Details of the model are provided in the supplementary material Section 1. From the numerical solution, average reaction rate in core is estimated.
Agglomerate analytical model.-Using a numerical agglomerate model in a complete cathode simulation is computationally expensive. An analytical or semi-analytical model is much preferred. The average reaction rate in the core can be expressed as 17 where R sur f O 2 is the reaction rate at the surface of agglomerate core (r = r agg ). The effectiveness factor E r is given as 18 where φ L is the Thiele modulus. Thiele modulus for a simple m th order reaction, or a Langmuir type of reaction has been provided in  literature; 19 however, the DT kinetics does not conform to either of these forms. If the DT kinetics can be approximated with an m th order reaction, the Thiele modulus can be obtained and the equation above used.

Results and Discussion
Fitting DT to m th order kinetics.-The DT reaction rate can be fitted to an m th order Tafel kinetics 7 as follows.
The reference concentration c re f O 2 is assumed to be the same as in DT kinetics. Taking the log of Eq. 6 yields which is linear with respect to log c O 2 . Using Eq. 3, i r xn is computed for different overpotential (η) and oxygen concentration (c O 2 ). Figure 2 shows the plot of i r xn against c O 2 at different overpotentials, and its fit to Eq. 6. The Tafel equation results in a good fit to the DT current. The slope of each line is estimated to find the reaction order m. The reaction order is a function of overpotential and varies between 0.45 and 0.52 (Figure 1, supplementary material). A half reaction order is assumed for further analysis. The intercept of each fitted line is calculated to get A(η). The term m log c re f O 2 is removed from A(η) to eliminate the c O 2 dependence (B(η) ≡ log i 0 − αF RT η). Figure 3 shows the profile of A(η) and B(η). The slope of B(η) is calculated to obtain α and the Tafel slope. The intercept of B(η) is calculated to compute the exchange current density, i 0 . As expected with DT, two distinct zones are observed with different Tafel slopes and exchange current densities. Table I gives the fitted kinetic parameters for the two Tafel zones.  Due to the implicit nature of the equations, the following iterative procedure is designed to obtain surface concentrations and reaction rates: 1. Obtain an initial estimate of rate constant from fitted Tafel kinetics parameters 2. Calculate φ L using the general Thiele equation, 19 Since c sur f O 2 is not yet known, assume first-order reaction for initial estimate  3. Calculate an adjusted effectiveness factor (assume a first order for initial estimate) where f is a fitted correction factor to account for the deviation from where i r xn is obtained from Eq. 3.
6. Repeat steps 2 and 3, with correct reaction order and surface concentration. 7. Find the average reaction rate in core (first DT estimate) [14] 8. Repeat steps 4 to 7 to get refined value ofR core O 2 . Using the above outlined procedure, the average agglomerate current was calculated and compared to fully numerical simulations (supplementary material, Figure 2). After two iterations, the DT estimates are within 2% of the numerically estimated values.
To ensure that the model can be used over different ranges of agglomerate sizes and film thicknesses, the effect of these parameters was studied. While changing the agglomerate radius keeps the accuracy of analytical estimate within 2% of the numerical solution (supplementary material, Figure 3), increasing the ionomer film thickness reduces the accuracy of the DT estimate (supplementary material, Figure 4). Although the analytical model is accurate at most film thicknesses; high overpotential and larger film thickness results in inaccurate predictions, thereby requiring a third iteration. Large thicknesses and high reaction rates limit the surface concentration.

Use of DT based agglomerate model for MEA analysis.-To
elucidate the importance of DT kinetics in PEFC simulation, the newly developed semi-analytical model was used in a cathode model of a 2-D cross-section of fuel cell. The cell was modeled in COMSOL using the approach of Zenyuk, et al. 16 (Details in supplementary section 4). Polarization curves are obtained using the model and compared to polarization curves obtained using Tafel kinetics. Figure 4 shows the comparison of DT and Tafel based polarization curves for a cell which is only limited by kinetic losses and provided with 100% O 2 feed. The DT based model has lower current than Tafel due to the lower reaction order. Due to the doubling of slope at higher currents for DT, the difference between DT and Tafel increases at higher currents. For an air feed, similar differences are observed between DT and Tafel (Supplementary material, Figure 5). Overall, it can be seen that choice of kinetics can affect the simulation results. Given that DT kinetics is better representative of experimental data, the new methodology is helpful in using DT kinetics in a cell model.

Summary
An iterative procedure was developed to calculate average agglomerate reaction rate using double trap kinetics. An m th order Tafel kinetics was fitted to the DT kinetics. The fitted Tafel parameters and reaction order are used to obtain general Thiele modulus and effectiveness factor, which are used to compute the average agglomerate reaction rate. The reaction rate from the analytical model is within 2% of the numerical results over a wide range of physical and operating parameters. The analytical model enables the use of DT in a full-cell model without being computationally expensive. Cell simulations with DT and Tafel kinetics show difference between the two models, demonstrating the need for an analytical agglomerate model with DT kinetics.