Introduction

Polymer electrolyte membrane fuel cells (PEMFCs) are environmentally friendly energy sources that can alleviate environmental problems because of their low emissions of environmental pollutant gases such as SOx, NOx, CO2, and CO1,2. PEMFCs have been used in various applications such as fuel-cell vehicles and power supplies (including portable power supplies) because, in addition to their ecofriendly benefits, PEMFCs can also generate high power densities and operate with short start-up times because of their low operating temperature1,3,4,5. In general, PEMFCs consist of membrane electrode assembly (MEA) layers, gas-diffusion layers, microporous layers, gas flow channels, and bipolar plates3. The MEA layers are especially important because the cell performance and durability of PEMFCs are strongly affected by the design and composition of their MEA layers4. MEA layers in a PEMFC system comprise catalyst layers (CLs) and a polymer membrane; the polymer membrane plays a critical role in transferring protons from the anode to the cathode in the process of generating electricity5. CLs also play a critical role because the electrochemical reactions related to energy conversion in a PEMFC, such as the hydrogen oxidation reaction or oxygen reduction reaction (ORR), occur in CLs. The structure of CLs includes a carbon matrix with a Pt catalyst (Pt/C) and proton-conducting ionomers. Importantly, the proton transfer performance of a CL is affected by its ionomer thin film on Pt/C because protons can directly reach the Pt surface through the hydrated ionomer thin film. Therefore, the composition and morphology of the ionomers strongly influence the performance of a PEMFC.

The proton-conducting ionomers in PEMFCs are categorized as a perfluorinated acid (PFSA), nonfluorinated hydrocarbon, or an acid–base complex1. Among these various polymer ionomers, PFSA ionomers such as Nafion (DuPont)6,7, Aquivion (Dow Chemical)8,9, Aciplex-S (Asahi Glass)10, and Flemion (Asahi Glass)10,11 have been extensively used in PEMFCs because of their excellent proton conductivity and good mechanical, chemical, and thermal stabilities12. PFSA-based ionomers are composed of main chains of polytetrafluoroethylene (PTFE) and side chains terminated by sulfonic acid groups. Both experimental13,14,15 and theoretical16,17,18,19,20,21,22,23,24 investigations have been performed to elucidate the proton transfer and O2 permeation mechanisms in PEMFC systems because the length of the side chains in PFSA ionomers affects PEMFC performance.

Garsany et al.13 and Siracusano et al.14 experimentally investigated PFSA ionomers—specifically, Aquivion, which has a short side chain, and Nafion, which has a long side chain—to improve PEMFC performance. They concluded that Aquivion exhibits better cell performance than Nafion in PEMFCs because the Aquivion ionomer in the cathode CLs has a lower proton transport resistance, lower charge transfer resistance for the ORR, and lower mass transport resistance than the Nafion ionomer. Baschetti et al.15 investigated gas permeation in Nafion and Aquivion ionomers at various temperatures and relative humidities. Humidity and temperature have especially strong effects on gas permeability, and Nafion 117 ionomer was found to exhibit greater O2 gas permeability than Aquivion at 50 ℃.

Several groups have investigated the relationship between the diffusion coefficients of water and hydronium ions and PFSA morphologies in systems with various water contents and at different temperatures using molecular dynamics (MD) simulations16,17,18,19,20,21. The diffusion coefficients of water and hydronium ions were found to increase with increasing PFSA water content and increasing temperature. In addition, the sulfur–sulfur interatomic distance in PFSA increased with increasing PFSA water content. MD simulations22,23,24 have also been performed to investigate O2 permeation in Nafion ionomers in CLs. Kurihara et al.22,23 investigated the permeation of O2 gas into a Nafion ionomer on a Pt surface using MD simulations because such simulations are useful for understanding the nanoscale structures in the CLs of PEMFCs. They concluded that the diffusion coefficient of O2 molecules increased and the solubility of O2 molecules decreased with increasing water content in the Nafion ionomer. Jinnouchi et al.24 also used MD simulations to investigate O2 permeation through a Nafion thin film on a Pt surface, where the water content of the Nafion film was varied. Their results indicated that O2 permeation in Nafion increased with increasing water content and that understanding the behavior of O2 in PFSA ionomers on a Pt surface is critical to understanding its permeation properties. Therefore, the aforementioned experimental results indicate that the length of the side chain in PFSA ionomers can affect both the performance of PEMFCs and the O2 permeation behavior. Therefore, studies comparing of the O2 permeation properties of Nafion and Aquivion are needed to elucidate the effect of side-chain length in PFSA ionomers in PEMFCs.

In the present study, computational simulations using the full atomistic MD simulation technique are carried out to obtain detailed molecular information for calculating the transport properties of hydrated PFSA ionomers with various water contents on a Pt surface. In addition, the O2 permeation properties of hydrated PFSA with different side-chain lengths were measured at the interfacial region on the Pt surface. Therefore, two types of PFSA ionomers—Nafion (longer side chain) and Aquivion (shorter side chain)—were prepared for measurement of the O2 permeability at various water contents, enabling the relationship between the hydrated PFSA structure and the O2 permeation properties to be elucidated. In addition, the distribution of O2 and water in PFSA ionomers on a Pt surface were also analyzed using density profiles and radial distribution functions (RDFs) with various water contents at the operating temperature of a PEMFC (353.15 K). We expect that the results of this study will provide detailed information about O2 permeability of water-containing PFSA ionomers on a Pt surface and can provide guidance for the design of PFSA ionomers for use in PEMFCs.

Computational details

Model preparation

Figure 1 shows the chemical structures of the Nafion and Aquivion ionomers. Each ionomer was composed such that each polymer had 10 repeat units with 10 sulfonic acid groups per polymer chain. The molecular weight of the Nafion and Aquivion polymers was 9969.83 g/mol and 8309.63 g/mol per polymer chain; equivalent weights (EWs) of ~ 1000 g/mol and ~ 830 g/mol were applied, respectively. Water, O2 molecules, and hydronium ions were prepared for constructing hydrated PFSA ionomers. The components of each PEMFC system are summarized in Table 1. Materials Studio visualization software was used for the atomic visualization25.

Figure 1
figure 1

Chemical structures of PFSA component at (a) ionized Nafion (EW = 1000), (b) ionized Aquivion (EW = 830), and (c) water, (d) hydronium ion, and (e) O2 molecules.

Table 1 Details of parameters for the MD simulations.

Force-field and MD parameters

To describe inter- and intramolecular interactions in the Nafion and Aquivion in PEMFC systems, we applied a modified DREIDING force field26 in our simulations. The DREIDING force field has been widely used to describe PEMFCs systems27,28,29,30,31,32,33. The force fields of water molecules and Pt atoms were applied using F3C force field34 and the embedded-atom method (EAM) force field, respectively35. The total potential energy Etotal in PEMFC systems can be calculated using to Eq. (1):

$$E_{{{\text{total}}}} = E_{{{\text{vdW}}}} + E_{{\text{Q}}} + E_{{{\text{bond}}}} + E_{{{\text{angle}}}} + E_{{{\text{torsion}}}} + E_{{{\text{inversion}}}} + E_{{{\text{EAM}}}}$$
(1)

where \({E}_{\mathrm{vdW}}\), \({E}_{\mathrm{Q}}\), \({E}_{\mathrm{bond}}\), \({E}_{\mathrm{angle}}\), \({E}_{\mathrm{torsion}}\), \({E}_{\mathrm{inversion}}\), and \({E}_{\mathrm{EAM}}\) are the van der Waals, electrostatic, bond-stretching, angle-bending, torsion, inversion, and the EAM energies, respectively. For calculating entire MD simulations for PEMFC systems, the large-scale atomic/molecular massively parallel simulator (LAMMPS) code36,37 from Plimpton at Sandia was used. All MD simulations were carried out using the velocity Verlet algorithm38 to integrate equations of atomic motion, with a time steps of 1 fs. The electrostatic interactions in our systems were calculated using the particle–particle particle–mesh method39. The charges of particles in Nafion and Aquivion were calculated via density functional theory (DFT) calculations using the Mulliken charge analysis method40 in the Materials Studio software25. All DFT calculations for charge analyses were carried out using the double numerical basis set with polarization (DNP) function and the generalized gradient approximation with the Perdew–Burke–Ernzerhof functional41.

Force-field parameters between PEMFC components and the Pt surface

We used the nonbonded interaction energies reported by Brunello et al.42 to describe the interactions of Pt atoms with Nafion, Aquivion, water, and hydronium ions. In addition, for the interaction energies between O2 and a Pt slab, we calculated van der Waals parameters via DFT calculations to describe detailed intermolecular interactions using a Pt (111) slab with three atomic layers with periodic boundary conditions (PBCs) of 8.324 × 8.324 × 25.000 Å3, as shown in Fig. 2. The DFT calculation details were the same as those used in the charge analysis (Sect. 2.2.1), and a semi-empirical dispersion correction (DFT-D) with the Tkatchenko–Scheffler scheme was additionally applied43. Band-structure calculations with k-points were performed with a 4 × 4 × 1 Monkhorst–Pack k-point mesh44.

Figure 2
figure 2

(a) The optimized structures of O2 molecules on the Pt (111) surface. (b) The change in potential energies according to changes in distance between O2 and the Pt (111) surface, as determined using DFT calculations, fitting data, and DREIDING force field.

Model construction

To generate the Pt (111) surface and hydrated PFSA ionomers, we constructed a Pt (111) slab with five atomic layers with PBCs of 38.446 × 38.845 × 300.000 Å3. A z-direction length of 300 Å was used to prevent interaction beyond the PBCs. To construct hydrated PFSA ionomers on a Pt (111) surface, the λ (water molecules per sulfonic acid group) values were 2.92, 6.15, 9.77, and 13.83 for hydrated Nafion and Aquivion ionomers. The PFSA ionomers were composed of six polymer chains with 60 hydronium ions to maintain electrical neutrality. The Monte Carlo simulation in the Amorphous Cell module of the Materials Studio software25 was used to construct the initial configuration of hydrated PFSA ionomers on a Pt (111) surface.

MD simulations

After the initial hydrated PFSA with Pt (111) surface models were constructed, the temperature was gradually increased to the operating temperature (353.15 K) of PEMFCs with canonical ensemble (NVT) for 1 ns. Consequently, 15 ns MD simulations were performed by NVT simulation to obtain the equilibrated structures. To simulate the permeation of O2 into PFSA ionomers, 735 O2 molecules with a thickness of 100 Å were placed on the equilibrated PFSA ionomers. After the O2 molecules were added, a total of 150 ns of NVT simulations at 353.15 K was performed. The last 10 ns of NVT simulation was used for data collection of the O2 permeation properties.

Results and discussion

van der Waals parameters for O2 and Pt

For a better description of the O2 permeation process inside a hydrated Nafion ionomer thin film on a Pt (111) surface, the atomic interaction curve between O2 and the Pt surface was reproduced by DFT under the framework of the DREIDING force field. Figure 2a shows O2 adsorbed onto the Pt surface, which was built for calculation of the adsorption energy as a function of the z-distance. The calculated adsorption energy as a function of distance was fitted to the Lennard–Jones potential in Fig. 2b, which well reproduced the results of DFT calculations. The Lennard–Jones potential function is shown in Eq. (2):

$$E=\varepsilon \left[{(\frac{{r}_{\mathrm{m}}}{r})}^{12}-2{(\frac{{r}_{\mathrm{m}}}{r})}^{6}\right]$$
(2)

where \(E\) indicates the potential energy with changing distance \(r\) and \(\varepsilon\) is the depth of the potential well at distance \({r}_{\mathrm{m}}\). The values of \(\varepsilon\) and \({r}_{\mathrm{m}}\) for the oxygen atoms in O2 molecules on a Pt surface are 4.070 kcal/mol and 2.338 Å, respectively. The fitted interaction well describes detailed interactions that the DREIDING force field cannot describe. Thus, the fitted interaction between O2 and Pt surfaces was used in the MD simulation to analyze the process of O2 permeation into a Nafion ionomer coated onto a Pt surface.

Equilibrated structure

Equilibrated structures in Fig. 3a–h were obtained from the data corresponding to the final 5 ns of the MD trajectories. Figure 3a–d show a hydrated Nafion thin film on a Pt surface, and Fig. 3e–h show a hydrated Aquivion thin film on a Pt surface, both under different hydration levels, λ. As λ increases, the Nafion and Aquivion films become gradually segregated into hydrophilic water clusters and hydrophobic regions with PTFE backbones of Nafion and Aquivion ionomer, respectively. Water molecules are predominantly adsorbed as a thin layer onto the Pt surface because they strongly interact with this surface30.

Figure 3
figure 3

Snapshots of equilibrated systems comprising the Pt (111) surface, hydrated Nafion ionomers, and O2 molecules with λ values of (a) 2.92, (b) 6.15, (c) 9.77, and (d) 13.83; snapshots of equilibrated systems comprising the Pt (111) surface, hydrated Aquivion ionomers, and O2 molecules with λ values of (e) 2.92, (f) 6.15, (g) 9.77, and (h) 13.83.

To quantitatively analyze the permeation of O2 molecules, we investigated the density profile of hydrated Nafion, water molecules with hydronium ions, and O2 molecules, as shown in Fig. 4a–d. Dotted lines represent the average density of a bulk Nafion membrane, as reported in our previous study20, which indicates that the density of hydrated Nafion on a Pt surface is not substantially different from that of a bulk membrane except for the density at the Nafion–Pt interface. As shown in the equilibrated structures in Fig. 3a–h, hydrated Nafion ionomers and water molecules with hydronium ions exhibit the highest density at the hydrated Nafion–Pt interface because of their strong attractive interaction. The purple line indicates the point of O2 solvation from the hydrated Nafion–gas interface, which is determined by the distance corresponding to the average number of total solvated O2 molecules. At this point, the density of O2 molecules abruptly decreases from its maximum value. These results are in good agreement with those of Jinnouchi et al.24, who reported that an energy barrier at the Nafion–gas interface dominates the solubility of O2 in hydrated Nafion. Inside hydrated Nafion, dissolved O2 molecules exhibit the highest density and water molecules with hydronium ions exhibit the lowest density. This trend becomes more discernible as the hydration level increases. These results suggest that O2 molecules are not preferentially positioned inside hydrophilic domains but rather at the interfacial regions between hydrophobic and hydrophilic regions.

Figure 4
figure 4

Density distributions of hydrated Nafion ionomer (with included water and hydronium ions; black lines), water and hydronium ions (blue lines), and O2 molecules (red lines) on the Pt (111) surface along the thickness direction, with λ values of (a) 2.92, (b) 6.15, (c) 9.77, and (d) 13.83. Dashed lines show the previously reported density for hydrated bulk Nafion20. The purple line indicates the point of O2 solvation occurs from the interface of gas/hydrated Nafion.

We also investigated the density profile of water molecules, hydronium ions, and O2 molecules in hydrated Aquivion; the results are shown in Fig. 5a–d. Like Nafion, hydrated Aquivion and water molecules with hydronium ions show the highest density at the Aquivion–Pt interface because of their strong attractive interaction. The hydrated Aquivion is thinner than the hydrated Nafion because of Aquivion’s shorter side chains and lower EW. The density of O2 molecules also abruptly decreases the maximum value at the distance indicated by the purple line, which represents the distance at which O2 solvation begins. Inside hydrated Aquivion, the dissolved O2 molecules exhibit the highest density and water molecules with hydronium ions exhibit the lowest density.

Figure 5
figure 5

Density distributions of hydrated Aquivion ionomer (with included water and hydronium ions; black lines), water and hydronium ions (blue lines), and O2 molecules (red lines) on the Pt (111) surface along the thickness direction, with λ values of (a) 2.92, (b) 6.15, (c) 9.77, and (d) 13.83. Dashed lines show the previously reported density for hydrated bulk Aquivion20. The purple line indicates the point of O2 solvation occurs from the interface of gas/hydrated Nafion.

O2 solubility and permeation

To quantify the solvation of O2 molecules by hydration level, we calculated the solubility on the basis of the average number of dissolved O2 molecules at the distance indicated by the purple vertical line in Figs. 4 and 5. The solubility of Nafion and Aquivion ionomers as a function of the hydration level is presented in Fig. 6a. The solubility of O2 decreases with increasing hydration level in both the Nafion and Aquivion ionomers. Because the interfacial region between hydrophilic and hydrophobic domains becomes relatively limited when the net ionomer volume increases under hydration, the number of sites available for dissolved O2 molecules also becomes limited. Therefore, the ionomers exhibit low O2 solubility at high hydration levels (λ). At the same λ, the O2 solubility is higher in the Nafion ionomer than in the Aquivion ionomer, which means that the solvation of O2 molecules in the Aquivion ionomer is restricted compared with that in the Nafion ionomer.

Figure 6
figure 6

O2 (a) solubility and (b) permeability in hydrated Nafion and Aquivion on a Pt (111) surface.

Using the solubility values, we derived the permeability coefficient of O2 molecules. Gases permeate through ionomers via a solution-diffusion mechanism, where the dissolved gas molecules diffuse into the ionomers45. Thus, the permeability coefficient P is described by the equation

$$P=D\times S$$
(3)

where D is the diffusion coefficient and S is solubility. The diffusion coefficient in our system was obtained from the self-diffusion coefficient of O2 molecules in the bulk structures of Nafion and Aquivion. The calculated permeability coefficients for O2 molecules in Nafion and Aquivion ionomer are shown in Fig. 6b; these values are consistent with those reported by Baschetti et al.15. With increasing hydration level λ, the permeability of O2 in both Nafion and Aquivion ionomers increases until λ = 9.77. At λ > 9.77, the permeability increases slightly because, despite the diminished solubility, the self-diffusion coefficient of O2 in a bulk ionomer membrane increases with increasing hydration level; thus, permeability increases with increasing hydration level. In our previous work20, we reported that O2 molecules exhibit greater diffusion in Aquivion ionomers than in Nafion ionomers because of the Aquivion ionomers’ better-developed water channels. However, at the same hydration level λ, both the permeability and solubility of O2 are higher in Nafion ionomer than in Aquivion ionomer. This result means that solubility is more critical to permeability, which is dominated by the availability of solvation sites for dissolved O2 molecules in the ionomer.

RDF analysis

The interface region between hydrophilic and hydrophobic domains in a hydrated PFSA ionomer, where O2 molecules preferentially dissolve, is most likely to be the side-chain part above the sulfonic acid groups. That is, O2 solvation mainly occurs at the side chains of the ionomers. Thus, the difference in solubility between Nafion and Aquivion is reasonably deduced to arise from their different side-chain structures. In this regard, we analyzed the structure between O2 and the main component of the Nafion and Aquivion side-chain structures to understand the difference in O2 solubility between them. As shown in Fig. 7a, we analyzed correlations between Carbon(center)–Oxygen(O2) and Oxygen(side chain)–Oxygen(O2) using RDFs. The RDF of each pair is described by the following equation:

Figure 7
figure 7

(a) Molecular model with the atomic definition for O2 location analysis using RDFs. RDFs at (b) C(center)–O(O2) and c) O(side chain)–O(O2) and the first coordination number of (d) C(center)–O(O2) and (e) O(side chain)–O(O2).

$$g_{A - B} \left( r \right) = \left( {\frac{{n_{B} }}{{4\pi r^{2} dr}}} \right)/\left( {\frac{{N_{B} }}{V}} \right)$$
(4)

where \({n}_{B}\) is the number of B particles located at a distance r in a shell of thickness \(dr\) from particle A, NB is the number of B particles in the system, and V is the total volume of the system; NB/V can be represented by the number density, ρ. The ρg(r) of each pair of Nafion and Aquivion ionomers is presented in Fig. 7b, c, respectively. The ρg(r) of both C(center)–O(O2) and O(side chain)–O(O2) shows the first peak at 7.87 Å and 8.51 Å, respectively. These values decrease with increasing hydration level. At the same hydration level, the ρg(r) of both C(center)–O(O2) and O(side chain)–O(O2) is greater in the Nafion ionomer than in the Aquivion ionomer. This result means that O2 has a more favorable correlation with C(center) and O(side chain) in the Nafion ionomer than with those in the Aquivion ionomer. The first coordination number (CN) at the distance of the first peak was calculated (Table 2). As shown in Fig. 7d, e, the CN of C(center)–O(O2) and O(side chain)–O(O2) pairs also decreases as the hydration level increases with increasing ρg(r). At the same hydration level, the Nafion ionomer shows a higher CN of each pair than the Aquivion ionomer. This result suggests that more O2 molecules are coordinated to the side chain of the Nafion ionomer than to that of the Aquivion ionomer, especially to the side chains’ oxygen atoms, which results in greater solubility of O2 in the Nafion ionomer. In addition, we note that Aquivion shows a similar CN as Nafion at λ = 2.92. As described in the previous section, we also observed that the difference in solubility and permeability is more discernible at a hydration level greater than λ = 2.92. At a low hydration level, the site for O2 solvation is not significantly different due to only a few water molecules located near the side chain leading to a lack of phase segregation, and O2 molecules can be easier located near the side chain of PFSA than other water contents. However, as the hydration level increases, the solvation sites of O2 molecules become limited to the side-chain region of the ionomer due to the water molecules formed the water cluster near sulfonic acid groups of PFSA. Although the water molecules formed the water cluster near sulfonic acid groups, O2 molecules can be located near the side chain of Nafion easier than Aquivion because the side chain of Nafion is longer than that of Aquivion, leading to higher CNs of Nafion. The CN of C(center)–O(O2) and O(side chain)–O(O2) pairs of Aquivion were sharply decreased at λ = 6.15 because we assume that the water cluster can be formed in Aquivion easier than Nafion due to higher sulfonic acid group concentration. After that, the CNs of Aquivion at 9.77 to 13.83 was gradually decreased with developing water clusters and channels near the side chain. Consequently, Nafion and Aquivion exhibit an observable difference in solubility at higher hydration levels.

Table 2 The coordination numbers of Nafion and Aquivion on a Pt surface at 353.15 K, where the water contents (λ) were varied.

Conclusion

We prepared two types of hydrated PFSA ionomers on a Pt (111) surface and varied their water contents (2.92 ≤ λ ≤ 13.83) to analyze the distribution morphologies and permeability properties of O2 molecules in detail using full atomistic MD simulations at the operating temperature of PEMFCs (353.15 K). The hydrated PFSA ionomer structures were prepared using Nafion (long side chain) and Aquivion (short side chain). In the MD simulations, O2 molecules were gradually permeated into hydrated PFSA ionomers on the Pt surface. Density profile analysis for hydrated PFSA ionomers was performed to calculate the O2 solubility in hydrated PFSA ionomers with various λ. The O2 solubility in hydrated Nafion ionomers on the Pt surface was higher than that in hydrated Aquivion ionomers on the Pt surface at the same λ values. The RDF analysis showed that the first CN of C(center)–O(O2) and O(side chain)–O(O2) pairs of hydrated Nafion ionomers was greater than that of hydrated Aquivion ionomers at the same water content. These results suggest that more O2 molecules are coordinated to the side chains of the Nafion ionomer than to those of the Aquivion ionomer, especially oxygen atoms, which results in greater O2 solubility in the Nafion ionomer.