Viscosity of Nafion Oligomers as a Function of Hydration and

Subscriber access provided by ONDOKUZ MAYIS UNIVERSITESI

Article

Viscosity of Nafion Oligomers as a Function Of Hydration and Counterion Type: A Molecular Dynamics Study Kevin Brendan Daly, Athanassios Z. Panagiotopoulos, Pablo G. Debenedetti, and Jay B Benziger J. Phys. Chem. B, Just Accepted Manuscript • Publication Date (Web): 12 Nov 2014 Downloaded from http://pubs.acs.org on November 13, 2014

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Viscosity of Nafion Oligomers as a Function of Hydration and Counterion Type: a Molecular Dynamics Study Kevin B. Daly, Athanassios Z. Panagiotopoulos, Pablo G. Debenedetti, and Jay B. Benziger∗ Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544 E-mail: [email protected]

∗

To whom correspondence should be addressed

1

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract The design of fuel cells and lithium-ion batteries is constrained, in part, by mechanical creep and perforation of the polymer electrolyte, a process that is poorly understood at the molecular level. The mechanical stiffness (quantified as shear viscosity) and structure of a widely used polymer electrolyte, Nafion, are studied in the limit of a low solvent volume fraction (≤ 26% v/v H2 O) using molecular dynamics simulations. The viscosity is shown to increase by up to 4 orders of magnitude in response to changes in composition representing as little as 2 wt% of system. Two types of compositional changes are considered: changes in solvent volume fraction and counterion type. A system with a counterion Xv+ for every v Nafion monomers and y water molecules is denoted as (RSO3 )v X · (H2 O)y . The following trend is observed in viscosity: (RSO3 )2 Ca > RSO3 Na > RSO3 H · (H2 O)3 > RSO3 H ≈ RSO3 H · (H2 O)10 . This trend correlates with changes in the strength of the SO3– /Xv+ /SO3– cross-links and the size of the cross-link networks. Counterion type is shown to strongly influence the morphology. The simulations are able to reproduce some important experimental trends without crystalline domains or high-MW effects like entanglements, providing a simplified understanding of the mechanical properties of Nafion.

Keywords: Nafion, fuel cells, polymer membranes, plasticizers, ionic crosslinks, mechanical properties

2


Page 2 of 37

Page 3 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


1

Introduction

Polymer electrolytes are an attractive alternative to conventional cation conductors in devices that store and release energy electrochemically. In lithium-ion batteries, thin-film polymer electrolytes take the place of toxic, flammable organic solvents. 1 In hydrogen fuel cells, polymer electrolytes permit operation at high power densities (∼0.7 W/cm2 ) and mild temperatures (< 373 K). 2 In these devices, the polymer chains have a MW that is high enough to make the material a solid, imparting mechanical integrity to the electrolyte that prevents direct contact between the electrodes. 1 Mechanical integrity is particularly important for fuel cells, in which the membrane is compressed between the electrodes to minimize contact resistance, leaving residual stresses of around 5 bar under typical operating conditions. 3 Under these conditions, fuel cell membranes must resist perforation, since pinhole leaks permit the cross-over of reactant gases. The resulting parasitic reactions would not only lower the fuel cell voltage, but also create a thermal hot spot. 4 This hot spot would locally accelerate the mechanical degradation of the membrane, starting a destructive cycle. 4 Mechanical degradation constrains the design of fuel cells in many ways. Ideally, fuel cell membranes should be as thin as possible to minimize ionic resistance and material costs. However, thinner membranes have been shown to fail sooner. For example, the lifetime of the most widely used membrane, Nafion, is O(10000 h) for a thickness of O(100 μm), but only O(1000 h) for a thickness of O(10 μm). 5 The lifetime also decreases with decreasing equivalent weight (EW), which is a measure of the mass of polymer per mole of ionizable sites. 5 A lower EW is desirable since it enhances water uptake, and therefore proton conductivity. 6 Particularly detrimental to mechanical integrity are temperatures above 373 K, where membranes dehydrate and lifetimes plummet to O(100 h) 5,7 . On the other hand, high temperatures are ideal for accelerating reaction kinetics and minimizing CO poisoning of the catalyst. 8,9 . At high temperatures, waste heat is also rejected more easily, because the rate of heat transfer is proportional to the difference between the ambient and operating 3



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

temperatures. 10 The importance of mechanical degradation has motivated fundamental research into the mechanical properties of Nafion. These properties exhibit a rich dependence on temperature, cation type, and water activity (aw = fw /fwsat , where fw is the fugacity of water and fwsat is the fugacity of pure water at saturation and at the same temperature). At aw = 0, the elastic modulus of dry Nafion is observed to drop by two orders of magnitude between 333 and 353 K. 11,12 This abrupt decrease has been attributed to an order/disorder transition involving the ionizable sulfonic acid groups of the chains. 12 At higher humidities, X-ray and neutron scattering indicate that these hydrophilic groups form nanoscale clusters with water molecules. These clusters are hypothesized to persist all the way to aw = 0, and to break apart at sufficiently high temperature. 12 Unfortunately, this hypothesis cannot be easily verified with X-ray or neutron scattering experiments on dry membranes, because of the low contrast in electron or hydrogen density, respectively. The lack of contrast makes the clusters difficult to detect at low aw . The contrast can be improved by replacing the H+ cation with Cs+ or N(CH3 )4+ , respectively, 13–16 but the impact of this change on the morphology of the clusters at low aw is unknown. The temperature range of the hypothesized order/disorder transition is close to an α relaxation peak at T = 383 K in the loss tangent (tan δ), a property measured using dynamic mechanical analysis (DMA). 13 This temperature is sometimes interpreted as the “Tg ”, although two other relaxation peaks exist at ∼173 K (γ peak) and ∼293 K (β peak). 13 These latter two peaks are attributed to transitions in the perfluorocarbon chains. 13 The α peak shifts to higher temperatures as water is added to membrane. 17,18 This peak also jumps to 523 K if the cation is changed from H+ to Na+ , indicating a strong stiffening effect. This effect can also be seen in the storage 19 and elastic 20 moduli, both of which increase by two orders of magnitude upon swapping the cation at T ≈ 373 K. It is unclear whether the increase gets even larger at T ≫ 373 K, since these temperatures lie outside the experimental data range for the H+ system (e.g. see Figure S9 in the Supporting Information). 19,20 At

4


Page 4 of 37

Page 5 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


T > 473 K, the storage modulus can be increased further by at least an order of magnitude if the cation is switched from Na+ to Ca2+ . 19 The dry conditions described above are in practice only encountered during the start-up and shut-down of the fuel cell. 21 At steady state, the fuel cell is operated at high humidity, since the proton conductivity of Nafion is known to increase by nearly five orders of magnitude from aw = 0 to aw = 1. 22 Adding water increases the elongational creep of the membrane at T . 323 K; 12 in other words, water behaves as a plasticizer, as would be expected for most polymer/solvent systems. 23 However, at T & 323 K, water stiffens the membrane at aw = 0 − 0.05, and then plasticizes the membrane at higher humidities. 12 Water is known to increase the size of hydrophilic clusters at higher humidity, 13 so it is hypothesized that at aw = 0 − 0.05, water pulls together sulfonate groups that would otherwise be homogeneously dispersed in the polymer. 12 These clusters of sulfonates are thought to act as ionic cross-links that stiffen the membrane. 12 This non-monotonic, stiffening/plasticizing phenomenon has also been observed in DMA experiments. 17,18 Intriguingly, no mechanical measurement to date has detected this behavior in the Na-form of Nafion. 20 For both practical and fundamental reasons, the mechanical properties of Nafion at low hydration are of particular interest, yet they cannot be easily interpreted in terms of microscopic structure using X-ray or neutron scattering. An alternative approach is molecular simulation, since in principle it can predict the self-assembled morphology purely from atomistic forces. However, only a small fraction of the published simulation literature has attempted to calculate either the mechanical properties of Nafion or the structure of Nafion at low humidity. One challenge is how to classically represent the sulfonic hydrogen. In the limit of a single chain at infinite dilution, the sulfonic acid group can be safely represented as a sulfonate interacting non-covalently with a hydronium ion. The assumption of complete deprotonation breaks down at hydrations of less than three waters per sulfonic acid group (λ = 3), according to small-scale ab initio simulations of hydrated CF3 SO3 H. 24 In these simulations, the SO3 H

5



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

group existed in a protonated instead of deprotonated form at λ < 3. 24 Ab initio calculations have also been performed on full Nafion chains to study proton transport at a wide range of water contents (λ ≥ 3), albeit with limited length scales (< 3 nm) and time scales (< 130 ps). 25,26 In one recent study, Devanathan et al. observed that only 6.6% of protons were bonded to SO3– groups at λ = 3, confirming that at this level of hydration, full deprotonation is a safe assumption. 26 75% of protons were found to exist as SO3– /H3 O+ contact ion pairs. This fraction dropped to 29% when λ was increased to 9, because higher hydration favored H5 O2+ , H7 O3+ , and H9 O4+ ions. These species have also been observed in other confined environments like carbon nanotubes. 26 The delocalized excess proton in Nafion can be more cheaply simulated using reactive molecular dynamics methods, for example the multistate empirical valence bond approach. 27 Existing studies of Nafion using the multistate empirical valence bond approach have focused on proton conduction at moderate to high hydration (λ ≥ 6). However, these studies only considered charge transfer between water molecules, assuming complete deprotonation of sulfonic acid groups. 28,29 The accuracy of assuming full deprotonation was also assessed using classical molecular dynamics (MD) calculations of the water sorption isotherm. 30 These calculations assumed full deprotonation at all hydration levels. At moderate hydrations, the simulations agreed very well with experiments, but at λ < 2, the simulations began to systematically underpredict the experimental water activity. This was attributed to water molecules interacting too strongly with the membrane because they were represented as hydronium ions. Despite this limitation, some structural data were reported. Daly et al. observed rod-like hydrophilic clusters that persist to moderate humidity (λ ≈ 4); this highly aspherical morphology is consistent with analysis of experimental proton conductivities based on percolation theory. 22,30 Devanathan et al. studied clustering of water and hydronium oxygens at λ = 1 − 20, and estimated the percolation threshold of these clusters to be λ = 5 − 6. 31 Interestingly, the percolation transition based on analysis of experimental conductivity data occurs at λ ≈ 3.3

6


Page 6 of 37

Page 7 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(10 vol.% H2 O), 22 demonstrating the ambiguity in defining this transition. Another challenge to simulating Nafion at low humidity is the long relaxation time of the polymer chains. Lucid et al. found that even short chains with a MW of ∼35000 Da could not be equilibrated with conventional MD at a moderate hydration level (λ = 6.5) on the time scale of the simulations (50 ns). 32 The authors found that lower energy configurations could be achieved with temperature-accelerated MD, 33 although there is still no guarantee that these chains were at their theoretical equilibrium state. Equilibration can also be improved by using coarse-grained techniques such as dissipative particle dynamics. 34 When applied to Nafion, this technique increases the allowable timestep by up to a factor of 1000. 35 On the other hand, the resolution was sacrificed to the extent that a CF2 SO3 · (H2 O)3 cluster was represented as a single particle. This low resolution made the technique unsuitable for hydrations where λ < 3. Moreover, at this resolution, the structure of the counterions is totally invisible. Nevertheless, the technique proved useful for studying the effect of varying the length of the side chain. Side chains with 7 consecutive carbons/ether groups were observed to have larger water clusters than side chains with 4 consecutive carbons/ether groups. 35 Combining low-humidity simulations with calculations of mechanical properties demands not only long times for equilibration, but also long times for production (analysis/property calculation). Long production times are necessary because even simple mechanical properties such as the shear viscosity are challenging to compute with low uncertainty. 36 Additionally, mechanical properties like elastic moduli and elongational creep are strong functions of the MW and crystallinity of the polymer backbone, parameters that correlate strongly with equilibration time. Therefore, measurements on experimental systems cannot be quantitatively predicted using molecular simulation unless an equilibrated morphology is available through some other technique. A system of short Nafion chains with a MW of ∼15000 Da was used by Ozmaian and Naghdabadi to estimate the Tg at a hydration of λ = 14. 37 To obtain the Tg , the system was cooled with NPT molecular dynamics at rates of 0.125-0.5

7



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

K/ps. 37 These cooling rates are roughly 13-14 orders of magnitude faster than experimental cooling rates, so mode-coupling theory was used to extrapolate to experimental time scales and compare glass transition temperatures. The observed Tg of 345 K falls in between the experimental temperatures for the β transition (∼293 K) and the α transition (∼383 K) in dry membranes. 13 In experiments, if water is added to dry membranes, the α transition shifts to higher temperatures 17,18 and the β transition shifts to lower temperatures. 38 Hence, this Tg from simulations does not currently have an obvious interpretation. To the best of our knowledge, no other attempts have been made to calculate mechanical properties from molecular simulations of Nafion. The objective of the present work is to better understand the correlation between experimentally observed trends in mechanical properties and the underlying behavior of the SO3– /Xv+ /SO3– cross-links. Properties of particular interest include the increase in stiffness as the cation is exchanged from H+ to Na+ to Ca2+ . Also of interest is the non-monotonic dependence of stiffness on humidity at T & 323 K with H+ as the cation. The strength and connectivity of the ionic cross-links are quantitatively measured using MD simulations with semi-atomistic resolution. These structural properties are compared to calculations of the shear viscosity in the limit of zero shear rate, a substitute for less accessible mechanical properties like the elastic and storage moduli. To make the simulations tractable, the MW of the chains is limited to 1181-3468 Da and crystalline polymer domains are not explicitly included in the initial configurations. Equilibration of the systems is accelerated using replica-exchange molecular dynamics (REMD). 39 Equilibrated configurations are then provided as input to NVT-MD simulations to calculate transport properties. Further details of the simulation approach are explained in the Methods section. The following section, Results and Discussion, reports calculations of the shear viscosity and the translational relaxation time of the chains. Also included are snapshots of the SO3– /Xv+ /SO3– cross-links as well as various quantitative measures of the strength and connectivity of the cross-links. Finally, the Conclusions section summarizes the observed trends in transport properties and their

8


Page 8 of 37

Page 9 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


molecular-level explanation.

2

Methods

2.1

Force field

Simulated systems consist of three components: polymer electrolyte chains (Nafion), solvent molecules (water), and counterions (H+ , Na+ , or Ca2+ ). Chains of two different molecular weights are considered: 1181 and 3468 Da. In general, the chain architecture can be partly described by the EW, which in this study is 1181 g/mol. The exact locations of ionizable side-chains along the polymer backbone is unknown in experimental chains, 40 so the sidechains are separated from one another by 15 carbons, and separated from the termini by 7-8 carbons, making the overall arrangement approximately symmetric. The chains are represented semi-atomistically, with united-atom particles for CFX groups and explicit particles for all remaining atoms. Particles can interact through electrostatic, Lennard-Jones, bondstretching, bond-bending, and torsional forces. Lennard-Jones interactions were shifted and truncated at 1.2 nm, with standard long-range corrections added to both the pressure and potential energy. Electrostatic interactions were computed using the Smooth Particle Mesh Ewald (SPME) method 41 with fourth-order interpolation, a short-ranged cut-off of 1.2 nm, a κ of 2.6 nm−1 , and a grid spacing of no more than 0.1 nm. The functional forms and parameters of these interactions, unless otherwise specified, can be found in Daly et al. 30 , which modified the force-field in Cui et al. 42 . This particular parameterization has already been used to study water sorption and morphology, 30 as well as bulk 30 and interfacial transport of water. 43 Chains in the (RSO3 )2 Ca, RSO3 Na, RSO3 H · (H2 O)3 , and RSO3 H · (H2 O)10 systems are not covalently bound to their counterions. In the RSO3 H · (H2 O)3 and RSO3 H · (H2 O)10 systems, the H+ cation exists as H3 O+ , which is a good approximation at these hydration levels, according to ab initio studies. 24 Moreover, this approximation, together with the parameters from Daly et al. 30 , is sufficiently accurate 9



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

to reproduce the water sorption isotherm from experiments at intermediate hydration levels (λ > 2). 30 In the (RSO3 )2 Ca and RSO3 Na systems, cations are described by parameters from Aaqvist 44 . Water molecules in all systems are modeled using TIP4P/2005. 45 Sulfonic acid groups in the RSO3 H (1181 Da) and (RSO3 H)3 (3468 Da) systems are assumed to be protonated, with covalently bound H+ atoms. The necessary, additional bond-stretching, bond-bending, and torsional parameters for the sulfonic acid group are taken from a recently developed force field for protonated triflic acid. 46 Charges for the entire side-chain are re-calculated using the Gasteiger-Marsili method; 47 electronegativity parameters not available in Gasteiger and Marsili 47 were taken from the Vega ZZ software package. 48,49 The Gasteiger-Marsili method had also been used to calculate charges in the original force-field for deprotonated chains. 30,42,50 The new charges for protonated chains are available in the Supporting Information. In general, partial charges on acidic hydrogens that are calculated using the Gasteiger-Marsili method exhibit a roughly linear correlation with pKa . 47 This trend is closely adhered to by the sulfonic hydrogen in Nafion, as shown in the Supporting Information.

2.2

NPT-REMD, NPT-MD, and NVT-MD

Unless otherwise specified, systems were equilibrated using replica-exchange molecular dynamics 39 in the isobaric-isothermal (NPT) ensemble 51 at a pressure of 1 bar. This technique consists of running multiple molecular dynamics simulations, or “replicas”, in parallel over a wide temperature range, and periodically attempting Monte Carlo moves that exchange configurations between replicas. Attempts are accepted with a probability that preserves the statistical-mechanical ensemble, in this case the NPT ensemble. 51 This technique allows configurations that would be kinetically trapped at low temperatures to anneal at higher temperatures, overcoming free energy barriers. The REMD implementation was GROMACS 4.6.3/4.6.5. 52 Replicas were kept at constant temperatures and pressures using the Nosé-Hoover ther10


Page 10 of 37

Page 11 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


mostat 53 and Parinello-Rahman barostat, 54 respectively. In a given REMD run, each replica i was assigned a temperature Ti such that Ti = T0 k i , where k is a constant. This distribution in temperature space is thought to keep the acceptance probability approximately uniform for exchange attempts. 55 A typical run would consist of 57 replicas spanning a temperature range starting at 200 K and ending at 500 K, so that the smallest temperature interval would be ∼3 K, and the largest interval would be ∼10 K. Exchanges between replicas were attempted with frequencies ranging from 10 to 1 ps–1 , which have been shown to work well for aqueous systems. 56 Between exchanges, MD was performed with a time step of 1-2 fs, with smaller time steps used for systems with H3 O+ because of the stiff bond-stretching potential in that molecule. Systems were equilibrated for 10-40 ns, with longer times used for slowly relaxing systems like (RSO3 )2 Ca. On this time scale, configurations make multiple round-trips through temperature space, as shown in Figure S3 in the Supporting Information. At very high temperatures (T > 500 K), most systems could reach equilibrium without the enhanced sampling from REMD. Therefore, to conserve computational resources, NPT-MD was carried out instead of NPT-REMD for the RSO3 Na and RSO3 H · (H2 O)3 systems at T > 500 K. After equilibration, all simulations were continued for another 10 ns to generate configurations for calculating structural properties. From these 10 ns production runs, initial configurations were selected for subsequent NVT-MD runs. Disabling the barostat and replica exchanges were necessary for calculating the viscosity and translational relaxation time (see Section 2.3). These transport properties are derived from time-dependent correlations that decay over very long time scales. Therefore, to obtain adequately low uncertainties for these properties, systems were simulated for 0.2-2.2 μs. On the other hand, the REMD production runs were far shorter because they were used for purely static properties. Moreover, REMD can generally explore configuration space more efficiently than MD, further reducing the need for long production runs. In the NVT-MD production runs, the temperature was kept fixed with the Nosé-Hoover

11



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 37

thermostat, 53 as in the NPT-REMD and NPT-MD simulations.

2.3

Analysis

The mechanical stiffness of the polymer electrolyte was quantified by computing the shear viscosity in the limit of zero shear rate. This mechanical property is a natural choice for the low-MW chains considered in this work, and is also a substitute for mechanical properties like the elastic and storage moduli that are less accessible to simulations. The underlying physics for all of these mechanical properties are very similar. For example, the following relationships hold for many polymer systems: 57

η(γ=ω) ˙ ≈

1p ′ [G (ω)]2 + [G′′ (ω)]2 ω

η¯(ǫ→0) ˙ ≈ 3η(γ→0) ˙

where η(γ=ω) ˙ is the shear viscosity at a shear rate of ω. G′ (ω) and G′′ (ω) are the storage and loss moduli, respectively, at a frequency of ω. η¯(ǫ) ˙ is the elongational viscosity at an elongation rate of ǫ. ˙ Although these relationships may not rigorously hold for Nafion, they suggest that we should expect at least qualitatively similar behavior among the various mechanical properties. The viscosity at zero shear rate can be denoted as η(γ=0), ˙ or more compactly, η0 . The viscosity was calculated from the pressure tensor using a Green-Kubo (GK) relation originally derived by Daivis and Evans: 58,59

η0

V = 10kB T

Z

t=∞

dt t=0

X

hPαβ (to )Pαβ (to + t)i

(1)

!

(2)

αβ

Pαβ = (παβ + πβα )/2 − δβα

X γ

παβ

πγγ

/3

" N # N N 1 XX 1 X mi viα viα + rijα fijβ = V i=1 2 i=1 j=1 12


(3)

Page 13 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where mi is the mass of particle i and viα is the α component of the velocity of particle i, and α is x, y, or z. rijα is the distance, along the α axis, between particle i and the minimum image of particle j with respect to particle i. fijβ is β component of the force acting between particles i and j, and includes both bonded and non-bonded interactions. δαβ is the Kronecker delta. The particular form of the right hand side of eq 3 is known as the “atomic virial form”. 60 The advantage of the GK relation in eq 1 over the more conventional version 61 is that the conventional version only incorporates off-diagonal components of the pressure tensor, and therefore utilizes 40% less data available in the pressure tensor. The correlation function hPαβ (to )Pαβ (to + t)i was computed every 10-50 fs, with shorter intervals used at higher temperatures where the correlation function decays faster. The time origin, to , can be defined at any point in the NVT-MD simulation, and multiple time origins were considered at intervals of 1-10 ps to increase sampling. The correlation function was integrated using the trapezoidal rule and inserted into eq 1 to estimate η0 . The integral cannot be carried out to t = ∞ for simulations of finite length, so the integration was truncated when a plateau in η0 (tf ) was observed, where η0 (tf ) is the expression on the righthand side of eq 1 integrated to t = tf instead of t = ∞. An example of η0 (tf ) over a range in tf can be found in Figure S4 in the Supporting Information. For slowly relaxing systems with very large values of η0 , the plateau was sometimes not easily detectable by visual inspection, so η0 (tf ) was fit to a saturation-growth-rate equation and extrapolated to t = ∞:

η0 (tf ) ≈ η0 (∞)

tf C + tf

(4)

η0 (∞) and C are fitting parameters that can be obtained by regressing 1/η0 (tf ) on 1/tf using a straight line. In practice, this approximation resulted in less truncation error at the expense of large uncertainties in η0 (∞) due to the extrapolation (for example see Figure S5 in the Supporting Information). The mechanical stiffness of the polymer electrolyte was also quantified by computing the

13



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 37

translational relaxation time, defined as follows:

τ := t such that hr2 (t)i=(1.5Rg )2

(5)

hr2 (t)i is the mean-squared displacement of the center of mass of a chain as a function of time, and Rg is the radius of gyration of the chain. The definition can be interpreted as the time it takes for a chain to diffuse a distance on the order of its own size. This type of relaxation time is commonly used in polymer science. 23 Figure S7 in the Supporting Information demonstrates how the translational relaxation time was calculated. A structural property relevant to mechanical stiffness is the number of sulfur atoms crosslinked together in a cluster or network. Two sulfur atoms are defined as cross-linked if there is at least one cation in the system that is less than a prescribed distance away from both sulfurs. This definition is illustrated in Figure 1. Two sulfurs belong to the same network if they are either directly cross-linked, or can be connected by a path of cross-links. A system can consist of multiple networks of cross-linked sulfurs, and the size of these networks is measured by computing a weighted arithmetic mean, MS , of the network size, N :

MS =

cN N 2 N P cN N

P

(6)

N

cN is the count of networks of size N , and was computed using the hierarchical clustering code in Scipy. 62 This quantity is analogous to the weight-average molecular weight of polydisperse polymer melts, 23 and can be interpreted as the average size obtained by picking sulfur atoms at random and measuring the sizes of the networks to which they belong.

14


Page 15 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1: Schematic of a hydrophilic cluster in a Nafion polymer matrix with an arbitrary cation type. Water molecules are represented with stick models, while other hydrophilic species are represented with space-filling models. Two sulfur atoms are considered crosslinked if r1 and r2 , the distances to an arbitrary cation in the system, are both less than rc (defined in Figure 5). Sulfur atoms are colored according to the cross-link network to which they belong.

3 3.1

Results Systems near equilibrium

The polymer electrolytes in this study approach equilibrium very slowly, despite their low MW. Their equilibration was substantially accelerated with NPT-REMD, especially at low temperatures, as shown in Figure 2. For example, at 200 K, chains in the RSO3 H system settled into a potential energy minimum that was roughly 10kB T lower than that accessible to NPT-MD. On the other hand, chains at 500 K could traverse phase space much more easily, so as expected, the chains reached the same potential energy, regardless of whether NPTMD or NPT-REMD was used. This example demonstrates that at low relative humidity, enhanced sampling techniques such as REMD are essential for studying Nafion with atomistic resolution at all but the highest temperatures.

15



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 37

Figure 2: Total energy per chain, in kB T , for the RSO3 H system as a function of equilibration time using REMD vs. MD.

The viscosities of equilibrated systems are shown in Figure 3. When viewed on an Arrhenius plot over a ∼800 K temperature range, these viscosities reveal clear trends in mechanical stiffness, despite the large uncertainties that are inevitable for these slowly relaxing systems. The uncertainties were estimated by first block averaging data that are adjacent in inverse temperature space: N 1 X log10 [η0 (1/T )] = log10 η0i (1/T i ) N i=1

1 T

N 1 X 1 = N i=1 T i

(7) (8)

The upper and lower error bounds on log10 [η0 (1/T )] were estimated, respectively, as the maximum and minimum values in the block of N adjacent data. As an added benefit, the blocking operation also smoothed the data. In Figure 3, log10 [η0 (1/T )] with N = 3 is shown, and the raw, unblocked data can be found in Figure S6 in the Supporting Information.

16


Page 17 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 3: Shear viscosity in the limit of zero shear rate as a function of inverse temperature. Different symbols denote systems with varying MW, hydration levels and cation types. Purple N: (RSO3 )2 Ca. Cyan : RSO3 Na. Black H: (RSO3 H)3 . Green : RSO3 H · (H2 O)3 . Red : RSO3 H · (H2 O)10 . Blue ⋆: RSO3 H. Remarkably, the data in Figure 3 span 4-5 decades at a single temperature, despite very modest changes in composition. For example, at T ≈ 580 K, η0 for (RSO3 )2 Ca and RSO3 H · (H2 O)10 is O(10 Pa·s) and O(10−4 Pa·s), respectively. At T ≈ 500 K, η0 differs by nearly four orders of magnitude between the RSO3 Na and RSO3 H systems, even though Na+ represents ∼2 wt% of the RSO3 Na system. This ratio of viscosities is approximately equal to the ratio of the viscosities of water and honey at room temperature. 63 The observed trend in viscosity, (RSO3 )2 Ca > RSO3 Na > RSO3 H, agrees qualitatively with experimental measurements of the storage modulus (see Figure S9 in the Supporting Information). When the cation is restricted to H+ and the MW to 1181 Da, Figure 3 shows a clear non-monotonic dependence of viscosity on humidity: RSO3 H · (H2 O)3 > RSO3 H ≈ RSO3 H · (H2 O)10 . The corresponding values of aw are approximately 0.3, 0, and 0.8, respectively, using previous simulation data obtained with the same force field. 30 This trend in viscosity qualitatively agrees with experimental measurements of the elongational creep at T = 323 K (see Figure S8 in the Supporting Information). 12 The viscosity is inversely proportional to the self-diffusivity for pure-component, smallmolecule Newtonian fluids. 64 This relationship breaks down for very high MW materials, 64 17



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

but it is unclear whether this criterion is met by concentrated systems of 1181-Da Nafion chains. The self-diffusivity of even short Nafion chains is not accessible to atomistic simulations, but a relaxation time is often used in lieu of a self-diffusivity, in which case τ ∼ η is sometimes assumed. 65 The definition of a relaxation time is ambiguous, 59 so a common definition from polymer physics was used in this study, as explained in Section 2.3. Clearly, τ exhibits the same trends with respect to aw and cation type as η0 , as shown in Figure 4. For RSO3 Na systems, η0 has an activation energy of 9100 ± 700 K over the available temperature range, while τ has an activation energy of 11500 ± 600 K. For RSO3 H · (H2 O)3 systems, the activation energies for η0 and τ are 4400 ± 300 and 4400 ± 100 K, respectively. The similar activation energies for the two properties suggest that τ ∼ η0 is a reasonable approximation for these systems. Importantly, η0 and τ are obtained from two very different types of simulation data: the pressure tensor and the mean-squared displacement of the chains, respectively. Across these two properties, the trends with respect to aw and cation type are consistent, suggesting that the trends are not very sensitive to the measure used to quantify mechanical stiffness.

Figure 4: Translational relaxation time (eq 5) as a function of inverse temperature. Different symbols denote systems with varying hydration levels and cation types. Purple N: (RSO3 )2 Ca. Cyan : RSO3 Na. Green : RSO3 H · (H2 O)3 . Red : RSO3 H · (H2 O)10 . Blue ⋆: RSO3 H. One objective of the present work is to better understand the correlation between η0 or τ 18


Page 18 of 37

Page 19 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and the behavior of the SO3– /Xv+ /SO3– cross-links. The strength of SO3– /Xv+ interactions has been shown by experimental measurements of conductivity 66 and IR spectra 13 to increase significantly when Xv+ is changed from Na+ to Ca2+ . In simulations, the strength of these interactions can be quantified by computing the Xv+ /SO3– pair correlation function, g(r). This function can be expressed in terms of the potential of mean force, w(r): 67

w(r) g(r) = exp − kB T

The potential of mean force can be interpreted as the change in free energy when two particles in a solvent (the system) are moved from infinite separation to a distance r. 67 Minima in w(r) correspond to maxima in g(r). The function g(r) at T = 353 K is shown in Figure 5. The g(r)’s of all systems include a broad, well-defined peak between 0.325-0.4 nm. The g(r)’s of the RSO3 H and (RSO3 H)3 systems have an additional peak at 0.215 nm that represents the intra-molecular sulfur that is covalently connected to the hydrogen through an oxygen (S-O-H). This intra-molecular peak will be present regardless of whether SO3– /Xv+ /SO3– cross-links exist in the system, and therefore does not reflect the strength of the cross-links. Instead, we define the strength of the cross-link as the maximum of g(r) (i.e. minima in w(r)) subject to the constraint that r > r∗ , where r∗ is set to 0.25 nm to exclude the intra-molecular peak.

19



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 37

Figure 5: Cation-sulfur pair correlation function at T = 353 K. Purple N: (RSO3 )2 Ca. Cyan : RSO3 Na. Black H: (RSO3 H)3 . Green : RSO3 H · (H2 O)3 . Red : RSO3 H · (H2 O)10 . Blue ⋆: RSO3 H. Black dotted line: distance (r∗ ) below which peaks are interpreted as intramolecular. Black dashed line: distance (rc ) below which sulfurs are interpreted as belonging to the first coordination shell.

The maxima of g(r) are shown in Figure 6, and the following trend is observed: (RSO3 )2 Ca > RSO3 Na > RSO3 H · (H2 O)3 > RSO3 H · (H2 O)10 > RSO3 H. Ca2+ interacts more strongly with SO3– than Na+ , as expected from experimental measurements of conductivity 66 and IR spectra. 13 Both of these cations form stronger cross-links than H+ . Adding water to the RSO3 H system appears to strengthen the cross-links, probably because the dipole-dipole interactions between SO3 H groups get replaced by stronger charge-charge interactions between SO3– and H3 O+ . Increasing the water content from RSO3 H · (H2 O)3 to RSO3 H · (H2 O)10 weakens the cross-links, although the position of the peak in Figure 5 does not change, suggesting that individual pairs of SO3– and Xv+ are left intact. Adjacent SO3– /Xv+ pairs become separated by the additional water, and this increased separation, in combination with dielectric screening from water, weakens the charge-charge interactions.

20


Page 21 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6: The maximum of g(r) in Figure 5 subject to the constraint that r > r∗ . Purple N: (RSO3 )2 Ca. Cyan : RSO3 Na. Black H: (RSO3 H)3 . Green : RSO3 H · (H2 O)3 . Red : RSO3 H · (H2 O)10 . Blue ⋆: RSO3 H. This dilution effect of the solvent can also be measured by calculating the number of sulfurs in the first coordination shell of the cation. The first coordination shell corresponds to the first peak in g(r), 67 so the occupancy of the shell can be obtained by integrating under that peak and multiplying by the density. As in Figure 6, we exclude intra-molecular peaks (r < r∗ ). We truncate the integration at rc = 0.45 nm, since this cut-off is sufficiently large to include the first intermolecular peak in all systems, as shown in Figure 5. Figure 7 indicates that there is roughly 1 sulfur in the first coordination shell of H3 O+ in the RSO3 H · (H2 O)10 system, confirming that most SO3– /Xv+ pairs are isolated from one another. As expected, the coordination shell of Ca2+ contains the most sulfurs, since it has a higher valence than the cations in the other systems. Interestingly, both the maxima in g(r) and the coordination number are relatively insensitive to temperature.

21



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 7: Number of intermolecular sulfurs in the first coordination shell of the cation, where the first coordination shell is defined as r < rc (see Figure 5). Purple N: (RSO3 )2 Ca. Cyan : RSO3 Na. Black H: (RSO3 H)3 . Green : RSO3 H · (H2 O)3 . Red : RSO3 H · (H2 O)10 . Blue ⋆: RSO3 H. The coordination number is a measure of nearest-neighbor connectivity; the long-range connectivity of SO3– /Xv+ /SO3– cross-links may also influence η0 and τ . Distinct networks of cross-links are shown schematically in Figure 1, and in actual configurations at T = 353 K in Figure 8. Most strikingly, the networks in the RSO3 Na system are completely percolated, despite the total absence of water molecules that are needed for RSO3 H systems to percolate. Surprisingly, the (RSO3 )2 Ca has smaller networks than the RSO3 Na system, despite having much stronger bonds and more nearest neighbors (see Figures 6 and 7). We speculate that the smaller networks in the (RSO3 )2 Ca system are a result of the electroneutrality constraint: there are half as many Ca2+ ions as Na+ ions. The greater number of Na+ ions may allow the SO3– groups to assemble into string-like clusters that percolate more easily, as opposed to the more spherical clusters seen in (RSO3 )2 Ca (see Figure 8).

22


Page 22 of 37

Page 23 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8: Networks of cross-linked sulfur atoms (defined in Figure 1) at an arbitrary point in time for systems at T = 353 K. Opaque colored blobs are Connolly surfaces 68 enclosing sulfurs belonging to distinct networks. Each color represents a separate network. Translucent gray blobs are Connolly surfaces enclosing all hydrophilic species (SO3– , H2 O, cations). Figure 8 illustrates that percolation of the hydrophilic clusters does not necessarily imply percolation of the SO3– /Xv+ /SO3– cross-links, even in the RSO3 H·(H2 O)3 system. Nevertheless, increasing the water content from RSO3 H to RSO3 H · (H2 O)3 results in enlarged crosslinked networks. On the other hand, increasing the water content further to RSO3 H·(H2 O)10 greatly shrinks the networks by diluting the SO3– /Xv+ pairs. Figure 8 shows a clear non-monotonic dependence of network size on water content; this dependence is quantified in Figure 9. Notably, the RSO3 H system never achieves a truly disordered state (MS ≈ 1) across the temperature range of 200-500 K, suggesting that

23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

increasing the temperature merely shifts the reaction equilibrium for clustering, rather than driving a sharp order/disorder transition. 12 On the other hand, the expected transition to monotonic behavior at T < 323 K 12 is not observed in η0 or τ , leaving open the possibility that by missing the sharp order/disorder structural transition, the simulations also miss the mechanical transition. More simulations with much higher MW chains are needed to check whether the mechanical transition involves a change in the crystalline regions or the conformations of the polymer backbone.

Figure 9: Weighted arithmetic mean (eq 6) of the size of the cross-linked networks (defined in Figure 1) as a function of temperature. Purple N: (RSO3 )2 Ca. Cyan : RSO3 Na. Black H: (RSO3 H)3 . Green : RSO3 H · (H2 O)3 . Red : RSO3 H · (H2 O)10 . Blue ⋆: RSO3 H. Black dashed line: the total number of sulfur atoms in the system. Blue dotted line: approximate location of the percolation transition in the RSO3 Na system. Another notable feature of Figure 9 is a percolation transition in the RSO3 Na system at T ≈ 620 K, a remarkably high temperature. The extremely high connectivity of the RSO3 Na system would explain why the corresponding values of η0 and τ are not substantially lower than those of (RSO3 )2 Ca, despite the fact that Ca2+ interacts much more strongly with SO3– groups, as shown in Figure 6. Similarly, in Figure 6, SO3– /Na+ interactions are only marginally stronger than SO3– /H3 O+ interactions in RSO3 H · (H2 O)3 , yet large differences between the two systems are observed in η0 and τ . The high connectivity of RSO3 Na appears to greatly increase the mechanical stiffness; this hypothesis will be discussed further

24


Page 24 of 37

Page 25 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


in Section 3.2. The RSO3 Na, (RSO3 )2 Ca, and RSO3 H systems have networks that span 1-2 orders of magnitude in size at a given temperature, demonstrating that morphology is a strong function of counterion type at low humidity. Counterintuitively, morphology has been suggested to be insensitive to counterion type at moderate to high humidity (aw & 0.3), a hypothesis based on the similarity of the small-angle X-ray scattering spectra (SAXS) for different cation-exchanged forms of Nafion. 69 At moderate to high humidity, there is enough contrast in electron density to detect clusters, whereas under dry conditions, the clusters may be invisible, depending on counterion type. 14 For example, clusters normally show up as an “ionomer peak” at a Bragg spacing of 3-5 nm in the SAXS spectrum; 13 no such peak was observed by Fujimura et al. in the Na-form of Nafion at aw = 0. 14 Fujimura et al. speculated that clusters exist in the Na-form, but are simply not detectable due to insufficient contrast in electron density. 14 This hypothesis is supported by the the larger clusters observed in the RSO3 Na system in Figure 9. The results in Figure 9 also imply that other counterions like Cs+ and N(CH3 )4+ cannot be assumed to have a negligible impact on the morphology at low humidities. These counterions are of practical interest because they enhance contrast in scattering and TEM experiments. 13–16 Cs+ has also been explored as a dopant that reduces methanol crossover in direct methanol fuel cells. 70 The size of cross-linked networks, MS , appears to be necessary but not sufficient for explaining the trends in η0 and τ . For example, the RSO3 H system has cross-link networks that are much larger than those of the RSO3 H · (H2 O)10 system, yet for the two systems, η0 and τ are nearly the same. The somewhat stronger SO3– /H3 O+ interactions in the RSO3 H · (H2 O)10 system, as observed in Figure 6, appear to compensate for the smaller networks of cross-links. Hence, it is plausible that both the strength of the cross-links and the size of the cross-link networks determine the mechanical stiffness. The dependence of the mechanical and structural data on the MW of the chains was examined by simulating a system of 3468 Da chains, denoted as (RSO3 H)3 . An increase

25



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

in MW stiffens the polymer electrolyte to a surprisingly large degree, as seen in Figure 3. To quantify this stiffening, we fitted each data set to an Arrhenius equation, η0 = C exp(−Ea /kB T ). From this fitting, we find that η0 ((RSO3 H)3 )/η0 (RSO3 H) ≈ 46 over 1/T = 2 − 3.5. The viscosity would be expected to increase by factors of 3 and 42 for simple polymer melts containing unentangled and entangled chains, respectively, assuming a constant monomer friction coefficient. 23 The structural differences between the (RSO3 H)3 and RSO3 H systems are minimal according to Figures 6,7, and 9. On the other hand, the density increase in Figure 10 is roughly 5-10%, a non-negligible amount. The dependence of viscosity on density is discussed in Section 3.2.

Figure 10: Mass density at p = 1 bar as a function of temperature. Purple N: (RSO3 )2 Ca. Cyan : RSO3 Na. Black H: (RSO3 H)3 . Green : RSO3 H·(H2 O)3 . Red : RSO3 H·(H2 O)10 . Blue ⋆: RSO3 H.

3.2

Systems away from equilibrium

Varying the composition of Nafion changes many structural properties simultaneously, as Figures 5-10 illustrate. In an effort to isolate the individual contributions of these properties, certain quantities can be fixed while varying others by taking systems out of equilibrium. This process was achieved by rapidly cooling and/or compressing equilibrated systems, and then performing NVT-MD as described in Section 2.2. The details of these out-of-equilibrium systems are summarized in 1. In the first run, an equilibrated RSO3 Na system at T = 687 26


Page 26 of 37

Page 27 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


K was simultaneously cooled and compressed over the course of 7 ps to the approximate temperature and volume of an equilibrated reference system at T = 469 K. The structural properties gmax and N (rc ) are nearly the same between the quenched system and equilibrated reference system at T = 469 K, suggesting that these properties are relatively insensitive to temperature, as shown in Figures 6 and 7. On the other hand, MS is 87% smaller in the quenched system, indicating that the networks of cross-links in the quenched system are largely unchanged from the initial configuration at T = 687 K. Clearly, the networks cannot rearrange themselves and percolate on the time scale of the simulation without an enhanced sampling technique like REMD. The translational diffusion time is 98% smaller in the quenched system, a decrease in stiffness that can be attributed to the absence of long-range connectivity in the quenched system. Table 1: Details of systems taken out of equilibrium by rapid cooling and/or compression.a System

Ti

Tf

ρi

ρf

T f /T ref

ρf /ρref

f ref gmax /gmax

N (rc )f /N (rc )ref

MSf /MSref

τ f /τ ref

RNa RH

687 500

469 500

1.70 1.87

2.00 2.07

1.00 1.00

0.99 1.10

1.04 0.99

1.08 1.10

0.13 1.20

0.02 2.73

a

Superscripts “i” and “f” denote properties of the system before and after quenching, respectively. Superscript “ref” denotes properties of the reference system at equilibrium, included for comparison.

In the second run, an equilibrated RSO3 H system at T = 500 K was compressed over the course of 2 ps to the approximate volume of the (RSO3 H)3 system at T = 500 K. Upon reaching this higher volume, the system was relaxed for an additional 5 ns before analysis of its structural and transport properties. The structural properties N (rc ) and MS in 1 are slightly larger than those of the RSO3 H reference system, but approximately equal to those of the (RSO3 H)3 system, according to Figures 7 and 9. τ increased by 273% relative to the equilibrated RSO3 H system, suggesting that under equilibrium conditions, the higher density of the (RSO3 H)3 system is partly responsible for the factor of 46 increase in η0 discussed in Section 3.1. Figures 3 and 4 have already shown that τ ∼ η0 is a reasonable approximation,

27



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

so we can estimate the isochoric ratio of viscosities for the RSO3 H and (RSO3 H)3 systems by dividing the factor of 46 by 2.73, obtaining 17. This factor is still larger than the factor of 3 we would expect for a system of unentangled chains with a constant monomer friction coefficient. 23 On the other hand, the assumption of a constant monomer friction coefficient is known to be less accurate for shorter chains. 23 In any case, the change in viscosity is significantly smaller than what we would expect for an entangled polymer melt (42), which strongly suggests that the chains in this study are below their critical entanglement length.

4

Conclusions

The viscosity and structure of Nafion are studied as a function of humidity and cation type using molecular dynamics simulations. Nafion is modeled as monodisperse, low-MW chains with no crystalline domains. This highly simplified representation succeeds in reproducing two important trends seen in experiments with much higher-MW chains: (1) at T & 323 K, water acts as an antiplasticizer at low humidity and a plasticizer at high humidity; (2) exchanging the H+ counterion with Na+ , or Na+ with Ca2+ , stiffens the membrane. This stiffening is particularly dramatic at T ≈ 580 K, where the viscosities of the RSO3 H·(H2 O)10 and (RSO3 )2 Ca systems differ by five orders of magnitude. In principle, the viscosities obtained in this study could be directly compared with experiments, provided that Nafion chains can be synthesized with a sufficiently low MW. The trends in viscosity are correlated with changes in the strength of the SO3– /Xv+ /SO3– cross-links and the size of the cross-link networks. Adding water to the RSO3 H system strengthens cross-links by abstracting protons to form H3 O+ ions, which are linked to surrounding SO3– groups through charge-charge interactions instead of weaker dipole-dipole interactions. In the RSO3 H · (H2 O)10 system, these stronger interactions compensate for the extremely disconnected networks of cross-links, making the system about as stiff as the RSO3 H system. Systems with Na+ form slightly stronger cross-links than those with H3 O+ ,

28


Page 28 of 37

Page 29 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


yet are significantly stiffer due to the complete percolation of their cross-links at T < 620 K. Surprisingly, the connectivity of the Na+ system is higher than even that of the Ca2+ system, despite the fact that the Ca2+ system has stronger cross-links and more sulfurs in its coordination shell due to its higher valence. Across all systems at aw = 0, the size of the cross-link networks span 1-2 orders of magnitude, demonstrating that morphology is very sensitive to cation type at low humidities. Therefore, caution should be used when drawing conclusions about the H-form of Nafion based on scattering and TEM experiments performed with contrast-enhancing agents such as Cs+ and N(CH3 )4+ .

Acknowledgements P. G. D. gratefully acknowledges the support of the National Science Foundation (Grant No. CHE-1213343). A. Z. P. would like to acknowledge support for this work from the Department of Energy, Office of Basic Energy Sciences, under grant DE-SC0002128.

Supporting Information Available: Partial charges for the RSO3 H and (RSO3 H)3 systems; details of the REMD simulations; details of the viscosity and relaxation time calculations; experimental measurements of mechanical properties. This material is available free of charge via the Internet at http://pubs.acs.org.

References (1) Hallinan Jr, D. T.; Balsara, N. P. Polymer Electrolytes. Annu. Rev. Mater. Res. 2013, 43, 503–525. (2) Costamagna, P.; Srinivasan, S. Quantum Jumps in the PEMFC Science and Technology from the 1960s to the Year 2000: Part II. Engineering, Technology Development and Application Aspects. J. Power Sources 2001, 102, 253–269.

29



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(3) Satterfield, M. B.; Majsztrik, P. W.; Ota, H.; Benziger, J. B.; Bocarsly, A. B. Mechanical Properties of Nafion and Titania/Nafion Composite Membranes for Polymer Electrolyte Membrane Fuel Cells. J. Polym. Sci. B Polym. Phys. 2006, 44, 2327–2345. (4) Collier, A.; Wang, H.; Zi Yuan, X.; Zhang, J.; Wilkinson, D. P. Degradation of Polymer Electrolyte Membranes. Int. J. Hydrogen Energy 2006, 31, 1838–1854. (5) Borup, R.; Meyers, J.; Pivovar, B.; Kim, Y. S.; Mukundan, R.; Garland, N.; Myers, D.; Wilson, M.; Garzon, F.; Wood, D. Scientific Aspects of Polymer Electrolyte Fuel Cell Durability and Degradation. Chem. Rev. 2007, 107, 3904–3951. (6) Jalani, N. H.; Datta, R. The Effect of Equivalent Weight, Temperature, Cationic R Forms, Sorbates, and Nanoinorganic Additives on the Sorption Behavior of Nafion .

J. Membr. Sci. 2005, 264, 167–175. (7) Lakshmanan, B.; Huang, W.; Olmeijer, D.; Weidner, J. W. Polyetheretherketone Membranes for Elevated Temperature PEMFCs. Electrochem. Solid-State Lett. 2003, 6, A282–A285. (8) Yang, C.; Costamagna, P.; Srinivasan, S.; Benziger, J.; Bocarsly, A. B. Approaches and Technical Challenges to High Temperature Operation of Proton Exchange Membrane Fuel Cells. J. Power Sources 2001, 103, 1–9. (9) Savadogo, O. Emerging Membranes for Electrochemical Systems: Part II. High Temperature Composite Membranes for Polymer Electrolyte Fuel Cell (PEFC) Applications. J. Power Sources 2004, 127, 135–161. (10) Zhang, J.; Xie, Z.; Zhang, J.; Tang, Y.; Song, C.; Navessin, T.; Shi, Z.; Song, D.; Wang, H.; Wilkinson, D. P. et al. High Temperature PEM Fuel Cells. J. Power Sources 2006, 160, 872–891.

30


Page 30 of 37

Page 31 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(11) Zhao, Q.; Benziger, J. Mechanical Properties of Perfluoro Sulfonated Acids: The Role of Temperature and Solute Activity. J. Polym. Sci. Part B: Polym. Phys. 2013, 51, 915–925. (12) Majsztrik, P. W.; Bocarsly, A. B.; Benziger, J. B. Viscoelastic Response of Nafion. Effects of Temperature and Hydration on Tensile Creep. Macromolecules 2008, 41, 9849–9862. (13) Mauritz, K. A.; Moore, R. B. State of Understanding of Nafion. Chem. Rev. 2004, 104, 4535–4586. (14) Fujimura, M.; Hashimoto, T.; Kawai, H. Small-Angle X-Ray Scattering Study of Perfluorinated Ionomer Membranes. 1. Origin of Two Scattering Maxima. Macromolecules 1981, 14, 1309–1315. (15) Fujimura, M.; Hashimoto, T.; Kawai, H. Small-Angle X-Ray Scattering Study of Perfluorinated Ionomer Membranes. 2. Models for Ionic Scattering Maximum. Macromolecules 1982, 15, 136–144. (16) Rubatat, L.; Rollet, A. L.; Gebel, G.; Diat, O. Evidence of Elongated Polymeric Aggregates in Nafion. Macromolecules 2002, 35, 4050–4055. (17) Kreuer, K.-D. The Role of Internal Pressure for the Hydration and Transport Properties of Ionomers and Polyelectrolytes. Solid State Ionics 2013, 252, 93–101. (18) Kreuer, K.-D. Ion Conducting Membranes for Fuel Cells and Other Electrochemical Devices. Chem. Mater. 2014, 26, 361–380. (19) Fan, Y.; Tongren, D.; Cornelius, C. J. The Role of a Metal Ion Within Nafion upon Its Physical and Gas Transport Properties. Eur. Polym. J. 2014, 50, 271–278. (20) Zhao, Q. Structures and Properties of Sulfonated Ionomers Probed by Transport and

31



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Mechanical Measurements: The Role of Solute Activity. Ph.D., Princeton University, United States – New Jersey, 2012. (21) Kusoglu, A.; Karlsson, A. M.; Santare, M. H.; Cleghorn, S.; Johnson, W. B. Mechanical Behavior of Fuel Cell Membranes Under Humidity Cycles and Effect of Swelling Anisotropy on the Fatigue Stresses. J. Power Sources 2007, 170, 345–358. (22) Wu, X.; Wang, X.; He, G.; Benziger, J. Differences in Water Sorption and Proton Conductivity Between Nafion and SPEEK. J. Polym. Sci., Part B: Polym. Phys. 2011, 49, 1437–1445. (23) Rubinstein, M.; Colby, R. H. Polymer Physics; Oxford University Press, USA, 2003. (24) Paddison, S. J. Proton Conduction in PEMs: Complexity, Cooperativity and Connectivity. Device and Materials Modeling in PEM Fuel Cells 2009, 385–412. (25) Choe, Y.-K.; Tsuchida, E.; Ikeshoji, T.; Yamakawa, S.; Hyodo, S.-a. Nature of Proton Dynamics in a Polymer Electrolyte Membrane, Nafion: A First-Principles Molecular Dynamics Study. Phys. Chem. Chem. Phys. 2009, 11, 3892. (26) Devanathan, R.; Idupulapati, N.; Baer, M. D.; Mundy, C. J.; Dupuis, M. Ab Initio Molecular Dynamics Simulation of Proton Hopping in a Model Polymer Membrane. J. Phys. Chem. B 2013, 117, 16522–16529. (27) Voth, G. A. Computer Simulation of Proton Solvation and Transport in Aqueous and Biomolecular Systems. Acc. Chem. Res. 2006, 39, 143–150. (28) Petersen, M. K.; Voth, G. A. Characterization of the Solvation and Transport of the Hydrated Proton in the Perfluorosulfonic Acid Membrane Nafion. J. Phys. Chem. B 2006, 110, 18594–18600. (29) Feng, S.; Voth, G. A. Proton Solvation and Transport in Hydrated Nafion. J. Phys. Chem. B 2011, 115, 5903–5912. 32


Page 32 of 37

Page 33 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(30) Daly, K. B.; Benziger, J. B.; Debenedetti, P. G.; Panagiotopoulos, A. Z. Molecular Dynamics Simulations of Water Sorption in a Perfluorosulfonic Acid Membrane. J. Phys. Chem. B 2013, 117, 12649–12660. (31) Devanathan, R.; Venkatnathan, A.; Rousseau, R.; Dupuis, M.; Frigato, T.; Gu, W.; Helms, V. Atomistic Simulation of Water Percolation and Proton Hopping in Nafion Fuel Cell Membrane. J. Phys. Chem. B 2010, 114, 13681–13690. (32) Lucid, J.; Meloni, S.; MacKernan, D.; Spohr, E.; Ciccotti, G. Probing the Structures of Hydrated Nafion in Different Morphologies Using Temperature-Accelerated Molecular Dynamics Simulations. J. Phys. Chem. C 2013, 117, 774–782. (33) Maragliano, L.; Vanden-Eijnden, E. A Temperature Accelerated Method for Sampling Free Energy and Determining Reaction Pathways in Rare Events Simulations. Chem. Phys. Lett. 2006, 426, 168–175. (34) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics. EPL 1992, 19, 155. (35) Wu, D.; Paddison, S. J.; Elliott, J. A. A Comparative Study of the Hydrated Morphologies of Perfluorosulfonic Acid Fuel Cell Membranes with Mesoscopic Simulations. Energy Environ. Sci. 2008, 1, 284–293. (36) Hess, B. Determining the Shear Viscosity of Model Liquids from Molecular Dynamics Simulations. J. Chem. Phys. 2002, 116, 209–217. (37) Ozmaian, M.; Naghdabadi, R. Molecular Dynamics Simulation Study of Glass Transition in Hydrated Nafion. J. Polym. Sci. Part B: Polym. Phys. 2014, 52, 907–915. (38) Yeo, S. C.; Eisenberg, A. Physical Properties and Supermolecular Structure of Perfluorinated Ion-Containing (nafion) Polymers. J. Appl. Polym. Sci. 1977, 21, 875–898.

33



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(39) Sugita, Y.; Okamoto, Y. Replica-Exchange Molecular Dynamics Method for Protein Folding. Chem. Phys. Lett. 1999, 314, 141–151. (40) Jang, S. S.; Molinero, V.; C ¸ agin, T.; Goddard III, W. A. Effect of Monomeric Sequence on Nanostructure and Water Dynamics in Nafion 117. Solid State Ionics 2004, 175, 805–808. (41) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577–8593. (42) Cui, S.; Liu, J.; Selvan, M. E.; Keffer, D. J.; Edwards, B. J.; Steele, W. V. A Molecular Dynamics Study of a Nafion Polyelectrolyte Membrane and the Aqueous Phase Structure for Proton Transport. J. Phys. Chem. B 2007, 111, 2208–2218. (43) Daly, K. B.; Benziger, J. B.; Panagiotopoulos, A. Z.; Debenedetti, P. G. Molecular Dynamics Simulations of Water Permeation Across Nafion Membrane Interfaces. J. Phys. Chem. B 2014, 118, 8798–8807. (44) Aaqvist, J. Ion-Water Interaction Potentials Derived from Free Energy Perturbation Simulations. J. Phys. Chem. 1990, 94, 8021–8024. (45) Abascal, J. L. F.; Vega, C. A General Purpose Model for the Condensed Phases of Water: TIP4P/2005. J. Chem. Phys. 2005, 123, 234505. (46) Sunda, A. P.; Venkatnathan, A. Molecular Dynamics Simulations of Triflic Acid and Triflate Ion/water Mixtures: A Proton Conducting Electrolytic Component in Fuel Cells. J. Comput. Chem. 2011, 32, 3319–3328. (47) Gasteiger, J.; Marsili, M. Iterative Partial Equalization of Orbital Electronegativity—a Rapid Access to Atomic Charges. Tetrahedron 1980, 36, 3219–3228. (48) Pedretti, A.; Villa, L.; Vistoli, G. VEGA: A Versatile Program to Convert, Handle and

34


Page 34 of 37

Page 35 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Visualize Molecular Structure on Windows-Based PCs. J. Mol. Graphics Modell. 2002, 21, 47–49. (49) Pedretti, A.; Villa, L.; Vistoli, G. VEGA – an Open Platform to Develop ChemoBio-Informatics Applications, Using Plug-in Architecture and Script Programming. J. Comput. Aided Mol. Des. 2004, 18, 167–173. (50) Vishnyakov, A.; Neimark, A. V. Molecular Dynamics Simulation of Nafion Oligomer Solvation in Equimolar Methanol−Water Mixture. J. Phys. Chem. B 2001, 105, 7830– 7834. (51) Okabe, T.; Kawata, M.; Okamoto, Y.; Mikami, M. Replica-Exchange Monte Carlo Method for the Isobaric–isothermal Ensemble. Chem. Phys. Lett. 2001, 335, 435–439. (52) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447. (53) Nosé, S. A Unified Formulation of the Constant Temperature Molecular Dynamics Methods. J. Chem. Phys. 1984, 81, 511. (54) Parrinello, M. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182. (55) Kofke, D. A. On the Acceptance Probability of Replica-Exchange Monte Carlo Trials. J. Chem. Phys. 2002, 117, 6911–6914. (56) Sindhikara, D. J.; Emerson, D. J.; Roitberg, A. E. Exchange Often and Properly in Replica Exchange Molecular Dynamics. J. Chem. Theory Comput. 2010, 6, 2804–2808. (57) Morrison, F. A. Understanding Rheology; Oxford University Press, 2001.

35



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(58) Daivis, P. J.; Evans, D. J. Comparison of Constant Pressure and Constant Volume Nonequilibrium Simulations of Sheared Model Decane. J. Chem. Phys. 1994, 100, 541–547. (59) Shi, Z.; Debenedetti, P. G.; Stillinger, F. H. Relaxation Processes in Liquids: Variations on a Theme by Stokes and Einstein. J. Chem. Phys. 2013, 138, 12A526. (60) Theodorou, D. N.; Boone, T. D.; Dodd, L. R.; Mansfield, K. F. Stress Tensor in Model Polymer Systems with Periodic Boundaries. Macromol. Theory Simul. 1993, 2, 191– 238. (61) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press, 2002. (62) Jones, E.; Oliphant, T.; Peterson, P.; et al., SciPy: Open Source Scientific Tools for Python; 2001; [Online; accessed 2014-08-11]. (63) Belitz, H.-D.; Grosch, W.; Schieberle, P. Food Chemistry; Springer Science & Business Media, 2004. (64) Cussler, E. L. Diffusion: Mass Transfer in Fluid Systems; Cambridge University Press, 2009. (65) Chen, S.-H.; Mallamace, F.; Mou, C.-Y.; Broccio, M.; Corsaro, C.; Faraone, A.; Liu, L. The Violation of the Stokes–Einstein Relation in Supercooled Water. PNAS 2006, 103, 12974–12978. (66) Gavach, C.; Pamboutzoglou, G.; Nedyalkov, M.; Pourcelly, G. AC Impedance InvesR Perfluorosulfonic Membranes. J. tigation of the Kinetics of Ion Transport in Nafion

Membr. Sci. 1989, 45, 37–53. (67) Chandler, D. Introduction to Modern Statistical Mechanics; Oxford University Press, 1987. 36


Page 36 of 37

Page 37 of 37

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(68) Connolly, M. L. Analytical Molecular Surface Calculation. J. Appl. Crystallogr. 1983, 16, 548–558. (69) Elliott, J. A.; Hanna, S.; Elliott, A. M. S.; Cooley, G. E. Interpretation of the SmallAngle X-Ray Scattering from Swollen and Oriented Perfluorinated Ionomer Membranes. Macromolecules 2000, 33, 4161–4171. (70) Tricoli, V. Proton and Methanol Transport in Poly(perfluorosulfonate) Membranes Containing Cs+ and H+ Cations. J. Electrochem. Soc. 1998, 145, 3798–3801.

Figure 11: Table of Contents Image

37


Viscosity of Nafion Oligomers as a Function of Hydration and

Recommend Documents