Effects of Polymer Morphology on Proton Solvation and Transport in

Article pubs.acs.org/JPCC

Effects of Polymer Morphology on Proton Solvation and Transport in Proton-Exchange Membranes Shulu Feng, John Savage, and Gregory A. Voth* Department of Chemistry, James Franck Institute, and Computation Institute, University of Chicago, 5735 S. Ellis Ave., Chicago, Illinois 60637, United States ABSTRACT: The effects of polymer morphology on proton solvation and transport in hydrated Nafion are investigated by using a novel reactive molecular dynamics approach. Three of the most significant morphological models of Nafion, the lamellar model, the cylinder model, and the cluster-channel model, are studied. The three models exhibit distinct proton transport (PT) patterns, which result in different proton diffusion rates. In both the lamellar and the cylinder models, the interaction between protons and the sulfonate groups is shown to be the key factor in determining PT behavior. For the cluster-channel model, the geometrical shape also plays an important role in influencing the PT behavior. The change in the excess proton solvation structure as a function of the distance between protons and sulfonate groups is also analyzed. It is found that the increase of the water cylinder radius or water layer height leads to the presence of more protons around the sulfonate groups, while, for the cluster model, an increase in the water sphere radius leads to the presence of fewer protons around the sulfonate groups. Furthermore, for the lamellar and cylinder models, the hydrated protons around the sulfonate groups consist of more Zundel-like (H5O2+) structures when the hydration levels decrease, which is also influenced by the different morphological structures of Nafion.

1. INTRODUCTION The proton-exchange membrane (PEM) is one of the primary components used in fuel cell technologies. Extensive experimental1−14 and computational studies15−57 have been conducted with a goal of helping to improve the performance of PEMs, including the operational temperature range, mechanical properties, water management, selectivity, gas permeation, and proton conductivity. For example, composite electrolytes were found to contribute to the improvement of proton-exchange membrane fuel cell (PEMFC) performance at high temperatures (130 °C).8 A highly ordered mesoporous Nafion membrane with remarkable water retention ability was synthesized via a micelle templating method.7 The selectivity of different composite membranes was investigated, with the silicalite/Nafion having the highest selectivity.6 The sulfonate Calprene H6120 hybrid membranes exhibit better proton conductivity than the Nafion 117 composites.33 The computational study of proton transport (PT) in PEMs has been complicated by the fact that both standard diffusive (vehicular) proton motion and Grotthuss proton shuttling (hopping) contribute to the PT behavior. To date, only very few studies21,54,55,58 have appeared in the literature in which both PT mechanisms have been simulated within a common, deterministic molecular dynamics (MD) framework at length and time scales relevant to the PT behavior in PEMs. Other MD studies have mostly utilized either “classical” hydronium cation models, in which Grotthuss shuttling is not included, computationally demanding ab initio MD (AIMD) studies on very small systems for short time scales,59,60 or ad hoc simulation methods with proton hopping somehow included into the MD algorithms as a stochastic event.29,47 © 2012 American Chemical Society

The morphology of the PEM determines the confinement and connectivity of the aqueous phase, thereby controlling the macroscopic proton transport.61,62 Many attempts have been made to characterize the effects of the chemical configuration of underlying polymers, such as backbone length, side-chain length, and equivalent weight, on the overall morphological characteristics of the PEM, such as phase separation and cluster size distribution.60,63−77 The side-chain length is particularly important in determining the overall PEM morphology. Molecular dynamics (MD) simulations have found that, with a longer side chain in perfluorosulfonic acid (PFSA) membranes, more clustering of the sulfonate groups and more local water−water clustering was present with a less-connected aqueous domain.68 A dissipative particle dynamics (DPD) study on the hydrated morphology of the 3M PFSA revealed that the longer side chain generally results in a stronger aggregation of water clusters and higher flexibility of the side chains.69 A comparison between Nafion and Hyflon revealed that the influence of the different sidechain lengths in the two polymers on the local atomic-level structure is relatively small.71 The phase separation in the PEMs is considered one of the major factors affecting their performance,60 and this phenomenon has been observed in a significant number of experimental and computational studies.56,57,72,73 An investigation of the morphology of cross-linked polytetrafluoroethylene (cPTFE) based radiation-grafted PEMs shows a phase separation with Received: May 16, 2012 Revised: August 19, 2012 Published: August 20, 2012 19104

dx.doi.org/10.1021/jp304783z | J. Phys. Chem. C 2012, 116, 19104−19116

The Journal of Physical Chemistry C

Article

one phase consisting of both water and poly(styrene sulfonic acid).74 It was also found that the water channels are narrower and the estimated percolation threshold is lower in sulfonated aromatic poly(ether ether ketone) (SPEEK) than in Nafion.75 A coarse-grained study of SPEEK and Nafion-based copolymers predicts a tricontinuous morphology, where the percolation of water and hydrophilic and hydrophobic blocks results in the formation of three cocontinuous phases.76 The influence of other factors on PEM morphology, such as water content, temperature, methanol concentration, and electric field, has been investigated. Mesoscale simulation of the proton-conducting membrane based on SPEEK shows that the microphase separation of hydrophilic and hydrophobic polymer chain units occurs even at small water contents.78 With an increase in methanol concentration, an increased solvent cluster size was observed in the hydrated Nafion.79 The presence of an electric field was observed to have an influence on the Nafion morphology as well.80 Several experimental and theoretical studies have also been devoted to providing a detailed relationship between the morphological features and proton conductivity.34,81−85 It was found that, at the same equivalent weight (EW), increasing the side-chain length would lead to higher proton diffusivity.81 The proton conductivity was found to increase linearly with an increase in the degree of sulfonation after a threshold of 52%.82 The water content and density of sulfonate groups are also suggested to be key factors for PT properties in sulfonic acid based PEMs.83 In an AIMD study of short-side-chain PFSA, the increase of side chain and backbone flexibility was found to be correlated with enhanced proton−water affinity.34 A coarsegrained study of the effects of molecular geometry on PT in PEM suggests that molecular geometry with hydrophilic blocks, encompassed with hydrophobic blocks, has relatively good compatibility with high proton conductivity and water resistance.84 With a coarse-grained methodology developed to describe mesoscopic PEM PT, it was shown that water content, network structure, local electronic field, and local PT behavior all modulate the overall proton conductivity of the PEM.85 With respect to the most studied PEM, Nafion, a reasonable degree of information regarding its morphological properties has been collected.70,86,87 However, the morphology of Nafion is still under debate, as a result of the fact that this ionomer has random chemical structures that can organize in the complex formation of crystalline and ionic domains.86 Large-scale MD simulations have been performed70 for the six most significant morphological models of Nafion, including the cluster-channel model (referred to as the cluster model in the following sections), the parallel cylinder model, the local order model, the lamellar model, the rod network model, and the random model. All of these models consist of water domains, surrounded by the Nafion polymer. A more detailed description of different models can be found in ref 70. The water part of the cluster-channel model consists of water spheres connected by water channels.56,70,86,88 The average water sphere and water channel diameters are ∼4 and ∼1 nm, respectively, and the mean distance between the centers of the two closest water spheres is ∼5 nm. The parallel cylinder model’s water domain consists of a series of aligned tubes filled with water; the water tubes are not connected by water bridges or rods.89 The proposed diameters of the water channels are distributed between 1.8 and 3.5 nm, with an average of 2.4 nm. This model is referred to as a cylinder model in the present study.

Figure 1. Molecular structure of Nafion used in the present study.

Figure 2. Upper left: snapshot of the cylinder-type model at hydration level 15; the axis of the cylinder is on the Z direction. Upper right: snapshot of the lamellar model at hydration level 15. Lower left and right: snapshots of the cluster model from two different perspectives; the clusters are connected on the Z axis. Red spheres represent the oxygen atoms, white spheres represent the hydrogen atoms, yellow spheres represent the sulfur atoms, and the rest represent the Nafion backbone atoms.

The water domains of the local order model consist of water spheres, also not connected by water rods or bridges.86,90,91 The size of the water spheres is similar to that of the clusterchannel model. The lamellar model consists of parallel slabs, which are filled with water and polymer molecules in an alternate manner.86,89,92,93 The water layer thickness ranges between ∼0.75 and ∼2.8 nm, while the thickness of the polymer layer is distributed between ∼2.1 and ∼3.5 nm. The water domains of the rod network model consist of random connected flexible tubes filled with water. 89,94,95 The diameters of the tubes are smaller than ∼3 nm. The random model does not involve any particular assumption for the Nafion structure.16 Interestingly, it was found that all nonrandom models exhibit the characteristic experimental ionomer scattering peak.70 In our previous study on proton solvation and transport in Nafion,54,55 the activation energy and PT as a function of the water loading levels was characterized by the self-consistent iterative multistate empirical valence bond (SCI-MS-EVB) reactive MD method.96 As a continuation of these earlier studies, the aim of the present research is to elucidate the influence of various possible Nafion morphologies on the PT behavior. By investigating PT in the lamellar, cylinder, and cluster-type models, it is shown in this work that the sulfonate 19105



Article

Table 1. Detailed Information for the Simulated Systemsa hydration level model water layer height of the lamellar models (nm)

6 0.8 [0.094 ± 0.021] 1.0 [0.068 ± 0.024]

water cylinder radius of the cylinder models (nm) spherical water cluster radius of the cluster models (nm)

1.05 [0.061 ± 0.015] 1.3 [0.0436 ± 0.0060]

10 1.0 [0.373 ± 0.071] 1.2 [0.228 ± 0.037] 1.0 [0.204 ± 0.050] 1.2 [0.239 ± 0.061] 1.35 [0.128 ± 0.029] 1.6 [0.195 ± 0.046]

15 1.0 1.2 1.0 1.2

[0.48 ± 0.098] [0.365 ± 0.082] [0.228 ± 0.048] [0.417 ± 0.099]

a

The excess proton diffusion constants (cm2/s) for each model are shown in the square brackets. All errors are estimated by the bootstrapping technique.

Figure 3. Distribution of the difference between the largest and the second-largest MS-EVB amplitudes Δc2 for excess protons in three distinct regions. For each proton, the shortest distance between that proton and all other sulfur atoms is l. The excess protons are then divided into three groups according to the value of l, such that l < 4.3 Å (black lines), 4.3 Å < l < 5.3 Å (red lines), and 5.3 Å < l (blue lines). (a) Hydration level 6 in the lamellar-type model; water layer height of 8 Å (solid lines) and 10 Å (dashed lines). (b) Hydration level 10 in the lamellar-type model; water layer height of 10 Å (solid lines) and 12 Å (dashed lines). (c) Hydration level 10 in the cylinder-type model; water cylinder radius of 10 Å (solid lines) and 12 Å (dashed lines). (d) Hydration level 15 in the cylinder-type model; water cylinder radius of 10 Å (solid lines) and 12 Å (dashed lines). (e) Hydration level 6 in the cluster model; water sphere radius of 10.5 Å (solid lines) and 13 Å (dashed lines). (f) Hydration level 10 in the cluster model; water sphere radius of 13.5 Å (solid lines) and 16 Å (dashed lines). 19106



Article

Figure 4. The change of Δc2 caused by increasing the water layer thickness or water cylinder radius, i.e., the dashed line minus the solid line in Figure 3. An identical coloring scheme to that in Figure 3 is adopted. (a) Hydration level 6 in the lamellar-type model. (b) Hydration level 10 in the lamellar-type model. (c) Hydration level 10 in the cylinder type-model. (d) Hydration level 15 in the cylinder-type model. (e) Hydration level 6 in the cluster model. (f) Hydration level 10 in the cluster model.

2. SIMULATION METHOD AND SETUP In previous studies, the SCI-MS-EVB methodology has been used to characterize the proton translocation behavior in the Nafion systems at different hydration levels and temperatures.21,54 In the present study, this method has been employed to further explore the proton solvation and transport (PS&T) phenomena of hydrated Nafion systems with different morphologies. The molecular mechanics potential for Nafion was adapted from Jang et al.16 The newly developed MS-EVB3 parameters,97 which include a more accurate underlying water model,98 have been adopted for the hydrated protons and water. The SCI-MS-EVB methodology was developed to effectively deploy the MS-EVB method in multiproton systems. In this approach, each proton and its associated water molecules are defined as an “EVB-complex”. A single proton MS-EVB problem is then solved in the effective fields created by all

group collective structures have a strong impact on the translocation behavior of protons. Hydrated excess protons in the lamellar-type model have higher diffusion rates than in the other two types of models, when the water layer thickness, water cylinder, and water sphere radius are small. Increasing the water layer thickness and water cylinder radius leads to the sulfonate groups being surrounded by more protons, while increasing the water sphere radius results in fewer protons being present around the sulfonate groups. Interestingly, the increases in the water layer thickness and the water cylinder radius have opposite effects on PT. The water sphere radius in the cluster model also changes the PT behavior, especially at low hydration levels. The remaining sections of this article are organized as follows: the simulation method and setup are provided in section 2, and the simulation results are discussed in section 3. Conclusions are drawn in section 4. 19107



Article

other complexes. Because of the interactions between complexes, the single proton problem is solved iteratively to generate the final consistent solution. The overall computational cost scales linearly with the number of excess protons.96 Nafion models with lamellar, cylinder, and cluster types at three distinct hydration levels were constructed (see Figure 2). Model cylindrical systems of three different radii and model lamellar systems of three different water slab thicknesses were characterized for the purpose of comparison. For the cluster models, spherical water clusters were connected by water cylinders with a diameter of 10 Å; two different water sphere sizes were constructed at each hydration level. More detailed information on these models may be seen in Table 1. The following procedures were employed in building and studying the model Nafion systems for all three types of systems. First, water domains with the desired geometric shapes were obtained from equilibrated simulations of water molecules. Forty distinct water molecules (64 for the cluster models with larger radii at each hydration level) were then randomly selected and converted to classical hydronium cations. A total of 40 Nafion 117 monomers (64 for the cluster models with larger radii at each hydration level), with their molecular structures (as shown in Figure 1), were then added to the selected water structures to form a simulation box. The sulfonate groups were arranged so that they were in close contact with the surfaces of the water domains. The simulation boxes were then equilibrated with standard MD to relax the Nafion molecules until the proper density was achieved. Position constraints were added to the water molecules and the hydronium cations in the relaxation simulation. In the second stage, constant NVT simulations of 8 ns were performed for each system, in which the C1 (Figure 1) atoms were fixed to maintain the proper geometric shape of the Nafion polymer, while still allowing reasonable fluctuations of the sulfonate groups. The additional restraint on the C1 atoms had a negligible impact on the proton translocation behavior for two reasons. First, the time scale of the polymer diffusion rate99 was much slower than the proton diffusion motion. Second, compared with the sulfonate groups, the positions of the C1 atoms were much farther from the hydrated protons in the models studied. All other atoms in the system were allowed to relax. The hydrated protonic species were represented by a classical hydronium force field during this relaxation phase. In the third, and final, data collection stage, the MS-EVB3 potential was employed to describe the Grotthuss shuttling hydrated protons in the production SCI-MS-EVB simulations and each system was equilibrated for 200 ps in the constant NVT ensemble, with the positions of the C1 atoms fixed. In each production run, 1 ns of SCI-MS-EVB simulation in the NVT ensemble at 300 K was performed for data collection and analysis (a total of 0.5 ns was collected for the cluster model at hydration level 10 with a water sphere radius of 16 Å). The van der Waal interaction cutoff was 10 Å.

Figure 5. (a) Proton diffusion coefficients for different morphological models. Lamellar model (solid lines), cylinder model (dashed lines), and cluster model (dotted lines), with hydration level 6 (black), hydration level 10 (red), and hydration level 15 (blue). (b) Meansquare displacements (MSDs) of hydration level 10. Lamellar model with a water layer height of 10 Å (black line), cylinder model with a water cylinder radius of 10 Å (red line), and cluster model with a water sphere radius of 13.5 Å (green line).

Δc2 = c21 − c22, as a function of the proton center of excess charge−sulfonate distance is plotted here in Figure 3. The center of excess charge (CEC) is used to characterize the location of the net positive charge defect associated with the hydrated excess proton (see, e.g., ref 97). The distribution Δc2 is scaled so that the relative amplitude corresponds to the relative number of protons. As noted,54,55 a difference Δc2 around 0 corresponds to the Zundel-like solvation structure and a difference larger than ∼0.4 corresponds to the Eigen-like structure. To elucidate the effects of increasing the water cylinder radius, water layer thickness, and water sphere radius, the difference in Δc2 distributions is depicted in Figure 4. This difference is defined as the dashed lines minus the solid lines in Figure 3. The black lines in Figure 4 represent the Δc2 distribution of an increased amount of protons in the neighborhood of the sulfonate groups. The positive values in the lamellar and cylinder models reveal that an increase in water cylinder radius and water layer height leads to more protons being present in the neighborhood of the sulfonate groups at all of the studied hydration levels. At a high hydration level of 15 in the cylinder models (Figure 4d), the peak of the black line shows that the increased amount of hydrated protons mainly consists of the Eigen-like structures. On the other hand, as the hydration level decreases to 10 (Figure 4c), the increased amount of protons consists of more Zundel-like structures.

3. RESULTS AND DISCUSSION 3.1. Proton Solvation Structure. The excess proton solvation structure was found to be influenced by the Nafion morphology adopted, and it was characterized by the MS-EVB amplitude distributions. These amplitudes provide information on the dominant structures of the hydrated protons, for example, the Eigen (H9O4+) or Zundel (H5O2+) cations in bulk water environments.54 The distribution of the difference between the largest two EVB amplitudes c21 and c22, that is, 19108



Article

Figure 6. The number of unique sulfonate groups visited by excess protons during the simulation: Panels (a) and (b) are the lamellar-type model, at hydration level 10, with an initial water layer thickness of 10 and 12 Å, respectively. Panels (c) and (d) represent the cylinder-type model, at hydration level 15, with an initial water cylinder radius of 10 and 12 Å, respectively.

Figure 7. The number of unique sulfonate groups visited by excess protons in the cluster model. Panels (a) and (b): hydration level 10 with an initial water sphere radius of 13.5 and 16 Å, respectively. Panels (c) and (d): hydration level 6, with a water sphere radius of 10.5 and 13 Å, respectively.

A similar pattern was found in the lamellar-type models as well (Figure 4a,b). In the lamellar model, a noticeable increase in

the hydronium-like structure was also found when the hydration level was decreased from 10 to 6 due to the high proton 19109



Article

concentration (as depicted by the peak ∼ 0.7 in Figure 4a). The hydronium-like structure is stabilized by several sulfonate groups nearby. Contrary to the lamellar and cylinder models, in the cluster model, the number of protons in the first solvation shell decreases as the radius of the sphere increases at both hydration levels 6 and 10 (Figure 4e,f). The proton solvation structure also changes when the water sphere radius increases. More protons are hydrated in the form of the Eigen-like structure in the region beyond the first solvation shell, as depicted by the sum of the green and red lines in Figure 4e,f. Therefore, in general, the difference in water−sulfonate morphology alters the number of hydrated protons around the sulfonate groups. Increasing the water cylinder radius and the water layer height results in more, not fewer, protons presenting in the neighborhood of the sulfonate groups and an increased proportion of Zundel-like hydrated structures at lower hydration levels. Increasing the water layer height also results in more hydronium-like hydrated structures presenting in the neighborhood of sulfonate groups at the lowest hydration level studied. On the other hand, the number of protons in the sulfonate groups’ first solvation shell decreases as the water sphere radius increases. The particular Nafion morphological character (i.e., lamellar, cylinder, cluster) influences the solvation structure change as well. The different proton distributions resulting from the changes in the water layer height, water cylinder, and water sphere radius can be explained by the competition between the attraction from the sulfonate groups and the geometric shape change of the water region. In all three models, increasing the water domain sizes leads to the water region outside of the first solvation shell of sulfonate groups occupying a greater volume of the whole water region. On the other hand, the excess protons are more easily trapped by the surrounding sulfonate groups with the increase in water domain size, since the interactions with other sulfonate groups is strongly weakened (i.e., screened by the additional exposure to water). Thus, it is more favorable for the protons to visit the water region beyond the sulfonate groups’ first solvation shell from the perspective of entropy, but it is more energetically favorable for the proton to associate with a particular sulfonate group. For the cluster models, the region beyond the sulfonate groups’ first solvation shell increases more rapidly than that for the other two models, therefore, leading to the decrease in proton population in the first solvation shell with an increase in the water cluster size. 3.2. Proton Diffusion. To characterize the PT behavior, excess proton CEC diffusion coefficients were first calculated, as shown in Figure 5a (the errors are listed in Table 1 so that the view of Figure 5a is not blocked). The excess proton diffusion coefficients were estimated as 1/6 of the mean-square displacement (MSD) slope from 50 to 600 ps. The MSD plot of hydration level 10 in three different models is depicted in Figure 5b. In the lamellar model, as the height of the water layer increases, the rates of proton diffusion decrease at all of the hydration levels examined. By contrast, the diffusion rates increase in the cylinder models as the radius of the water cylinder increases. Furthermore, as shown in Figure 5, the hydrated protons exhibit faster diffusion in the lamellar model with a thickness of 10 Å, than in the cylinder model with a radius of 10 Å. At hydration levels of 10 and 15, as the water layer thickness and the cylinder radius increase, the proton diffusion rates converge to the same value; that is, the cylinderand lamellar-type models tend to have the same diffusion rates

Figure 8. Trajectory of a single excess proton motion during 500 ps of SCI-MS-EVB MD simulation, by using the identical coloring scheme to that in Figure 2. Panels (a) and (b) are typical proton transport patterns in the lamellar-type model, at hydration level 10, with a water layer thickness of 10 and 12 Å, respectively. Panel (c) is a typical proton transport pattern in the cylinder model at hydration level 10, with the water cylinder radius being 10 Å. The white dots represent the beginning of the simulation trajectory, and deep blue dots represent the end of the simulation trajectory.

at a water cylinder radius of 12 Å and a water slab thickness of 12 Å. The diffusion rates in the cluster model are, however, different from the lamellar and cylinder models. With an increase in the water sphere size, the diffusion rates decrease at hydration level 6, but they actually increase at hydration level 10. This diverse behavior of the diffusion rates indicates that the PT mechanisms are rather different at the two hydration levels. The differences in the diffusion rates were found to be strongly connected with the PT patterns. Lamellar and cylinder models exhibit different PT patterns, which leads to the different behaviors of the diffusion constant change as the water layer height/water cylinder radius increases (Figure 5a). The complex behavior of diffusion coefficients in the cluster model is correlated with the PT pattern change from low to high hydration levels. A more detailed insight into the relationship between the diffusion coefficient change and the PT pattern can be gained from the PT patterns analysis, as discussed in section 3.3. 19110



Article

Figure 9. Trajectory of a single excess proton’s motion in the cluster model, by using the identical coloring scheme to that in Figure 2. Upper left and right are typical proton transport patterns during 1 ns at hydration level 6, with a water sphere radius of 10.5 and 13 Å, respectively. Lower left and right are typical proton transport patterns at hydration level 10, with a water sphere radius of 13.5 Å (1 ns) and 16 Å (500 ps), respectively. In each trajectory, the color gradually changes from pink to white to deep blue with time. The red circle depicts the approximate center of the cluster domain.

3.3. Proton Transport Patterns. 3.3.1. Lamellar Model. The PT pattern has been analyzed to explain the difference in proton diffusion coefficients in water slabs of different thicknesses. Figure 8a shows the 500 ps trajectory of a single excess proton CEC, in the lamellar model with an initial water layer thickness of 10 Å, in which both Nafion/water interfaces (i.e., the water layer’s two surfaces are in contact with a Nafion layer) are visited by this single protonic charge. The middle of the water slab is visited by the proton as well, during the transport process. However, as the water slab thickness increases, the probability for one single proton to visit the central bulk water region is significantly reduced. Instead, the protonic charge has the tendency to transport along the surface of Nafion, where the sulfonate groups are located, as depicted in Figure 8b.

To investigate further the proton diffusion behavior at each hydration level during the 1 ns simulation time, the numbers of different sulfonate groups visited by protons are plotted in Figures 6 and 7. The definition of visiting is the same as described in our previous study.54 At each hydration level, the number of sulfonate groups visited is correlated with the proton diffusion rates shown in Figure 5; that is, the higher the proton diffusion rates, the more sulfonate groups visited. This result indicates a strong correlation between the distribution of the sulfonate groups and proton mobility, which is consistent with the conclusion from our previous study that proton mobility was the major factor in deciding the numbers of sulfonate groups visited.54 In the next section, the role of the sulfonate groups in the PT events, in the different Nafion models, is discussed to provide a more detailed picture for the PT−sulfonate group relationship. 19111



Article

regions are more localized and less connected. Localization effects of protons near one or more sulfonate groups are observed, which leads to the slow diffusive motion. Further insight into the PT in the lamellar-type model can be gained by considering the sulfonate group distribution. In terms of the geometric factor, the constraints on the directions perpendicular to the water layer (z axis) are reduced as the water layer thickness increases. The energy factor exhibits the opposite effect, which leads to more well-connected regions in the model with lower water layer heights and less-connected regions in the model with greater water layer heights, as shown in Figure 10. 3.3.2. Cylinder Model. As opposed to the lamellar model, in the cylinder model, excess protons have a tendency for translocation at the Nafion−water interface; that is, the hydrated protons have a smaller probability of visiting the central water region, and they remain close to the neighborhood of the sulfonate groups. To provide a physical picture of the proton translocation, a 500 ps trajectory of a single excess proton is shown in Figure 8c. Similar to the lamellar model, the energy factor and the geometric factor are employed to explain the change in proton diffusion when altering the water cylinder radius. As the radius of the water cylinder increases, the constraints in the plane perpendicular to the cylinder principal axis (XY plane as defined in the upper left panel of Figure 2) are reduced. However, more protons are present in the neighborhood of the sulfonate groups and a more localized and less-connected two-dimensional distribution is observed in the XY plane (Figures 11a and 11b). Therefore, the faster diffusive motions of protons in the model with a larger water cylinder radius are not caused by the reduced geometric constraints in the XY plane. A further examination of the proton distribution in the XZ plane (Figure 11c,d) shows better-connected regions in the model with larger water cylinder radii. This leads to the faster diffusive motion for the protons. The underlying mechanism for these better-connected regions can be rationalized by noting the closer distance between the sulfonate groups, as depicted in the radial distribution functions (RDFs) (Figure 12). Similar to the situation in the lamellar model, the closer distance between the sulfonate groups in the hydrophilic region surface helps facilitate the proton transport from one sulfonate group to another by reducing the corresponding energy barriers caused from the electrostatic interactions. With this information in hand, we can now explore the reasons behind the qualitative difference in the diffusion behavior between the cylinder and lamellar models. From a starting point of a cylinder radius of 10 Å or a lamellar height of 10 Å, we can see that the protons are more free to explore the center of the hydrophilic region in the lamellar model for a given hydration level. Since protons are able to cross from interface to interface and are less hindered by electrostatic attraction of the sulfonates when they are in the center, we see a higher diffusion rate in the lamellar model. Increasing the interfacial spacing produces two different outcomes in the two models. As the protons in the cylindrical model are already trapped at the interface, increasing the radius has very little effect on their diffusion rate away from the interface, but the closer spacing of the sulfonates on the surface increases the probability of proton exchange between sulfonates on the interface, thus increasing proton mobility. Increasing the water layer height in the lamellar model, however, serves to limit the excess protons' presence in the center of the hydrophilic region, thus trapping the protons to a given interface. Any gain in diffusion rate from the closer spacing of sulfonate groups is

Similar to the translocation of ions in confined environments,100,101 the dynamics of protons in hydrated Nafion can be characterized by (1) the morphology of Nafion, (2) the interaction energy between the Nafion and excess protons, and (3) the interaction between the protonic excess charges. In the simulation of lamellar models with two different water layer heights, the proton concentrations are at an identical level. Therefore, the change in proton−proton interaction in different water layer thicknesses is largely induced by the change of geometric shape and the Nafion/proton potential energy interactions. Therefore, the change in dynamical properties caused by changing the water slab thickness can be rationalized by changes in geometric shape and Nafion/proton potential energy interactions, which are referred to as “geometric factor” and “energy factor” in the following sections. Figure 10 shows the two-dimensional probability distribution of protons in each system. In the lamellar model with an initial water slab thickness of 10 Å (Figure 10a), the geometric regions preferred by protons are well-connected, which facilitates the proton transport between different regions, thereby leading to higher rates of diffusive motion. With an increase of the water layer height (Figure 10b), the preferred

Figure 10. Two-dimensional probability distributions of the excess protons in the lamellar-type models in the XZ plane, with directions defined as in Figure 2. Panels (a) and (b) represent the lamellar model with an initial water layer thickness of 10 and 12 Å, and both distributions are calculated at hydration level 10. 19112



Article

Figure 11. Two-dimensional excess proton probability distributions in the cylinder models in the XY plane (a, b) and the XZ plane (c, d) with directions defined in Figure 2. Panels (a) and (b) show the cylinder model at hydration level 15, with an initial water cylinder radius of 10 and 12 Å, respectively. Panels (c) and (d) are the cylinder model at hydration level 15, with a starting water cylinder radius of 10 and 12 Å, respectively.

3.3.3. Cluster Model. The PT in the cluster models exhibits different behaviors at the two studied hydration levels. At low hydration level 6, with an increase in the water sphere radius, the diffusion rate decreases (Figure 5), similar to the lamellar system’s behavior. A detailed analysis of the PT trajectory shows that, in the model with a radius of 10.5 Å, the protons can move across the central water regions more frequently than those in the model with a radius of 13 Å. This is similar to the PT in the lamellar model, as discussed previously. The interaction with the sulfonate groups on the opposite side of the sphere facilitates the proton moving out of the first solvation shell of the closest sulfonate group and crossing the accessible central water regions. When the water sphere radius increases to 13 Å, the interaction between the proton and sulfonate groups on the other side is weakened. Therefore, the protons tend to move along the sulfonate groups’ surface, as seen in the upper right panel of Figure 9. At hydration level 10 with an increase in the water sphere radius, the diffusion rate increases (Figure 5), which is now similar to the behavior in the cylindrical system. The probability of visiting the central water region is low for both of the radii,

Figure 12. Radial distribution functions (RDFs) between the sulfur atoms at hydration levels 10 (black lines) and 15 (red lines) in the cylinder model. The solid lines and dashed lines represent systems with water cylinder radii of 10 and 12 Å, respectively.

more than counteracted by the reduction in proton mobility within the center of the hydrophilic region, and so we see a decrease in diffusion rate. 19113



Article

rates. For clusters at hydration level 10, the protons mainly transport along the sulfonate group surfaces and a larger radius leads to faster proton diffusion. The results presented in this paper, taken together, suggest that there may be significant payoffs from synthetic efforts devoted to modulating and optimizing the morphology of proton-conducting membrane materials, subject to other constraints, such as thermal and mechanical stability.

since the cluster size is large enough to reduce the interaction between the proton and the sulfonate groups on the opposite side of the spherical cluster, again much like how the protons are trapped at the interface in the cylindrical model. The main transport pathway is along the inner surface of the sphere. When the spherical cluster radius increases, the spacing between the sulfonates is reduced. Therefore, the proton can more easily traverse the interface, leading to a higher diffusion rate. The apparent complex behavior of the cluster model can now be seen to reduce to the same factors as we have discussed before in the previous two models. Diffusion rates are low if the proton is trapped on the interface, but diffusion rates increase either if the proton can cross the hydrophilic region or if the spacing between the sulfonates on the interface decreases. One factor that is noticeably different in the cluster model, however, is the overall lower rates of diffusion compared with those of the other two models (Figure 5). This can be rationalized by the fact that the water region in this model is effectively an isolated sphere, which limits the distances the protons can travel compared to the other models, which are each essentially infinitely long in at least one direction. It was also observed that, when the proton is transported to the channel region of the cluster model, it could be trapped locally, which further reduces the diffusion rate of the system as a whole.

■

AUTHOR INFORMATION

Corresponding Author

*Fax: 773-795-9106. Phone: 773-702-9092. E-mail: gavoth@ uchicago.edu. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

This research was supported by the Department of Energy (DOE), Office of Basic Energy Sciences (BES), Division of Chemical Sciences, Geosciences, and Biosciences, through Grant No. DE-FG02-10ER16171.

(1) Giffin, G. A.; Piga, M.; Lavina, S.; Navarra, M. A.; D’Epifanio, A.; Scrosati, B.; Di Noto, V. J. Power Sources 2012, 198, 66−75. (2) Yoonessi, M.; Heinz, H.; Dang, T. D.; Bai, Z. W. Polymer 2011, 52, 5615−5621. (3) Buquet, C. L.; Hamonic, F.; Saiter, A.; Dargent, E.; Langevin, D.; Nguyen, Q. T. Thermochim. Acta 2010, 509, 18−23. (4) Hongsirikarn, K.; Mo, X. H.; Goodwin, J. G.; Creager, S. J. Power Sources 2011, 196, 3060−3072. (5) Kishi, A.; Inoue, M.; Umeda, M. J. Phys. Chem. C 2010, 114, 1110−1116. (6) Krivobokov, I. M.; Gribov, E. N.; Okunev, A. G.; Spoto, G.; Parmon, V. N. Solid State Ionics 2010, 180, 1694−1701. (7) Lu, J. L.; Lu, S. F.; Jiang, S. P. Chem. Commun. (Cambridge, U.K.) 2011, 47, 3216−3218. (8) Matos, B. R.; Santiago, E. I.; Rey, J. F. Q.; Ferlauto, A. S.; Traversa, E.; Linardi, M.; Fonseca, F. C. J. Power Sources 2011, 196, 1061−1068. (9) Molla, S.; Compan, V. J. Power Sources 2011, 196, 2699−2708. (10) Silberstein, M. N.; Pillai, P. V.; Boyce, M. C. Polymer 2011, 52, 529−539. (11) Sugawara, T.; Kawashima, N.; Murakami, T. N. J. Power Sources 2011, 196, 2615−2620. (12) Vogel, C.; Komber, H.; Quetschke, A.; Butwilowski, W.; Pötschke, A.; Schlenstedt, K.; Meier-Haack, J. React. Funct. Polym. 2011, 71, 828−842. (13) Zhao, Q. A.; Majsztrik, P.; Benziger, J. J. Phys. Chem. B 2011, 115, 2717−2727. (14) Zhao, Y.; Jiang, Z.; Xiao, L.; Xu, T.; Wu, H. J. Power Sources 2011, 196, 6015−6021. (15) Eikerling, M.; Paddison, S. J.; Zawodzinski, T. A., Jr. J. New Mater. Electrochem. Syst. 2002, 5, 1−9. (16) Jang, S. S.; Molinero, V.; Cagin, T.; Goddard, W. A. J. Phys. Chem. B 2004, 108, 3149−3157. (17) Paul, R.; Paddison, S. J. J. Phys. Chem. B 2004, 108, 13231− 13241. (18) Paddison, S. J.; Elliott, J. A. J. Phys. Chem. A 2005, 109, 7583− 7593. (19) Paul, R.; Paddison, S. J. J. Chem. Phys. 2005, 123, 224704. (20) Petersen, M. K.; Wang, F.; Blake, N. P.; Metiu, H.; Voth, G. A. J. Phys. Chem. B 2005, 109, 3727−3730. (21) Petersen, M. K.; Voth, G. A. J. Phys. Chem. B 2006, 110, 18594− 18600.

4. CONCLUSIONS In the present work, the influence of Nafion morphology on PS&T was investigated in detail. The solvation structure change as a function of the distance to the sulfonate groups was analyzed at different hydration levels and for different morphological structures. Increasing the water cylinder radius or water layer height leads to more protons being present around the sulfonate groups. Furthermore, at lower hydration levels, these extra protons around the sulfonate group tend to consist of more Zundel-like structures, which is also influenced by the distinct morphological structures of Nafion. In the cluster models, the number of protons in the sulfonate groups’ first solvation shell decreases as the water sphere size increases. The proton distribution difference in the three types of models can be rationalized by the competition of attraction from sulfonate groups and changes in the geometric shape of the water region. Interestingly, the lamellar-, cylinder-, and cluster-type models exhibit distinct PT patterns and resulted in different proton diffusion rates as the particular Nafion morphology is modulated. In all models, the electrostatic interaction between the sulfonate groups and protons is likely the key factor in determining the change of PT behavior. Particularly for the lamellar-type models, the electrostatic interaction change caused by the distance between the sulfonate groups across the hydrophilic region is a vital factor in determining the proton diffusion rates. For the cylinder-type models, the reduced distance between the sulfonate groups on the amphiphilic interface upon increasing the cylinder radius results in a lower energy barrier for PT between sulfonate groups, thereby leading to higher proton diffusion rates. With respect to the cluster models, the proton transport pattern correlates to cluster size. For the smallest cluster studied (hydration level 6 with radius of 10.5 Å), the protons can cross the central water region more frequently than the model with a larger radius at the same hydration level, thereby leading to slightly higher diffusion 19114



Article

(22) Devanathan, R.; Venkatnathan, A.; Dupuis, M. J. Phys. Chem. B 2007, 111, 8069−8079. (23) Devanathan, R.; Venkatnathan, A.; Dupuis, M. J. Phys. Chem. B 2007, 111, 13006−13013. (24) Glezakou, V.-A.; Dupuis, M.; Mundy, C. J. Phys. Chem. Chem. Phys. 2007, 9, 5752−5760. (25) Venkatnathan, A.; Devanathan, R.; Dupuis, M. J. Phys. Chem. B 2007, 111, 7234−7244. (26) Petersen, M. K.; Hatt, A. J.; Voth, G. A. J. Phys. Chem. B 2008, 112, 7754−7761. (27) Cui, S.; Liu, J.; Selvan, M. E.; Paddison, S. J.; Keffer, D. J.; Edwards, B. J. J. Phys. Chem. B 2008, 112, 13273−13284. (28) Hristov, I. H.; Paddison, S. J.; Paul, R. J. Phys. Chem. B 2008, 112, 2937−2949. (29) Devanathan, R.; Venkatnathan, A.; Rousseau, R.; Dupuis, M.; Frigato, T.; Gu, W.; Helms, V. J. Phys. Chem. B 2010, 114, 13681− 13690. (30) Yana, J.; Nimmanpipug, P.; Chirachanchai, S.; Gosalawit, R.; Dokmaisrijan, S.; Vannarat, S.; Vilaithong, T.; Lee, V. S. Polymer 2010, 51, 4632−4638. (31) Hofmann, D. W. M.; Kuleshova, L. N.; D’Aguanno, B. J. Power Sources 2010, 195, 7743−7750. (32) Choe, Y. K.; Tsuchida, E.; Ikeshoji, T.; Ohira, A.; Kidena, K. J. Phys. Chem. B 2010, 114, 2411−2421. (33) Fernandez-Carretero, F. J.; Suarez, K.; Solorza, O.; Riande, E.; Compan, V. J. New Mater. Electrochem. Syst. 2010, 13, 191−199. (34) Ahadian, S.; Mizuseki, H.; Kawazoe, Y. J. Membr. Sci. 2011, 369, 339−349. (35) Dorenbos, G.; Morohoshi, K. J. Chem. Phys. 2011, 134, 044133. (36) Francia, C.; Ijeri, V. S.; Specchia, S.; Spinelli, P. J. Power Sources 2011, 196, 1833−1839. (37) Ban, S.; Huang, C.; Yuan, X.-Z.; Wang, H. J. Phys. Chem. B 2011, 115, 11352−11358. (38) Berning, T.; Odgaard, M.; Kaer, S. K. J. Power Sources 2011, 196, 6305−6317. (39) Eikerling, M. H.; Berg, P. Soft Matter 2011, 7, 5976−5990. (40) Habenicht, B. F.; Paddison, S. J. J. Phys. Chem. B 2011, 115, 10826−10835. (41) Ilhan, M. A.; Spohr, E. J. Phys.: Condens. Matter 2011, 23, 234104. (42) Keffer, D. J.; Selvan, M. E.; Calva-Munoz, E. J. Phys. Chem. B 2011, 115, 3052−3061. (43) Ohkubo, T. O. T.; Kidena, K.; Takimoto, N.; Ohira, A. J. Mol. Model. 2011, 17, 739−755. (44) Parsafar, G.; Abroshan, H.; Akbarzadeh, H.; Taherkhani, F. Mol. Phys. 2011, 109, 709−724. (45) Sagarik, K.; Phonyiem, M.; Chaiwongwattana, S.; Lao-ngam, C. Phys. Chem. Chem. Phys. 2011, 13, 10923−10939. (46) Sebastiani, D.; Luduena, G. A.; Kuhne, T. D. Chem. Mater. 2011, 23, 1424−1429. (47) Selvan, M. E.; Keffer, D. J.; Cui, S. J. Phys. Chem. C 2011, 115, 18835−18846. (48) Tavakoli, B.; Roshandel, R. Renewable Energy 2011, 36, 3319− 3331. (49) Zhang, Z.; Promislow, K.; Martin, J.; Wang, H.; Balcom, B. J. J. Power Sources 2011, 196, 8525−8530. (50) Allahyarov, E.; Taylor, P. L.; Lowen, H. J. Phys.: Condens. Matter 2011, 23, 455102. (51) Kim, Y. G.; Bae, Y. C. J. Polym. Sci., Part B: Polym. Phys. 2011, 49, 1455−1463. (52) Ohkubo, T.; Iwadate, Y.; Kim, Y. S.; Henson, N.; Choe, Y. K. Theor. Chem. Acc. 2011, 130, 555−561. (53) Yu, T. H.; Sha, Y.; Liu, W.-G.; Merinov, B. V.; Shirvanian, P.; Goddard, W. A. J. Am. Chem. Soc. 2011, 133, 19857−63. (54) Feng, S.; Voth, G. A. J. Phys. Chem. B 2011, 115, 5903−5912. (55) Feng, S.; Voth, G. A. J. Phys. Chem. B 2011, 115, 10570−10570. (56) Elliott, J. A.; Wu, D.; Paddison, S. J.; Moore, R. B. Soft Matter 2011, 7, 6820−6827.

(57) Clark, J. K., II; Paddison, S. J.; Eikerling, M.; Dupuis, M.; Zawodzinski, T. A., Jr. J. Phys. Chem. A 2012, 116, 1801−1813. (58) Spohr, E.; Commer, P.; Kornyshev, A. A. J. Phys. Chem. B 2002, 106, 10560−10569. (59) Choe, Y.-K.; Tsuchida, E.; Ikeshoji, T.; Yamakawa, S.; Hyodo, S.-a. J. Phys. Chem. B 2008, 112, 11586−11594. (60) Park, C. H.; Lee, C. H.; Sohn, J. Y.; Park, H. B.; Guiver, M. D.; Lee, Y. M. J. Phys. Chem. B 2010, 114, 12036−12045. (61) Kong, C. S.; Kim, D. Y.; Lee, H. K.; Shul, Y. G.; Lee, T. H. J. Power Sources 2002, 108, 185−191. (62) Nieh, M. P.; Guiver, M. D.; Kim, D. S.; Ding, J. F.; Norsten, T. Macromolecules 2008, 41, 6176−6182. (63) Idupulapati, N.; Devanathan, R.; Dupuis, M. J. Phys. Chem. B 2011, 115, 2959−2969. (64) Abroshan, H.; Akbarzadeh, H.; Taherkhani, F.; Parsafar, G. Mol. Phys. 2010, 108, 3393−3404. (65) Mahajan, C. V.; Ganesan, V. J. Phys. Chem. B 2010, 114, 8357− 8366. (66) Mahajan, C. V.; Ganesan, V. J. Phys. Chem. B 2010, 114, 8367− 8373. (67) Lins, R. D.; Devanathan, R.; Dupuis, M. J. Phys. Chem. B 2011, 115, 1817−1824. (68) Liu, J. W.; Suraweera, N.; Keffer, D. J.; Cui, S. T.; Paddison, S. J. J. Phys. Chem. C 2010, 114, 11279−11292. (69) Wu, D. S.; Paddison, S. J.; Elliott, J. A.; Hamrock, S. J. Langmuir 2010, 26, 14308−14315. (70) Knox, C. K.; Voth, G. A. J. Phys. Chem. B 2010, 114, 3205− 3218. (71) Karo, J.; Aabloo, A.; Thomas, J. O.; Brandell, D. J. Phys. Chem. B 2010, 114, 6056−6064. (72) Brandell, D.; Karo, J.; Thomas, J. O. J. Power Sources 2010, 195, 5962−5965. (73) Wu, D. S.; Paddison, S. J.; Elliott, J. A. Macromolecules 2009, 42, 3358−3367. (74) Sawada, S.; Yamaki, T.; Ozawa, T.; Suzuki, A.; Terai, T.; Maekawa, Y. Kobunshi Ronbunshu 2010, 67, 224−227. (75) Komarov, P. V.; Veselov, I. N.; Chu, P. P.; Khalatur, P. G.; Khokhlov, A. R. Chem. Phys. Lett. 2010, 487, 291−296. (76) Komarov, P. V.; Veselov, I. N.; Chu, P. P.; Khalatur, P. G. Soft Matter 2010, 6, 3939−3956. (77) Dorenbos, G.; Morohoshi, K. Energy Environ. Sci. 2010, 3, 1326−1338. (78) Komarov, P. V.; Veselov, I. N.; Khalatur, P. G. Polym. Sci., Ser. A 2010, 52, 191−208. (79) Abroshan, H.; Akbarzadeh, H.; Taherkhani, F.; Parsafar, G. Mol. Phys. 2011, 109, 709−724. (80) Allahyarov, E.; Taylor, P. L.; Lowen, H. Phys. Rev. E 2010, 81, 031805. (81) Allahyarov, E.; Taylor, P. L. J. Polym. Sci., Part B: Polym. Phys. 2011, 49, 368−376. (82) Paik, Y.; Chae, S. A.; Han, O. H.; Hwang, S. Y.; Ha, H. Y. Polymer 2009, 50, 2664−2673. (83) Marschall, R.; Tolle, P.; Cavalcanti, W. L.; Wilhelm, M.; Kohler, C.; Frauenheim, T.; Wark, M. J. Phys. Chem. C 2009, 113, 19218− 19227. (84) Kinjo, T.; Murayama, Y.; Imai, K.; Munekata, T.; Hyodo, S. Kobunshi Ronbunshu 2010, 67, 187−191. (85) Jorn, R.; Voth, G. A. J. Phys. Chem. C 2012, 116, 10476−10486. (86) Mauritz, K. A.; Moore, R. B. Chem. Rev. (Washington, DC, U.S.) 2004, 104, 4535−4585. (87) Bass, M.; Berman, A.; Singh, A.; Konovalov, O.; Freger, V. J. Phys. Chem. B 2010, 114, 3784−3790. (88) Hsu, W. Y.; Gierke, T. D. J. Membr. Sci. 1983, 13, 307−326. (89) Schmidt-Rohr, K.; Chen, Q. Nat. Mater. 2008, 7, 75−83. (90) Rollet, A. L.; Gebel, G.; Simonin, J. P.; Turq, P. J. Polym. Sci., Part B: Polym. Phys. 2001, 39, 548−558. (91) Gebel, G.; Moore, R. B. Macromolecules 2000, 33, 4850−4855. (92) Haubold, H. G.; Vad, T.; Jungbluth, H.; Hiller, P. Electrochim. Acta 2001, 46, 1559−1563. 19115



Article

(93) Krivandin, A. V.; Solov’eva, A. B.; Glagolev, N. N.; Shatalova, O. V.; Kotova, S. L. Polymer 2003, 44, 5789−5796. (94) Kreuer, K. D. J. Membr. Sci. 2001, 185, 29−39. (95) Kim, M. H.; Glinka, C. J.; Grot, S. A.; Grot, W. G. Macromolecules 2006, 39, 4775−4787. (96) Wang, F.; Voth, G. A. J. Chem. Phys. 2005, 122, 144105. (97) Wu, Y. J.; Chen, H. N.; Wang, F.; Paesani, F.; Voth, G. A. J. Phys. Chem. B 2008, 112, 467−482. (98) Wu, Y. J.; Tepper, H. L.; Voth, G. A. J. Chem. Phys. 2006, 124, 024503. (99) Bucknall, D. G.; Butler, S. A.; Higgins, J. S. Macromolecules 1999, 32, 5453−5456. (100) Levitt, D. G. Biophys. J. 1991, 59, 278−288. (101) Levitt, D. G. Biophys. J. 1991, 59, 271−277.

19116


Effects of Polymer Morphology on Proton Solvation and Transport in

Recommend Documents