A Multiparadigm Modeling Investigation of Membrane Chemical Degradation in PEM Fuel Cells

The details of the computational approach based on CGMD simulations are explained elsewhere^43–45,50 and is developed in two major steps. In the first step, Nafion chains, water and hydronium molecules are replaced by corresponding spherical beads with predefined sub-nanoscopic length scale. In the second step, parameters of renormalized interaction energies between the distinct beads are specified.

We consider four main types of spherical beads: polar, nonpolar, apolar, and charged beads.⁵⁴ Clusters including a total of four water molecules or three water molecules plus a hydronium ion are represented by polar beads of radius 0.43 nm. The simulation box contained 72 coarse-grained Nafion chains each consisting of 20 monomers and 20 side chains. Each monomeric unit is represented by two apolar beads for backbone (red) and one single polar bead (green) for the entire sidechain (including etheric group) as depicted in Figure 1. We adopted the same coarse graining strategy as in our previous work⁴⁴ which was also suggested earlier in Ref. 55. A side chain unit has a molecular volume of 0.306 nm, equivalent to the molecular volume of a four-monomeric unit of polytetrafluoroethylene (PTFE), of size 0.325 nm. Thus in our coarse-graining, a monomeric backbone unit, i.e.,–CF₂–CF₂–CF₂–CF₂–CF₂–CF₂–CF₂–CF₂–, is represented by two beads and a perfluorinated ether sulfonic sidechain is represented by one single bead. Our coarse-graining strategy requires all beads have identical volume set at 0.315 nm³. The selected box size does not impact the cell performance since the pore structure remains unchanged for box volumes larger than 50 × 50 × 50 nm³ and we do not expect significant changes in the O₂ crossover or in the proton conductivity.

The interactions between non-bonded beads are modeled by the Lennard-Jones (LJ) potential

where the effective bead diameter (r_ij) is 0.43 nm for side-chain, backbone and water beads. The strength of interaction (D⁰) is limited to five possible values ranging from weak (1.8 kJ/mol) to strong (5 kJ/mol) beads.⁵⁴ The electrostatic interactions between charged beads are described by the Coulombic interaction

with relative dielectric constant ɛ_r = 20 in order to include screening. The effect of solvent is incorporated by changing ɛ_r as well as by varying the degree of dissociation of Nafion side chains. Interactions between chemically bonded beads (in Nafion chains, for example) are modeled by harmonic potentials for the bond length and bond angle

where the force constants are K_bond = 1250 kJ/mol.nm and K_angle = 25 kJ/mol.radian² respectively. r₀ and θ₀ are the equilibrium bond length and angle. The size of the simulation box can vary from 50 × 50 × 50 nm³ to 500 × 500 × 500 nm³, depending on the system and the composition. We conducted an annealing procedure over a period of 50 ps by increasing the temperature from 298 to 398 K, followed by a short MD simulation for 50 ps in a NVT ensemble, followed by a cooling procedure down to 298 K.⁴⁴ We did not observe any drift to a more ordered state after the equilibration procedure. For analysis purposes, however, final trajectories after full equilibration were used.

The degradation process is simulated as the following: first, the initial morphologies are generated using similar set of parameters and process that are explained elsewhere.⁴⁴ In order to mimic the degradation process, Nafion sidechains are randomly detached at different percentage from the backbone under various water contents. In our simulations, the detached sidechains are assumed to remain in the water phase inside the Nafion pore as "dissolved" sulfonic anions, thus the system remains electro-neutral.

We have used several tools to calculate structural properties of the degraded membrane to investigate the impact of side chain losses on the microstructural properties of the membrane. These techniques include Radial Distribution Function (RDF), Pore Side Distribution, Cluster Size analysis, and Pore Network Analysis.^43,44,50

General framework

The performance model is a multi-paradigm multiscale single-cell model implemented within our in house simulation package MS LIBER-T (see Figure 3).^56–58 This is a software coded in a modular framework on an independent C/Python language basis, highly flexible and portable with multiple application domains already demonstrated.^59,60 Similar to the previous models developed by Franco et al., the single-cell PEMFC model in MS LIBER-T represents explicitly the physical mechanisms at different scales as nonlinear sub-models in interaction (modularity)⁶¹ and it is designed to calculate electrochemical signals (e.g. polarization curves, cell potential vs. time, etc.) from the chemical and structural properties of the materials.^62–65

Zoom In Zoom Out Reset image size

For simplicity reasons, in this paper we restrict ourselves on the cell model based on the following overall assumptions:

• isothermal conditions;
• anode operating with water saturated H₂ and cathode operating with water saturated O₂. Under this assumption, the PEM hydration maintains at λ = 19. From this, only the CGMD database corresponding to this water content is used for the performance decay calculations shown in this paper. We underline that the model is general and other water contents and/or water content transients may be also studied, provided that H₂O transport is resolved across the PEM.Indeed, the mechanisms described are
• H⁺ transport across the membrane electrodes assembly within a 1-D approach;
• e⁻ transport across the CLs and gas diffusion layers within a 1-D approach;
• the presence of H₂ and O₂ in the on-catalyst ionomer film inside the CLs within the approach reported in Ref. 67. In the case of the anode, a 0D mass balance equation between the O₂ crossover rate from the PEM to the anode and the O₂ consumption rate on the catalyst, is solved to calculate the evolution of the resulting O₂ concentration in the on-catalyst ionomer film (Figure 4b);
• the interfacial nanoscale electrochemical mechanisms at the vicinity of the catalyst including both elementary kinetics and electrochemical double layer effects within an on-the-fly coupled KMC/continuum approach.^66,67 The impact of the surface roughness on the electrochemical double layer structure, as addressed in some publications,^68,69 is not taken into account in the present model.

We use the version of our model which treats the elementary kinetic reactions in the CLs through the Kinetic Monte Carlo Electrochemical Variable Size Method (KMC E-VSSM) that we have introduced and discussed in Ref. 67 for the cathode CL. This approach resolves the adsorption/desorption, adspecies surface diffusion and reactions on the catalyst surface during the PEMFC operation, and allows calculating the electrodes and cell potential.

In general, the performance model can account for the complete cycle of the PEM chemical aging, i.e.: H₂O₂ formation in the anode from the O₂ crossover from the cathode, H₂O₂ formation in the cathode through the ORR, the diffusion of H₂O₂ in the PEM and H₂O₂ decomposition into OH° and OOH° radicals in presence of Fe²⁺/Fe³⁺ Fenton's cations, the degradation removal of the side chains and the impact on the PEM proton conductivity and cell performance decay (Figure 4).

For demonstration purposes in this paper, the degradation kinetics is assumed to be of first order toward the OH° concentration, i.e.:

where R and R^' refers to the fresh and aged polymer respectively, SC refers to a single side chain and D to the detached side chain. The effect of hydroxide radical in the CGMD simulations is captured by presence of water and hydronium ion and therefore is not explicitly added to the CGMD simulations.

The model implements the CGMD database with spatially resolved 1D resolution across the PEM thickness. The PEM model, comprising the mathematical descriptions of the relevant transport and chemical reaction mechanisms described in the following sections, is coded in Python and coupled as an additional module in MS LIBER-T.

O₂ transport model in the PEM

The transport of O₂ is assumed to be simply governed by Fickean diffusion and it is spatially resolved in 1D through the PEM thickness (the H₂ transport was neglected due to its high reactivity at the high potentials typically found in the cathode):

H₂O₂ production models

The HOR in the anode is modeled in competition with a H₂O₂ production reaction through our KMC E-VSSM approach with DFT-calculated kinetic activation energies. The elementary reaction steps considered are detailed in Table I.⁶² Parameter values of the associated electrochemical double layer sub-model and of the surface diffusion processes are identical to the ones used in our previous work.⁶⁷

Table I. Elementary kinetic model for the H₂O₂ production in the anode CL.

	Activation energies (kJ.mol−1)
Reaction	Eactf	Eactb
	68	∞
	34	34
	53	53
	—	∞
	41	39
	35	68
	19	∞

For the H₂O₂ production within the ORR pathway in the cathode CL, we have considered the elementary kinetic reaction in Table II, still within our KMC E-VSSM approach, with the same parameters values for the activation barriers, diffusion barriers and electrochemical double layer sub-model that in Ref. 67. The reaction steps and the activation barriers values were determined by DFT calculations carried out on Pt(111) surfaces.^62,70 The modular character of the model allows considering here other DFT databases, with different H₂O₂ production rates, corresponding to Pt bimetallic catalysts. This will be the subject of a future publication.

Table II. Elementary kinetic model for the H₂O₂ production in the cathode CL.

	Activation energies (kJ.mol−1)
Reaction	Eactf	Eactb
	—	∞
	30	149
	38	43
	1	158
	88	94
	19	80
	24	47
	—	∞

For the adsorption steps we follow the same approach as in our previous publication:⁶⁷ according to the collision theory, the kinetic parameters are given by

where P and m are the partial pressure at x = L (that can be related to the concentration of the dissolved reactant in the electrolyte) and the atomic mass of a reactant respectively. sc is the sticking coefficient estimated from published values.^71,72 For the reaction steps, according to the extended transition state theory, the kinetic parameters are given by

where k_ij and E_ij,act are the kinetic rate constant and activation energy along the minimal energy path for going from state i to state j. The activation energies E_ij,act are written as

In Equation 16, f(σ) is the effect of the interfacial electric field, function of the catalyst charge density σ, onto the effective kinetics.⁶² This quantity together with the concentration of protons in the electrolyte at the vicinity of the reaction sites, is calculated by means of our non-equilibrium electrochemical double layer model in our work.⁶⁶

The KMC E-VSSM allows then the calculation of the evolution of the adspecies coverage on the anode and cathode catalyst surfaces.

H₂O₂ and iron ions reaction-transport models

The transport of H₂O₂ and iron ions is assumed to be diffusive and they are spatially resolved in 1D through the PEM thickness from the solution of

where the source/sink term is related to the Fenton reactions describing the H₂O₂ decomposition in OH°/OOH° in presence of the Fenton's cations. In this paper, for demonstration purposes, only the presence of Fe²⁺ and the Fenton reactions reported in Table III were considered.

Table III. Elementary kinetic model for the Fenton's reactions. The kinetic parameter values are from Ref. 34.

		Kinetic parameter (s−1)
Reaction	Kinetic rate	ki
		63.10−3
		3.3.105

The concentration of the reaction intermediates, reactants and products are calculated from the mass balances:

Membrane chemical aging model

The side chain concentration is resolved from the following balance equation:

The bulk electrolyte potential at the cathode is given by

where is the effective proton conductivity calculated from the CGMD database following Equation 6.

Some of the parameters values used are reported in Table IV. The values of the parameters not reported here are the same as in our previous publications.⁶²

Table IV. Parameters and numerical values.

Parameter	Numerical value (units)	Source
Temperature (T)	353 K	assumed
kDEG	100	Ref. 35
Selectrode	2.5 × 10−5 m2	Ref. 8
Fe2 +flux	2.7 × 10−10 kmol/sec.m3	assumed
γ (Figure 4b)	10−5	assumed
	1 × 10−9 m/sec	assumed
	4 × 10−8 m/sec	assumed
Membrane thickness	2.5 × 10−5 m	assumed
(BOL)	2 × 10−4 S/m	assumed
λ	15	assumed

PEM simulations were performed with λ = 2, 4, 9, 15. The mesoscopic structure of the hydrated membrane at  λ = 9 is visualized in Figure 5, revealing a sponge-like structure, similar to structures obtained by other meso-scale simulations.^75,73 Water beads together with hydrophilic beads of sidechains form clusters, which are embedded in the hydrophobic phase of the backbones. The detailed structural analysis indicates that the hydrophilic subphase is composed of a three dimensional network of irregular channels. The typical channel sizes are from 1 nm, 2 nm, and 4 nm at λ = 2, 4, 9, 15. This corresponds roughly to linear microscopic swelling. The site-site radial distribution function (RDF), obtained from CGMD simulations, matches very well to those from the atomistic MD simulations.⁴³ The RDF between the sidechain beads and the other components of the mixture show that side chains are surrounded with water and hydrated protons. The autocorrelation functions exhibit similar dependences on bead separation at all λ, indicating a strong clustering of sidechains due to the aggregation and folding of polymer backbones.⁷⁴ The degree of ordering of water near polymer|water interfaces decreases with increasingλ.

Zoom In Zoom Out Reset image size

So far, coarse-grained approaches offer the most viable route to the molecular modeling of self-organization phenomena in hydrated ionomer membranes.

Admittedly, the coarse-grained treatment implies simplifications in structural representation and in interactions, which can be systematically improved with advanced force-matching procedures. Moreover, it allows simulating systems with sufficient size and sufficient statistical sampling. Structural correlations, thermodynamic properties, and transport parameters, can be studied.

Figure 6 illustrates tortuosity factor as a function of the sidechain loss (%). Tortuosity factor is defined as a ratio of the geometrical pore length and pore axis. The insert depicts average pore size in nm as function of the degradation at different water contentment (λ = 2, 4, 9, 15). The diameters of water channels vary in the range of 1–7 nm, exhibiting a roughly linear increase from low to high water content. The average separation of side chains increases as well with water content, which indicates a continuous structural reorganization of polymer aggregates upon water uptake.⁴⁴ This could involve backbones sliding along each other in order to adopt more stretched conformations. The sidechain separation varies in a range of 1nm or slightly above. The network of aqueous domains exhibits a percolation threshold at a volume fraction of ∼10%, which is in line with the value determined from conductivity studies. This value is similar to the theoretical percolation threshold for bond percolation on an fcc lattice. It indicates a highly interconnected network of water nanochannels.

Zoom In Zoom Out Reset image size

In Figure 7a we can see the side chain separation (nm) as a function of the degree of side chain degradation (%) and water content (λ). Overall, by increasing the hydration level, side chains move apart, with their mean separation at Beginning Of Life (BOL) increasing from 0.92 nm at λ = 2 to 1.48 nm for λ = 15. Figure 7a also shows that by increasing degradation (percentage of side chains detached from backbone), side chain separation decreases. Side chain separations on a single ionomer chain are between 1.5–1.7 nm.⁷⁵ Charge distribution and the magnitude of electrostatic interactions between side chains determine how ionomer backbones are assembled into fibrils or lamellae, where, the net density of side chains increases. Loss of sidechains reduces electrostatic hindrances between side chains on the surface of fibrils. The latter facilitates easier backbone folding, leading to lower side chain separations. Thus, the effective packing density of side chains due to better polymer degradation indicates that side separations and aggregate sizes decreases with side chain loss. We speculate that decrease in sizes of backbone aggregates corresponds to a decrease in electrostatic interaction between remaining side chains attached to the backbone, which causes increases in packing density of backbones,

Zoom In Zoom Out Reset image size

In Figure 7b we illustrate normalized average of sidechain separation (d_ss) for λ = 9 with respect to that for the Beginning of Life (BOL) membrane model where no sidechain losses are imposed (%deg = 0). Using a random-walk model for proton diffusion based on the Einstein-Smoluchowski equation^39,76 one can estimate the relationship between proton conductivity, proton diffusivity and mean step distance,

where z is a constant dependent upon the dimensionality of random walk (6 for three-dimensional walk), l is the mean step distance, i.e. side chain separations (d_ss) and τ_D is the mean time between successive steps. Notice that the above relationship does not necessarily mean protons transfer via a "hopping" Grotthuss mechanism³⁹ and can similarly describe the vehicular mass diffusions. Overall, the normalized relationship between side chain separation and degradation of sidechains depicted in Figure 7b indicates that proton conductivity is declined by increasing sidechain losses in accordance to the above equation.

The CGMD-generated database has been used in the performance model following the algorithm presented in Overall methodology section.

We carried out four different simulation cases: OCV, 0.1, 0.3 and 0.5 A/cm² applied current densities. For all the cases, the systems were initialized without considering the PEM degradation along the first 0.1 seconds of simulated time, until the steady state is reached: then the degradation process is "activated" by switching on a constant input flux of Fe^{2 +} (10⁻¹⁴ kmol/m³.sec). As it is explained in Figure 4, the model allows capturing the HOR and ORR intermediates coverage evolution in both anode and cathode as well as other outputs (such as faradaic current, H₂O₂ production rate, ...) not shown here. The PEM simulation sub-model calculates the evolution of the porosity/tortuosity, the SO⁻₃ concentration, conductivity and the associated electrostatic potential evolution across the PEM along the simulation by help of Equation 20.

In Figure 8a, we report the calculated SO⁻₃ concentration evolution in time for the four cases investigated. The initial SO⁻₃ concentration was assumed to be 1.2 kmol/m³ for all the cases. In all the curves in Figure 8a the concentration decays in time with an overall rate which depends on the value of the applied current density. Indeed, the SO₃⁻ degradation rate increases as the applied current density decreases, which is in agreement with the experimental knowledge. This trend is because of the faradaic current present in the anode, and thus the HOR overall rate (in spite of the H₂O₂ formation), increases with the imposed total current density as the faradaic current evolves to balance it. Briefly, the less H₂O₂ is formed in the anode, the less SO⁻₃ is degradated.

Zoom In Zoom Out Reset image size

In Figure 8b we report the corresponding normalized evolutions for all the cases under investigation. The SO⁻₃ consumption generates a decrease on the tortuosity (τ) (see Figure 6) and an increase on the porosity (ε). As one can see, under OCV conditions, this factor varies from 1 to ∼3 after 24 hours of simulated time which is in agreement with the experimental known fact that OCV conditions are the most damageable for the PEM.

In Figure 8c, the PEM normalized conductivity evolution on time is shown, with a trend which correlates with the increase on the factor and the loss of SO⁻₃.

We underline that, according to our model, the increase over time of the factor during the PEM degradation, leads to a monotonous increase of the effective O₂ diffusion coefficient. The resulting increase of the O₂ crossover leads to the increase of the O₂ adsorption and consequent O formation in the anode catalyst which favors the ORR and the H₂O production. In Figure 9 we present, as an example, the evolution of the adspecies coverage on the anode catalyst under OCV conditions. During the first ∼25 hours of simulated operation, the HOR governs the anode operation with a high coverage of H (∼0.5 ML) and other species such as H₂ (∼0.2 ML) and H₂O₂ (∼0.05 ML). Later, once the O₂ crossover becomes sufficiently high, the anode becomes "cathodized" depleting the overall cell potential. This process is self-maintained because of the retro-feedback with the O₂ crossover: in other words, the more H₂O₂ is produced, the more sidechain loss is induced and finally the more O₂ crossover results, leading again to more H₂O₂ production and so on. Not surprisingly, in Figure 10, we can observe that the PEM degradation does not affect the calculated ORR adspecies coverage for the cathode under OCV conditions.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 11, we report the cathode ORR coverage evolution for the case of an applied current density of 0.5 A/cm². As it was previously predicted⁶⁷ the H₂O is the dominant adspecie (with an initial coverage value of ∼0.6 ML for H₂O, 0.3 ML for O₂ and 0.03 ML for O). The PEM degradation process at non-zero current leads to the cathode coverage evolution. Indeed, the cathode responds by increasing the water production as a consequence of the electrostatic potential decrease over time through the PEM (conductivity depletion).

Zoom In Zoom Out Reset image size

In Figure 12, two different simulated time ranges are presented for the HOR coverage evolution in the anode at 0.5 A/cm². As one can see around ∼33 hours, the anode starts to generate more H₂O, O₂H and OH, which leads to the cathode potential (the cell potential) decrease over time (Figure 13). This occurs because at ∼33 hours the porosity/tortuosity factor suddenly starts to increase (see Figure 8b), affecting the effective O₂ diffusion coefficient and favoring an undesired ORR in the anode side.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 13 we summarize the cell potential evolution for the four studied cases. As known from experimental knowledge, under the assumption that PEM aging is the only degradation mechanism involved, the cell performance degradation rate increases as the applied current density decreases. For the three cases where the applied current density is non-zero, the potential decay clearly splits into two regions: first a region until ∼25 hours where the cell potential decays roughly linearly; then a second region where the cell potential decays in a non-linear fashion. According to our model, at non-zero applied current densities, the first region is governed by the conductivity decrease due to the decrease in the side chain separation (see Figure 7 and Equations 20 and 21); and the cell potential decay in the second region is due to the drastic O₂ crossover increase. At OCV conditions, the overall potential decay is related to the O₂ crossover increase.